###### Abstract

The thermal restoration of chiral symmetry in QCD is known to proceed by an
analytic crossover, which is widely expected to turn into a phase transition
with a critical endpoint as the baryon density is increased. In the absence of
a genuine solution to the sign problem of lattice QCD, simulations at zero and
imaginary baryon chemical potential in a parameter space enlarged by a
variable number of quark flavours and quark masses constitute a viable way to
constrain the location of a possible non-analytic phase transition and its
critical endpoint. In this article I review recent progress towards an
understanding of the nature of the transition in the massless limit, and its
critical temperature at zero density. Combined with increasingly detailed
studies of the physical crossover region, current data bound a possible
critical point to $\mu_{B}$ $>$$\sim$ $3T$.

###### keywords:

QCD phase diagram; chiral symmetry restoration; lattice QCD

1 1 0 Academic Editor(s): Angel Gómez Nicola and Jacobo Ruiz de Elvira 22
September 2021 26 October 2021 https://doi.org/ Lattice Constraints on the QCD
Chiral Phase Transition at Finite Temperature and Baryon Density Lattice
Constraints on the QCD Chiral Phase Transition at Finite Temperature and
Baryon Density Owe Philipsen ${}^{\orcidA{}}$ Owe Philipsen Philipsen, O.

##  1 Introduction

Many salient features of the spectrum of hadrons and their interactions
observed in nature are determined by the near-chiral symmetry of the QCD
Lagrangian and its spontaneous breaking by the QCD vacuum. The fact that
nature breaks chiral symmetry also explicitly by non-vanishing quark masses
can be taken into account in terms of chiral perturbation theory Gasser and
Leutwyler (1984); Leutwyler (1994), due to the smallness of the $u$\- and
$d$-quark masses.

In a hot and/or dense medium, chiral symmetry is expected to be gradually
restored once the temperature exceeds $T$ $>$$\sim$ $160$ MeV or the baryon
chemical potential $\mu_{B}$ $>$$\sim$ $1$ GeV. Associated with this change of
the realised symmetry, one also expects a change of dynamics and its
underlying degrees of freedom which, asymptotically, should become quarks and
gluons. Similar to the vacuum properties of the theory, one also expects the
properties of this thermal transition to be closely related to those of the
theory in the massless limit.

The low energy scales of the transition region demand a non-perturbative first
principles approach like lattice QCD. Unfortunately, a severe sign problem
prohibits simulations by importance sampling for non-vanishing chemical
potential, for introductions see Gattringer and Lang (2010); Philipsen (2010);
Aarts (2016). Despite tremendous efforts over several decades, no genuine
solution to this problem is available to date, and knowledge of the QCD phase
diagram by direct calculation remains scarce. Nevertheless, during the last
years there has been considerable progress in studying the properties of the
thermal QCD transition in an enlarged parameter space, where the number of
flavours and quark masses is varied away from their physical values. Besides
giving theoretical insight into the interplay of symmetries and dynamics in
controlled situations, such studies also provide constraints on the physical
phase diagram, which are beginning to become phenomenologically relevant. In
this article, I collect some recent results allowing to better understand and
remove lattice artefacts, resulting in a consistent picture that also
accommodates apparently conflicting statements based on earlier simulations.
In particular, I will focus on developments suggesting a modified version of
the Columbia plot and constraints on the location of a possible critical
point.

##  2 Some Lattice Essentials

For readers less familiar with lattice QCD, it may be useful to recall a
number of basic definitions and issues that are important for the
interpretation of the following lattice results. For details, see the
corresponding textbooks Gattringer and Lang (2010); Montvay and Munster
(1997); DeGrand and Detar (2006). Conventionally one considers Euclidean
space-time discretised by hypercubic lattices with $N_{s}^{3}\times N_{\tau}$
points separated by a lattice spacing $a$. This represents a box with spatial
length $L=N_{s}a$ and Euclidean time extent $aN_{\tau}=1/T$, which corresponds
to inverse temperature as in continuum thermal field theory, once
(anti-)periodic boundary conditions for (fermionic) bosonic fields have been
applied. When a quantum field theory is formulated on the lattice, its action
(and all observables) differ from those in the continuum. As an example,
consider Wilson’s $SU(N_{c})$ pure gauge action

$S_{YM}=\beta\sum_{p}\Big{[}1-\frac{1}{N_{c}}{\rm
Tr}U(p)\Big{]}\;,\quad\beta=\frac{2N_{c}}{g^{2}}\;,$ (1)

which is constructed from so-called plaquette variables $U(p)$. These are the
fundamental squares of the lattice consisting of four gauge links
$U_{\mu}(x)=\exp(-igaA_{\mu}(x))$, which are group elements of $SU(N_{c})$.
For sufficiently small lattice spacing the plaquette can be expanded,

$U(p)=1+ia^{2}gF_{\mu\nu}-\frac{a^{4}g^{2}}{2}F_{\mu\nu}F_{\mu\nu}+O(a^{6}),$
(2)

and the classical action (1) is seen to reproduce the continuum Yang-Mills
action in the limit $a\rightarrow 0$. However: for any finite lattice spacing
$a$ the action differs by terms of $O(a^{2})$ from the continuum action. Such
terms disappearing in the continuum are called “irrelevant”, at finite lattice
spacing they constitute lattice artefacts. An improvement of lattice actions
in the sense of reducing such artefacts Symanzik (1983) can be achieved by
adding irrelevant terms to the lattice action, which subtract those already
present. Quite generally, lattice actions are then not unique but may differ
by arbitrary irrelevant terms. All such differences between valid actions and
observables must disappear in the continuum limit.

Approaching the continuum limit is a highly non-trivial affair, as one has to
take $a\rightarrow 0$ while keeping other dimensionful quantities, like masses
or temperature, fixed along a so-called line of constant physics in the
lattice parameter space. The lattice spacing is tuned indirectly through the
running gauge coupling evaluated at the cutoff scale, $g=g(\Lambda\sim
a^{-1})$, so that by asymptotic freedom $g(a\rightarrow 0)\rightarrow 0$. This
implies that the lattice gauge coupling $\beta(a\rightarrow
0)\rightarrow\infty$. For thermodynamics, keeping temperature constant during
the continuum approach implies $N_{\tau}\rightarrow\infty$ in the continuum
limit. Thus, while the precise value of the lattice spacing is determined by
some form of scale setting, finer thermodynamics lattices imply larger
$N_{\tau}$, and lattice artefacts can be parametrised in units of temperature
as powers of the dimensionless $aT=N_{\tau}^{-1}$.

The most problematic lattice artefacts for QCD by far are those afflicting the
fermion formulation. In particular, it is impossible to put chiral fermions on
a lattice without artificial doubler degrees of freedom while maintaining
locality Nielsen and Ninomiya (1981). It is clear that this implies severe
caveats for any investigation of chiral symmetry breaking. In a nut shell, the
current choices are:

  * •

The Wilson formulation: By adding irrelevant mass terms $\sim a^{-1}$, the
doubler degrees of freedom become heavy and decouple in the continuum limit.
However, for any finite lattice spacing chiral symmetry is broken completely
by these terms.

  * •

The staggered formulation: Spin and flavour degrees of freedom are distributed
to different lattice sites, which can be done so as to reduce the number of
doublers. The remaining ones are removed by taking an appropriate root of the
fermion determinant, which can only be fully valid in the continuum. In this
formulation the original chiral symmetry is reduced from $SU(2)_{L}\times
SU(2)_{R}\rightarrow O(2)$.

  * •

Formulations with a lattice version of full chiral symmetry exist in the form
of domain wall or overlap fermions, and are expected to eventually supersede
the previous formulations. However, they require complicated non-local
constructions and currently are more expensive to simulate by over an order of
magnitude.

Due to the enormous numerical cost of thermodynamical investigations for light
fermions, most lattice results covered here are based on different versions of
the staggered or Wilson discretisations. In order to avoid the damage these do
to chiral symmetry and its spontaneous breaking, it is mandatory to first take
the continuum limit to remove the lattice artefacts, and only then approach
the chiral limit.

A general difficulty in investigating phase diagrams is the fact that non-
analytic phase transitions only exist in infinite volume Yang and Lee (1952);
Lee and Yang (1952), while numerical simulations are obviously limited to
finite boxes. A change of dynamics associated with a phase transition is
signalled by a rapid change of the expectation values of suitable observables
$O({\bf x})$, accompanied by maximal fluctuations, i.e., peaks in their
susceptibilities,

$\chi_{O}=\int d^{3}x\;\Big{(}\langle O({\bf x})O(\mathbf{0})\rangle-\langle
O({\bf x})\rangle\langle O(\mathbf{0})\rangle\Big{)}\;.$ (3)

This allows to determine the location of a transition as a function of the
bare parameters of the theory. In a finite box this is always an analytic
crossover, and the behaviour on different volumes has to be compared to see
whether a non-analytic phase transition builds up in the thermodynamic limit.
The peak height will diverge at a second-order transition with $\sim
V^{\sigma}$, where $\sigma$ is a combination of critical exponents appropriate
for the observable and the universality class of the transition. For a first
order transition the divergence is $\sim V$, whereas for a crossover the peak
height will saturate at a finite value. More intricate (and reliable) finite
size scaling studies are based on the first three generalised cumulants
($n=2,3,4$ ) of the fluctuation of an observable, $\delta O=O-\langle
O\rangle$,

$B_{n}(\delta O)=\frac{\langle(\delta O)^{n}\rangle}{\langle(\delta
O)^{2}\rangle^{n/2}}\;,$ (4)

which besides the peak height also quantify the skewness and the tails of the
statistical distribution of the observable, thus giving access to various
critical exponents.

Finally, the infamous sign problem at finite baryon density appears because of
a property of the fermion determinant, defined with complex chemical potential
for baryon number (in continuum notation),

$\displaystyle\det(\cancel{D}+m-\mu\gamma_{0})$ $\displaystyle=$
$\displaystyle{\det}^{*}(\cancel{D}+m+\mu^{*}\gamma_{0})\;.$ (5)

Thus, the fermion determinant is complex for real chemical potential, but
stays real for purely imaginary chemical potential. For real chemical
potential, all imaginary parts cancel out of the partition function exactly,
but the real part is oscillatory and averages to the so-called average sign,
evaluated by the partition function with $\mu=0$,

$Z(\mu)=Z(0)\left\langle\frac{\det(\mu)}{\det(0)}\right\rangle_{0}=Z(0)\exp\Big{(}-V(f(\mu)-f(0))/T\Big{)}\;.$
(6)

This is also the simplest reweighting factor multiplying any observable, when
a $\mu=0$ ensemble is reweighted to $\mu\neq 0$. Its size is governed by the
free energy difference of the systems with and without chemical potential, and
is exponentially suppressed for large volumes. Hence all expectation values of
observables are lost in the statistical errors sooner or later. It must be
stressed that this is only a problem for numerical evaluations by importance
sampling, other evaluations and the physics itself have no issue with chemical
potential. For more details, see introductions like Gattringer and Lang
(2010); Philipsen (2010); Aarts (2016). So far there is no genuine algorithmic
cure for the sign problem, for an overview of earlier attempts see de Forcrand
(2009) and the proceedings of the annual lattice conferences for updates.

##  3 The Columbia Plot at Zero Baryon Density

As a point of departure, consider the thermal QCD transition at zero baryon
density, $\mu_{B}=0$. The nature of this transition for physical quark masses
has been known for some time to be an analytic crossover Aoki et al. (2006).
Away from the physical point, the order of the QCD thermal transition with
$N_{f}=2+1$ quarks as a function of quark masses is usually specified in a so-
called Columbia plot Brown et al. (1990), for which two possible realisations,
to be discussed below, are shown in Figure 1.

###  3.1 The Deconfinement Transition

Even though it is far from the physical point or the chiral limit, the heavy
mass corner is interesting in its own right and promises insights, which will
also be useful for light quarks. In the quenched limit QCD reduces to $SU(3)$
Yang-Mills theory in the presence of static quarks, whose propagators are
Polyakov loops which represent an order parameter for the global $Z(3)$ center
symmetry related to confinement McLerran and Svetitsky (1981). Its spontaneous
breaking at some critical temperature proceeds by a first-order deconfinement
phase transition Boyd et al. (1996), whose latent heat has recently also been
computed Shirogane et al. (2016). In the presence of dynamical quarks, the
center symmetry is explicitly broken by $1/m_{q}$ and the first-order phase
transition weakens with diminishing quark masses, until it disappears along a
second-order critical line in the 3d $Z(2)$ universality class. Contrary to
the chiral limit, this entire region can be simulated directly and a complete
non-perturbative understanding should be attainable in the near future. The
location of the critical line is under investigation Saito et al. (2011);
Ejiri et al. (2020); Cuteri et al. (2021); Kiyohara et al. (2021). While
continuum extrapolations are not yet available, the presently finest lattices
predict the pseudo-scalar mass evaluated on the critical point for $N_{f}=2$
to be about $m^{c}_{PS}\sim 4$ GeV Cuteri et al. (2021). Together with the
latent heat, these quantities should permit a detailed understanding of the
confinement-deconfinement transition and its interplay with screening by
dynamical quarks, which will be valuable for the physics interpretation of the
crossover region at the physical point. Finally, these quantities can serve as
benchmarks to test and tune truncations involved in functional renormalisation
group or Dyson-Schwinger approaches Fischer et al. (2015), or effective
lattice theories for the heavy quark region, which are employed to describe
finite density physics Fromm et al. (2012, 2013).

Figure 1: Possible scenarios for the nature of the thermal QCD transition as a
function of the quark masses, with a first-order (left) or a second-order
(right) transition in the chiral limit of the light quarks. Every point of the
plot represents a phase boundary, with an implicitly associated
(pseudo-)critical temperature $T_{c}(m_{u,d},m_{s})$. (Here and in the
following, “phase boundary” refers to a (pseudo-)critical combination of
parameters irrespective of the nature of the transition, which can be first
order, second order or crossover.)

###  3.2 The Chiral Transition at Zero Baryon Density

The situation is far more difficult in the opposite limit of massless quarks,
because it cannot be simulated directly. For several decades expectations have
been based on an analysis of three-dimensional sigma models, augmented by a ’t
Hooft term for the $U(1)_{A}$ anomaly, which represent Landau-Ginzburg-Wilson
effective theories for the chiral condensate as the order parameter of the
transition. The renormalisation group flow based on the epsilon expansion
Pisarski and Wilczek (1984) predicts the chiral phase transition to be first-
order for $N_{f}\geq 3$, whereas the case of $N_{f}=2$ is found to crucially
depend on the fate of the anomalous $U(1)_{A}$ symmetry: If the latter remains
broken at $T_{c}$, the chiral transition should be second order in the
$O(4)$-universality class, whereas its effective restoration would enlarge the
chiral symmetry and push the transition to first-order. For non-zero quark
masses, chiral symmetry is explicitly broken. A first-order chiral phase
transition then weakens to disappear at a $Z(2)$ second-order critical
boundary, just like the deconfinement transition in the heavy mass regime,
whereas a second-order transition disappears immediately. This results in the
two scenarios for the Columbia plot depicted in Figure 1. The results of the
epsilon expansion were essentially confirmed by numerical simulations of
three-dimensional sigma models Gausterer and Sanielevici (1988) and a
perturbative expansion in fixed dimensions Butti et al. (2003). A later high-
order perturbative analysis for the case with an effectively restored
$U(1)_{A}$ at the transition temperature finds a possibility for second-order
transitions also in this case, but with a symmetry breaking pattern
$U(2)_{L}\otimes U(2)_{R}\rightarrow U(2)_{V}$ signalling a different
universality class Pelissetto and Vicari (2013). Finally, a recent functional
renormalisation group analysis applied to $U(1)_{A}$ restoration in QCD, which
is in good agreement with lattice susceptibilities at finite quark masses
Braun et al. (2020a), favours the $O(4)$ scenario in the chiral limit Braun et
al. (2020b).

A large number of lattice simulations has been devoted to disentangle which of
these situations is realised in QCD. When the Columbia plot is considered on
the lattice with different lattice spacings, its parameter space is enlarged
by an additional axis. The chiral critical line then traces out a chiral
critical surface whose shape is discretisation-dependent, while for all valid
discretisations it must of course emerge from the same critical line in the
continuum limit. The following subsections collect the current evidence on the
location of the chiral boundary line in the Columbia plot.

###  3.3 $N_{\text{f}}=2$ and $N_{\text{f}}=2+1$

Early simulations (without sophisticated finite size scaling analyses) on
coarse $N_{\tau}=4$ lattices were fully consistent with the expected scenarios
shown in Figure 1: A first-order region could be clearly seen for
$N_{\text{f}}=3$, whereas the smallest available masses were consistent with a
continuous transition or crossover for $N_{\text{f}}=2$, both for unimproved
staggered Brown et al. (1990) and unimproved Wilson Iwasaki et al. (1996)
fermions. More recently and using finite size scaling of cumulants, a narrow
first-order region was also identified for $N_{f}=2$, again for unimproved
staggered Bonati et al. (2014); Cuteri et al. (2018) and unimproved Wilson
Philipsen and Pinke (2016) fermions on $N_{\tau}=4$. However, the location of
the $Z(2)$-boundary varies widely between these, indicating large cutoff
effects.

Recent investigations using the Highly Improved Staggered Quark (HISQ) action
start at the $N_{f}=2+1$ physical point and then gradually reduce the light
quark mass value, until they correspond to $m_{PS}\approx 45$ MeV, on lattices
$N_{\tau}=6,8,10$ Ding et al. (2019); Kaczmarek et al. (2020). No sign of a
first-order transition is detected. Thus either a $Z(2)$-critical point
bounding a very narrow first-order region is approached, or a second-order
transition in the chiral limit. The analysis employs a renormalisation group
invariant combination of chiral condensates as order parameter representing
the magnetisation-like variable, and the light quark masses $m_{l}=m_{u,d}$ in
units of the strange quark mass as symmetry breaking field,

$M=\Delta_{ls}=2(m_{s}\langle\bar{\psi}\psi\rangle_{l}-m_{l}\langle\bar{\psi}\psi\rangle_{s})/f_{K}^{4},\quad
H=m_{l}/m_{s}\;.$ (7)

Near a critical point the magnetic observables are dominated by universal
scaling functions

$M(t,h)=h^{1/\delta}f_{G}(z)+\ldots\;,\quad\chi_{M}(t,h)=\frac{\partial
M}{\partial H}=h_{0}^{-1}h^{1/\delta-1}f_{\chi}(z)+\ldots\;,$ (8)

with a scaling variable $z=t/h^{1/\beta\delta}$ expressed in terms of the
reduced temperature and external field,
$t=t_{0}(T-T^{0}_{c})/T^{0}_{c},h=H/h_{0}$, which contain the unknown critical
temperature in the chiral limit, $T_{c}^{0}$, and two non-universal parameters
$t_{0},h_{0}$.

Figure 2: (Left): Chiral susceptibility at the physical strange quark mass for
a range of decreasing light quark masses on $N_{\tau}=8$ lattices with HISQ
action, from Ding et al. (2019). (Right): Magnetic equation of state at the
pseudo-critical temperature approaching the chiral limit. Lines represent fits
to $O(N)$\- ($N=2,4$ second-order scenario) or $Z(2)$-scaling with a finite
critical quark mass (first-order scenario), from Kaczmarek et al. (2020).

Figure 2 (left) shows the chiral susceptibility as a function of temperature
for a series of decreasing light quark masses on $N_{\tau}=8$. A considerable
reduction of the pseudo-critical temperature, defined by the peak location of
the susceptibility, is apparent as the light quark mass is reduced. Moreover,
the peak height is growing as expected when approaching a true phase
transition. To test for the universality class of the approached transition,
the ratio $M/\chi_{M}$ is shown in Figure 2 (right), which near a critical
point is dominated by universal scaling functions. While it is difficult to
distinguish between $O(4)$ and $Z(2)$, it is apparent that any finite critical
quark mass bounding a first-order region has to be excessively small to be
consistent with the observed behaviour.

There is a similar scaling expression for the approach of the pseudo-critical
crossover temperature to the critical temperature in the chiral limit,

$\displaystyle T_{pc}(H)$ $\displaystyle=$ $\displaystyle
T_{c}(m_{l}=0)\Big{(}1+\frac{z_{X}}{z_{0}}H^{1/\beta\delta}\Big{)}+\mbox{sub-
leading}\;,$ (9) $\displaystyle T_{c}(m_{l}=0)$ $\displaystyle=$
$\displaystyle 132^{+3}_{-6}\;\mathrm{MeV}\;.$ (10)

In Ding et al. (2019) the variation between the possible sets of critical
exponents is observed to be very small, so that an extrapolation makes sense
even without knowledge of the true universality class. Moreover, the
extrapolations were checked to be stable under an exchange of the order of the
continuum and chiral extrapolations, leading to the critical temperature as a
first result for the chiral phase transition in the massless limit.

Simulations with Wilson fermions do not yet reach such small quark masses, but
are fully consistent with this picture down to the physical pion mass. In
particular, extrapolations of the pseudo-critical temperature to the chiral
limit have also been attempted by an investigation using $N_{\text{f}}=2+1+1$
$O(a)$-improved twisted mass Wilson fermions, with strange and charm quark
masses close to their physical values Kotov et al. (2021). In these
simulations lattice spacings are held fixed in a range
$a/\mathrm{fm}\in[0.062,0.082]$, and temperature is varied by
$N_{\tau}\in[12,20]$. The pseudo-critical temperature is defined in three
different ways: By the maximum of the chiral susceptibility $\chi$ as above,
and as the inflection point of polynomial fits to the chiral condensate
($\Delta$) and a subtracted version without additive renormalisation
$(\Delta_{3})$. The resulting temperatures as a function of pseudo-scalar
(pion) mass are shown in Figure 3 (left), together with extrapolations to the
chiral limit, employing either $O(4)$ exponents or $Z(2)$-exponents and a
critical pseudo-scalar mass up to $m_{\pi}\sim 100$ MeV. Again, it is not
possible to distinguish between these scenarios. As in the previous case, the
extrapolated critical temperature in the chiral limit is therefore robust
under changes of the critical exponents and quoted as

$T_{c}^{0}=134^{+6}_{-4}\;\mathrm{MeV},$ (11)

in remarkable agreement with the staggered result.

Figure 3 (right) shows an investigation of sections of the chiral critical
line using $O(a)$ clover-improved Wilson fermions Nakamura et al. (2019).
Starting point are the data for $N_{f}=3$ to be discussed separately below,
and on $N_{\tau}=6$ further points at larger strange quark masses have been
added. The critical line is then fitted assuming a tricritical strange quark
mass as explained in Section 3.5 plus polynomial corrections. Note that this
discretisation features a much wider first-order region, which even contains
the physical point on the coarser lattices. This must be a lattice artefact,
and the first-order region rapidly shrinks as $N_{\tau}$ is increased.

Figure 3: (Left): Pseudo-critical temperature of the crossover defined by the
chiral susceptibility $\chi$, the inflection point of the chiral condensate
$\Delta$ or an additively renormalised chiral condensate $\Delta_{3}$, for
$N_{\text{f}}=2+1+1$ twisted mass Wilson fermions close to the continuum.
Lines represent chiral extrapolations according to the $O(4)$ second-order or
finite critical $Z(2)$-mass scenario. From Kotov et al. (2021). (Right):
Columbia plot expressed in $\eta,\pi$-masses in units of the Wilson flow
parameter $t_{0}$. Critical points have been determined using an
$O(a)$-improved Wilson action. The first-order region includes the physical
point on coarse lattices, but shrinks drastically as $N_{\tau}$ is increased.
From Nakamura et al. (2019).

Several conclusions can be drawn from these results. Firstly, the width of a
potential first-order region as in Figure 1 (left) is bounded to a small
fraction of the physical light quark (or pion) masses. Second, the numerical
proximity of the critical exponent combinations $1/(\beta\delta)$ for the 3D
$O(2),O(4)$ and $Z(2)$ universality classes appears to allow for a robust
extrapolation of the chiral transition temperature to the massless limit with
remarkably small uncertainties. Conversely this statement implies, however,
that it is impossible to firmly identify the universality class in this way,
which would require exponentially accurate data. This problem might be avoided
by looking at the scaling of energy-like variables, which are governed by the
critical exponent $\alpha$ that changes sign between the $O(2),O(4)$ and the
$Z(2)$ universality classes. It was shown that the Polyakov loop behaves as an
energy-like observable, but unfortunately a firm distinction between
universality classes would require a further substantial reduction of the
light quark mass Clarke et al. (2021). Finally, note that the value of
$T_{c}(m_{l}=0)$ is $\sim 25$ MeV lower than the pseudo-critical temperature
at the physical point. This constrains the possible location of a chiral
critical point at finite baryon density, as will be discussed in Section 5.1.

###  3.4 $N_{\text{f}}=3$ and $N_{\text{f}}=4$

The $N_{\text{f}}=3$ theory, represented by the diagonal of the Columbia plot,
Figure 1, is particularly interesting because a sufficiently large first-order
transition and its associated $Z(2)$-critical endpoint are accessible by
direct simulations, rather than by extrapolation. However, the situation has
become more complicated since the early results, as shown by the compilation
of critical pseudo-scalar masses obtained by different actions and lattice
spacings in Table 3.4. The largest value in the table differs by almost an
order of magnitude from the lowest bound! Unless any of the employed lattice
actions is fundamentally flawed, all of them must eventually converge to the
same continuum limit, which is not in sight yet. However, the general trend is
for the critical mass to shrink when either a lattice with fixed action is
made finer (larger $N_{\tau}$), or when improved actions are employed. This
points to enormous cutoff effects, which quite generally appear to increase
the first-order region, i.e., make the transition stronger. Another
observation is for the first-order region to be considerably larger in the
case of Wilson fermions, which is consistent with stronger cutoff effects due
to the complete violation of chiral symmetry.

[H] Summary of previous studies (continued from Varnhorst (2015); de Forcrand
et al. (2007)) of the $N_{f}=3$ chiral critical point at $\mu_{B}=0$ . Action
$N_{\tau}$ $m_{PS}^{c}$ Ref. Year unimproved staggered 4 $\sim 290$ MeV Karsch
et al. (2001) 2001 p4 staggered 4 $\sim 67$ MeV Karsch et al. (2004) 2004
unimproved staggered 6 $\sim 150$ MeV de Forcrand et al. (2007) 2007 HISQ
staggered 6 $\lesssim 45$ MeV Bazavov et al. (2017) 2017 stout staggered 4-6
$\sim 0$? Varnhorst (2015) 2014 Wilson-$O(a)$-impr. 6-8 $\sim 300$ MeV Jin et
al. (2015) 2014 Wilson-$O(a)$-impr. 4-10 $\lesssim 170$ MeV Jin et al. (2017)
2017 Wilson-$O(a)$-impr. 4-12 $\lesssim 110$ MeV Kuramashi et al. (2020) 2020

As an example, consider Figure 4 ( left), where the critical pseudo-scalar
mass is shown from a comprehensive long term study with $O(a)$-improved Wilson
fermions for different lattice spacings. The critical meson mass shrinks by
nearly a factor of two when the lattice is changed between $N_{\tau}=6$ and
$N_{\tau}=12$. Its continuum value depends strongly on the extrapolation
function and could even be zero. The same effect is observed for
$N_{\text{f}}=4$ unimproved staggered fermions de Forcrand and D’Elia (2017),
as shown in Figure 4 (right). One expects the transition in this case to be
more strongly first-order than for $N_{\text{f}}=3$, since the spontaneously
breaking symmetry is larger. This is borne out both for Wilson Ohno et al.
(2018) and staggered fermions. Note that for $N_{\text{f}}=4$ staggered
fermions no rooting is applied, which therefore cannot cause any unphysical
behaviour. Moreover, the effect of reducing the lattice spacing is similar for
$N_{\text{f}}=2,3,4$ staggered fermions, i.e., with and without rooting.

We can then conclude that Wilson and staggered discretisations show the same
qualitative behaviour, with a monotonic decrease of the critical pseudo-scalar
mass $m_{PS}^{c}(N_{\tau})$ as the continuum is approached by increasing
$N_{\tau}$. Quantitative contradictions between different actions are only
avoided, if the final answer $m_{PS}^{c}(\infty)$ is below the lowest
available value $m_{PS}(N_{\tau})$. In summary, since finer lattices and
improved actions have become available, the situation for $N_{\text{f}}=3,4$
is quite similarly ambiguous to the previously discussed $N_{\text{f}}=2$ and
$N_{\text{f}}=2+1$.

Figure 4: Critical pseudo-scalar mass in units of temperature as a function of
$N_{\tau}^{-2}=(aT)^{2}$.
(Left): $N_{f}=3$ $O(a)$-improved Wilson fermions with different continuum
extrapolations, from Kuramashi et al. (2020). (Right): $N_{f}=3,4$ unimproved
staggered fermions, from de Forcrand and D’Elia (2017).

2

###  3.5 Three-State Coexistence and Tricritical Scaling

The uncertainty about the size of a potential first-order region can be
removed by a different analysis, which leads to a consistent picture covering
the results of all discretisations reported here. Starting point is the
observation that a change of the chiral phase transition in the massless limit
from first order to second order, brought about by variation of a continuous
parameter, necessarily proceeds by a tricritical point Rajagopal (1995). To
see this, assume the scenario in Figure 1 (right), and consider the
corresponding phase diagram for a fixed $m_{s}<m_{s}^{\mathrm{tric}}$, Figure
5 (left). For non-vanishing quark masses, the first-order transition weakens
to disappear in a critical endpoint, and the phase diagram is symmetric under
a chiral rotation $m_{u,d}\rightarrow-m_{u,d}$. Hence, in the massless limit
the critical temperature marks a triple point of three-state coexistence: Two
states with broken chiral symmetry, $\pm\langle\bar{\psi}\psi\rangle\neq 0$,
and one with chiral symmetry restored, $\langle\bar{\psi}\psi\rangle=0$. When
$m_{s}$ is increased, these points trace out a triple line (corresponding to
the green $m_{u,d}=0$ line in Figure 1) along which the latent heat of the
first-order transition decreases, until it vanishes in a tricritical point. At
the same time, the two wings of finite mass first-order transitions shrink, so
that the tricritical point marks the confluence of two critical endpoints.

This opens an additional path of investigation. Quite generally, tricritical
points come with their own set of critical exponents, which take mean field
values since their upper critical dimension is three. For a general
introduction to the theory of tricritical scaling, see Lawrie and Sarbach
(1984). In particular, the functional form of the critical wing lines entering
the tricritical point as a function of the symmetry breaking field, which is
the bare light quark mass in our case, is governed by a critical exponent,

$m_{s}^{c}(m_{l})=m_{s}^{\mathrm{tric}}+\mathcal{A}_{1}\cdot
m_{l}^{2/5}+\mathcal{O}\big{(}m_{l}^{4/5}\big{)}\;.$ (12)

This allows for an extraction of the location of a tricritical point with
fixed, known exponents, which is much easier than having to determine the
exponents themselves. Knowing the location of the tricritical point then fixes
the order of the transition on both sides of it, albeit without information
about the universality class. Attempts in this direction were made in de
Forcrand and Philipsen (2007); Nakamura et al. (2019), but they are not yet
conclusive because the lattices are coarse, whereas very high precision data
very close to the chiral limit would be required on sufficiently fine lattices
in order to unambiguously distinguish tricritical scaling from an ordinary
polynomial. However, this problem can be circumvented by a change of variables
in the bare parameter space of the theory.

Figure 5: (Left): Phase diagram with a first-order chiral phase transition.
For $T<T_{c}$ the $T$-axis corresponds to a coexistence line of
$\pm\langle\bar{\psi}\psi\rangle\neq 0$, and $T_{c}(m_{u,d}=0)$ represents a
triple point. (Right): Columbia plot for mass-degenerate quarks. Again every
point represents a phase boundary and has an implicitly associated
$T_{c}(m,N_{\text{f}})$.

###  3.6 The Chiral Phase Transition as a Function of $N_{\text{f}}$ and
$N_{\tau}$

In Cuteri et al. (2018) it was suggested to consider a different version of
the Columbia plot with degenerate quarks only, in which the variable strange
quark mass $m_{s}$ is traded for a variable $N_{\text{f}}$ analytically
continued to non-integer values. Starting with integer $N_{\text{f}}$, after
Grassmann integration the QCD partition function reads

$Z(N_{\text{f}},g,m)=\int{\cal D}A_{\mu}\;(\det
M[A_{\mu},m])^{N_{\text{f}}}\;e^{-\mathcal{S}_{\mathrm{YM}}[A_{\mu}]}\;,$ (13)

which can now be formally viewed as a statistical system of gauge field
variables depending on a continuous parameter $N_{\text{f}}$. In the lattice
formulation with rooted staggered fermions, whose determinant is raised to the
power $N_{\text{f}}/4$ in order to describe $N_{\text{f}}$ mass-degenerate
quarks, this is implemented straightforwardly. The Columbia plot Figure 1
(right) then translates to the version shown in Figure 5 (right), where the
tricritical strange quark mass is replaced by a tricritical number of
flavours, $2<N_{\text{f}}^{\mathrm{tric}}<3$, and the $N_{\text{f}}$-axis to
the right of it corresponds to the new triple line. The crucial advantage in
this modified parameter space is that, since there is no chiral transition for
$N_{\text{f}}=1$, a tricritical point $N_{\text{f}}^{\mathrm{tric}}>1$ is
guaranteed to exist as soon as there is a first-order region for any
$N_{\text{f}}>1$. At finite lattice spacing there is an additional parameter,
and the $Z(2)$-critical boundary traces out a chiral critical surface which
terminates in a tricritical line in the chiral limit of the lattice theory.
Since there is ample evidence for first-order transitions at coarse lattice
spacings, the task is now reduced to track their boundaries in the lattice
bare parameter space and determine the location of the tricritical line in the
lattice chiral limit.

Numerical results obtained with unimproved staggered fermions for a wide range
of flavours and $N_{\tau}=4,6,8$ Cuteri et al. (2018, 2021) are shown in
Figure 6. The left plot represents the lattice version of the Columbia plot
Figure 5 (right). One observes a summary of the previous discussion, namely a
first-order transition getting stronger with $N_{\text{f}}$ and with
increasing lattice spacing. The cutoff effects are growing with
$N_{\text{f}}$, for $N_{\text{f}}=7$ the critical bare quark mass is reduced
by a stunning factor of five when the lattice spacing is reduced ($N_{\tau}$
is increased) by a factor of two! No continuum critical mass $m_{c}$ is
discernible yet for any $N_{\text{f}}$. On the other hand, the intercepts of
the critical lines at $m=0$ constitute the sought-after tricritical line
$N_{\text{f}}^{\mathrm{tric}}(N_{\tau})$. Tricritical scaling can be
appreciated in the rescaled Figure 6 (right), where the data approach a
leading order scaling relation Cuteri et al. (2018). This allows for an
extrapolation

$\displaystyle am_{c}\big{(}N_{\text{f}}(N_{\tau}),N_{\tau}\big{)}$
$\displaystyle=\mathcal{D}_{1}(N_{\tau})\big{(}N_{\text{f}}-N_{\text{f}}^{\mathrm{tric}}(N_{\tau})\big{)}^{5/2}+\ldots\;.$
(14)

Figure 6: The chiral critical surface for unimproved staggered fermions,
projected onto the $(m/T,N_{\text{f}})$ plane Cuteri et al. (2021) (left) and
onto the $(am,N_{\text{f}})$ plane for $N_{\tau}=4$ Cuteri et al. (2018)
(right). Every point of the plane represents a phase boundary with an
implicitly tuned $\beta(am,N_{\text{f}},N_{\tau})$. Regions above the
$Z(2)$-lines represent crossover, those below first-order transitions. Fits on
the left show next-to-leading order tricritical scaling, on the right leading-
order tricritical scaling.

The scaling relation has been inverted here, because in simulations
$N_{\text{f}}$ is fixed whereas the critical quark mass is computed with
statistical errors. Note also that the $N_{\text{f}}=2$ data point in Figure 6
(right) has been obtained by a tricritical extrapolation in imaginary chemical
potential at fixed $N_{\text{f}}=2$ Bonati et al. (2014). This is an
independent confirmation of the bare quark mass as a tricritical scaling field
near its chiral limit. Scaling fits in the left figure take into account the
next-to-leading order terms to predict
$N_{\text{f}}^{\mathrm{tric}}(N_{\tau}=4)\approx 1.71(3)$ and
$N_{\text{f}}^{\mathrm{tric}}(N_{\tau}=6)\approx 2.20(8)$. While no
extrapolation for $N_{\tau}=8$ is available yet, the critical line further
flattens and appears to shift to the right, implying a line
$N_{\text{f}}^{\mathrm{tric}}(N_{\tau})$. One concludes that, for unimproved
staggered fermions, the $N_{\text{f}}=2$ massless theory shows a first-order
transition on $N_{\tau}=4$, but a second-order transition on all finer
lattices and in the continuum. The question now is what happens to the
$N_{\text{f}}\geq 3$ theories.

This can be seen in a different variable pairing in Figure 7 (left), where the
same critical bare quark masses are rescaled by the tricritical exponent, and
plotted as a function of $N_{\tau}$ for different fixed $N_{\text{f}}$-values.
Only a slight curvature is exhibited by those $N_{\text{f}}$-lines with three
data points, which thus are compatible with next-to-leading order scaling and
a tricritical point at some finite $N_{\tau}^{\mathrm{tric}}(N_{\text{f}})$,
obtained by the corresponding extrapolations.

###  3.7 Tricritical Scaling for Wilson Fermions

What about Wilson fermions, which exhibit the widest first-order region among
the discretisation schemes discussed here? In this case the conceptual
situation is more difficult, because chiral symmetry is broken completely so
that the bare quark mass and chiral condensate receive additive
renormalisation as well as the ordinary multiplicative one. Moreover, the bare
parameter phase diagram is expected to show lattice artefacts such as parity
broken phases Aoki (1984, 1986) or first-order bulk transitions Sharpe and
Singleton (1998); Farchioni et al. (2005); Burger et al. (2013), whose
location depends on the details of the lattice action and $N_{\tau}$. This may
considerably complicate the bare parameter phase diagram while, of course, all
actions should eventually agree in the continuum limit.

Here we consider an RG-improved Wilson gauge and non-perturbatively
O(a)-improved Wilson fermion action, for which the Aoki phase is limited to
the strong coupling regime and no other unphysical phases have been observed
Ali Khan et al. (2000), and which was used for the sequence of
$N_{\text{f}}=3$ simulations Jin et al. (2015, 2017); Kuramashi et al. (2020)
displayed in Figure 4 (left). These data have recently been reanalysed in
order to test for tricritical scaling Cuteri et al. (2021), in analogy to the
staggered case above. To this end, one can exploit that in chiral perturbation
theory

$m_{PS}^{2}\propto m_{q}\;,$ (15)

which also holds for Wilson fermions provided that an additively renormalised
quark mass is used. The data from Figure 4 are replotted in Figure 7 (right),
with the vertical axis rescaled to be proportional to the renormalised quark
mass as scaling variable. Perfect triritical scaling is observed, with
extrapolations that are robust under variation between leading and next-to-
leading order, as well as when using only the coarsest lattices
$N_{\tau}=4,6,8$.

Figure 7: The chiral critical surface projected onto the
$(am,aT=N_{\tau}^{-1})$ plane for unimproved staggered fermions Cuteri et al.
(2021) (left) and for $O(a)$-improved Wilson fermions Kuramashi et al. (2020)
(right), with tricritical scaling fits to both. From Cuteri et al. (2021).

###  3.8 Conclusions for the Continuum Limit

Knowing the structure of the bare parameter lattice phase diagram, its
implications for continuum physics can be stated. The continuum limit is
represented by the origin of the $(am,N_{\tau}^{-1})$ plane, which can be
approached along different lines of constant physics, representing theories
with different values for, e.g., the vacuum pseudo-scalar mass. For any
$N_{\text{f}}$ with a finite $N_{\tau}^{\mathrm{tric}}(N_{\text{f}})$, the
parameter region of first-order transitions is not continuously connected to
the continuum limit, and thus represents a mere cutoff effect. In other words,
for all $N_{\tau}>N_{\tau}^{\mathrm{tric}}$, simulations will only show
crossover behaviour for any vacuum pion mass, such that the transition in the
continuum chiral limit is approached from the crossover, and corresponds to an
isolated second-order point. By contrast, a first-order transition in the
chiral continuum limit implies a finite $m_{c}$ in physical units, which would
require the $Z(2)$-critical line $am_{c}(N_{\tau})$ to terminate in the origin
of the plot. Moreover, there would be no tricritical point and no associated
tricritical scaling, so that $am_{c}(N_{\tau})$ should have an ordinary Taylor
expansion about zero with the usual cutoff corrections $\sim a,a^{2}$ etc. In
Cuteri et al. (2021) it was found that polynomial next-to-leading order fits
are unable to accommodate the $N_{\tau}=4,6,8$ data for $N_{\text{f}}=5,6,7$
unimproved staggered fermions. For the $N_{\text{f}}=3$ $O(a)$-improved Wilson
fermions such fits are possible but have worse $\chi^{2}/\mathrm{dof}$ than
the tricritical fits. We must then conclude that the staggered action predicts
all massless theories with $N_{\text{f}}\leq 7$ to have second-order
transitions in the continuum limit, and that the Wilson clover-action agrees
with this for $N_{f}=3$. Note that this conclusion is perfectly consistent
with the absence of any first-order region in all simulations with improved
staggered actions so far.

In summary:

  * •

Qualitatively, the discretisation effects on the chiral phase transition are
the same for unimproved staggered fermions and either unimproved or
$O(a)$-improved Wilson fermions, making the transition stronger compared to
the continuum.

  * •

In both discretisations, the boundaries of the first-order transitions at
finite lattice spacings exhibit tricritical scaling and extrapolate to a
finite $N_{\tau}^{\mathrm{tric}}(N_{\text{f}}=3)$. This implies that the
first-order region is not connected to the continuum limit. Thus, when the
continuum limit is taken before the chiral limit, as is necessary to avoid
lattice artefacts, both predict a second-order transition.

  * •

Quantitatively, the cutoff effects are larger for Wilson fermions, resulting
in a larger $N_{\tau}^{\mathrm{tric}}(N_{\text{f}})$ than in the case of
staggered fermions. This might be explained by the fact that the respective
discretisations break chiral symmetry fully or only partially.

If these conclusions withstand the test of further investigations, they
require a modification of the Columbia plot, as shown in Figure 8, with a
second-order transition all along the $m_{u,d}=0$ limit, independent of the
strange quark mass. Note that the universality class of these second-order
transitions remains open. Since the chiral symmetry is different for the
$N_{\text{f}}=2,3$ boundary cases, one might expect critical exponents to
smoothly cross from one set of values to the other along the second-order
line.

Figure 8: The Columbia plot in the continuum, as suggested by the tricritical
scaling analyses of unimproved staggered and $O(a)$-improved Wilson fermions.
From Cuteri et al. (2021).

Numerical simulations are no mathematical proof, and the reported results are
in contradiction with the predictions from Pisarski and Wilczek (1984); Butti
et al. (2003). Given the presently available data, what would it take to avoid
this conclusion and render the chiral phase transition for $N_{\text{f}}\geq
3$ first order? This would require future data points for $am_{c}(N_{\tau})$
(or $am_{PS}(N_{\tau})$) on larger $N_{\tau}$ to extrapolate to zero as a
polynomial in $N_{\tau}^{-1}$, i.e. without tricritical scaling, independent
of the discretisation scheme. The tricritical scaling observed for the
presently available $N_{\tau}$ would then have to be “accidental”. On the
other hand, the question arises whether the three-dimensional sigma models
investigated in Pisarski and Wilczek (1984); Butti et al. (2003) represent the
most general case compatible with QCD. It is well-known that
$\phi^{4}$-theories do not permit tricritical points, which requires
$\phi^{6}$-terms. In three dimensions, $\phi^{6}$ theories are renorrmalisable
and moreover show stable infrared fixed points in non-perturbative studies
Tetradis and Litim (1996); Del Debbio and Keegan (2011), which could reconcile
effective theories with the current lattice results. It should be possible to
resolve these questions in the near future and finally settle the fate of the
thermal QCD phase transition with massless quarks.

##  4 The Columbia Plot with Chemical Potential

Once a chemical potential $\mu$ for quark (or $\mu_{B}=3\mu$ for baryon)
number is switched on, the Columbia plot in the continuum (or at fixed lattice
spacing viz. $N_{\tau}$) gets extended into a third dimension. In this case
the $Z(2)$-critical boundary lines sweep out critical surfaces, as sketched
schematically in Figure 9 (left). Labelling the third axis by $(\mu/T)^{2}$,
both real and imaginary chemical potential can be discussed. Imaginary
chemical potential, $\mu=i\mu_{i}$ with $\mu_{i}\in\mathbb{R}$, is unphysical,
but it does not induce a sign problem and thus can be simulated without
difficulty Alford et al. (1999). This can be exploited to extract various
properties of the phase diagram at sufficiently small real $\mu$ by analytic
continuation de Forcrand and Philipsen (2002); D’Elia and Lombardo (2003). Two
exact symmetries are important here. Because of $CP$-invariance, the QCD
partition function is an even function of the chemical potential
$Z(\mu)=Z(-\mu)$. Furthermore, for arbitrary fermion masses it is periodic in
imaginary chemical potential because of the global Roberge-Weiss (or center)
symmetry Roberge and Weiss (1986),

$Z\left(T,i\frac{\mu_{i}}{T}\right)=Z\left(T,i\frac{\mu_{i}}{T}+i\frac{2\pi
n}{N_{c}}\right)\;.$ (16)

Thus, as imaginary chemical potential is increased, the partition function
periodically cycles through the $N_{c}$ center sectors, which are
distinguishable by the phase of the Polyakov loop, whereas the thermodynamic
functions are invariant under a shift between sectors. The boundaries of these
sectors, located at $(\mu/T)_{c}=\pm i(2n+1)\pi/3$ with $n=0,1,2,\ldots$, are
marked by first-order transitions for high temperatures and crossover for low
temperatures, see Figure 9 (right). Varying $T$ at fixed imaginary chemical
potential, the analytic continuation of the QCD thermal transition is crossed,
whose order depends on $N_{f}$ and the quark masses, as discussed in the
previous sections. For a first-order chiral or deconfinement transition, the
transition lines meet up in a triple point, which marks a three-phase
coexistence between two phases of the Polyakov loop and their average. For a
thermal crossover the Roberge-Weiss transition ends in a critical endpoint
with 3D Ising universality. The boundary between these cases, corresponding to
specific quark mass values, is marked by a tricritical point. In the Columbia
plot for the first Roberge-Weiss plane at $\mu/T=i\pi/3$, these trace out
tricritical lines, in which the deconfinement critical surface ends in the
Roberge-Weiss plane de Forcrand and Philipsen (2010). These structures have
been established explicitly for unimproved staggered de Forcrand and Philipsen
(2010); Bonati et al. (2011) as well as unimproved Wilson Philipsen and Pinke
(2014) fermions and generalise to any discretisation respecting the center
symmetry.

Figure 9: (Left): Columbia plot with chemical potential as observed on coarse
lattices. The bottom plane corresponds to the first center-transition.
(Right): The QCD phase diagram at imaginary chemical potential. Vertical lines
mark first-order transitions between different center sectors, the dotted
lines are the analytic continuation of the thermal transitions at real $\mu$,
whose nature depends on the quark masses.

###  4.1 The Deconfinement Transition

Once again it is instructive to also consider the more easily accessible heavy
mass corner, and to study how the deconfinement transition changes in the
presence of a baryon chemical potential. Moreover, for heavy quarks one can
expand the fermion determinant in inverse quark mass and, after an additional
resummed strong coupling expansion, an effective lattice theory can be
constructed that reproduces the $Z(2)$-boundary on $N_{\tau}=4$ lattices at
$\mu=0$., but in addition can be solved at finite chemical potential, both
imaginary and real Fromm et al. (2012). The tricritical line in the Roberge-
Weiss plane dictates the deconfinement critical surface to emanate from it by
tricritical scaling,

$\frac{m_{c}(\mu)}{T}=\frac{m_{c}(0)}{T}+K\Big{[}\Big{(}\frac{\pi}{3}\Big{)}^{2}+\Big{(}\frac{\mu}{T}\Big{)}^{2}\Big{]}^{2/5}\;,$
(17)

where now the deviation of the chemical potential from the Roberge-Weiss plane
is the center symmetry breaking scaling field. This scaling behaviour is
confirmed numerically and, even at leading order, found to extend far into the
range of real chemical potentials. The region of first-order deconfinement
transitions is thus found to continuously shrink with real chemical potential
de Forcrand and Philipsen (2010); Fromm et al. (2012). This observation
implies that the chiral and deconfinement transition regions remain separated
for all chemical potentials, i.e., the respective critical surfaces do not
join. The chiral critical surface must then form a closed volume with the
boundaries of the quark mass and chemical potential parameter space, since one
cannot get from a first-order transition to a crossover without passing
through a critical point.

###  4.2 The Chiral Transition

For unimproved Wilson and staggered discretisations on $N_{\tau}=4,6$, the 3D
Columbia plot looks as in Figure 9 (left), with the region of chiral phase
transitions getting wider in the imaginary $\mu$ direction. Expanding the
$Z(2)$-critical light mass for small chemical potential for any fixed $m_{s}$,

$\frac{m^{c}_{l}(\mu)}{m^{c}_{l}(0)}=1+c_{1}\Big{(}\frac{\mu}{T}\Big{)}^{2}+O\Big{(}\Big{(}\frac{\mu}{T}\Big{)}^{4}\Big{)}\;,$
(18)

one may conclude that $c_{1}<1$ and by analytic continuation the first order
region shrinks in the real-$\mu$ direction. As an independent check on these
calculations at imaginary chemical potential, the curvature of the chiral
critical surface can be computed directly at $\mu=0$. For staggered fermions
one finds $c_{1}<0$ both at $m_{s}=m_{l}$ and at
$m_{s}=m_{s}^{\mathrm{phys}}$, and the next coefficient is negative as well de
Forcrand and Philipsen (2007); de Forcrand et al. (2007); de Forcrand and
Philipsen (2008). Thus the chiral transition strengthens with imaginary and
weakens with real chemical potential. This is opposite to a scenario with a
chiral critical point close to the temperature axis, which would require the
chiral transition at the physical point to strengthen with real $\mu$. The
only investigation finding a positive curvature so far is for $O(a)$-improved
Wilson fermions Jin et al. (2015). However, we have seen in the previous
sections the enormous cutoff effects on the location of the chiral critical
surface at $\mu=0$, and the question is how these vary with $\mu$ for a given
discretisation. The fact that for $\mu=0$ the critical surface moves to the
chiral limit makes it practically impossible to compute its curvature on fine
lattices.

Investigations in the Roberge-Weiss plane on finer lattices reveal the same
trend as seen at $\mu=0$, namely the chiral tricritical line moving towards
smaller quark masses, both for unimproved staggered Philipsen and Sciarra
(2020) and Wilson Czaban et al. (2016) quarks, Figure 10 (left). For stout-
smeared staggered Bonati et al. (2019) and HISQ Goswami et al. (2019) actions
on $N_{\tau}=4$, even the larger first-order region in the Roberge-Weiss plane
cannot be detected when starting from the physical point and reducing the pion
masses down to $m_{\pi}\approx 50$ MeV, as Figure 10 (right) shows by
demonstrating second-order scaling for the Roberge-Weiss endpoint.

$m_{PS}^{\mathrm{tric}}$/MeV

Figure 10: (Left): Tricritical pseudo-scalar mass values delimiting the first-
order chiral region in the RW-plane Philipsen and Sciarra (2020). (Right):
Finite size scaling with 3d $Z(2)$ exponents for stout-smeared staggered
fermions with $m_{PS}\approx 50$ MeV and $(m_{l}/m_{s})_{\mathrm{phys}}$ on
$N_{\tau}=4$ in the RW-plane Bonati et al. (2019).

Note that the order parameter for the center sectors is the imaginary part of
the Polyakov loop, even in the light quark regime. It is then interesting to
study the interplay between the center and chiral dynamics in the Roberge
Weiss plane, since they are independent symmetries. In particular, in the
chiral limit there must be a true phase transition for any chemical potential,
since the latter does not break chiral symmetry. Figure 11 shows the
renormalisation group invariant chiral condensate from Equation (7) and the
disconnected part of the corresponding susceptibility as functions of
temperature, as observed with HISQ fermions on $N_{\tau}=4$ Goswami et al.
(2019). Indeed, as the quark mass decreases a chiral transition is building
up, just as at $\mu=0$, and it is particularly interesting where the
transition happens. The location of the emerging peak in the susceptibility
is, at the current accuracy, consistent with the yellow bar in the figure
marking the critical temperature of the Roberge-Weiss endpoint, so that the
transition lines in the chiral limit would indeed be connected as in Figure 9
(left). Note that there is no a priori reason for this, since we have
independent center and chiral symmetries, with a $Z(2)$ universality of the
Roberge-Weiss endpoint and a larger, $O(N)$-like universality along the chiral
transition line. If these lines really do meet in the chiral limit, their
junction should be multicritical.

Figure 11: (Left): Chiral condensate as defined in Equation (7) in the RW
plane from HISQ fermions on $N_{\tau}=4$. (Right): Disconnected chiral
susceptibility at fixed quark mass ratio for various lattice sizes. Yellow
bands show the location of the RW endpoint transition temperature, extracted
from the peak location of $\chi_{L}$. From Goswami et al. (2019).

In summary, from the $\mu=0$ and imaginary chemical potential results together
we conclude that the entire chiral critical surface is shifting drastically
towards the chiral limit as the lattice spacing is decreased, and it is an
open question whether any first-order transition remains in the continuum
limit at any $\mu$. At the same time, this implies that the crossover at the
physical point for $\mu/T=\mu_{B}/(3T)\raise
1.29167pt\hbox{$<$\kern-7.5pt\raise-4.73611pt\hbox{$\sim$}}1$, which is where
an analytic continuation of Equation (18) holds (see Section 5.1), is very
soft and does not display any singular behaviour due to the chiral transition.

##  5 QCD at the Physical Point

Having discussed the order of the thermal QCD transition in a generalised
parameter space, the aim is to see how this constrains the physical quark mass
configuration. The qualitative relation between the phase diagram in the
chiral limit of the light quarks and their physical values has been worked out
some time ago by mean field methods Halasz et al. (1998); Hatta and Ikeda
(2003), as shown in Figure 12. In the chiral limit, there must be a non-
analytic phase transition for any value of the baryon chemical potential,
since the latter respects chiral symmetry. We have seen that lattice results
are indeed converging towards a second-order transition at $\mu_{B}=0$. By
contrast, the first-order nature of the chiral transition at $T=0$ is based on
expectations from NJL, Gross-Neveu and sigma models in mostly mean field
investigations, with no firm QCD predictions available yet. If this assumption
holds, there must be a tricritical point marking the change between first-
order and second-order behaviour, which might even have an influence on
physical QCD if the light quark masses happen to be in the scaling window of
the tricritical point Stephanov et al. (1998). At finite quark masses chiral
symmetry is broken explicitly, the second-order line turns into crossover and
the first-order transition weakens to disappear in a $Z(2)$-critical line,
which emanates from the tricritical point by tricritical scaling Hatta and
Ikeda (2003),

$\displaystyle\mu_{B}^{c}(m_{u,d})$ $\displaystyle=$
$\displaystyle\mu_{B}^{\mathrm{tric}}+A\;m_{u,d}^{2/5}+O(m_{u,d}^{4/5}),$
$\displaystyle T_{c}(m_{u,d})$ $\displaystyle=$ $\displaystyle
T^{\mathrm{tric}}+B\;m_{u,d}^{2/5}+O(m_{u,d}^{4/5})\;.$ (19)

This implies an ordering of the critical temperatures to be exploited below,

$T_{c}(m_{u,d}=0,\mu_{B}=0)>T_{\mathrm{tric}}(m_{u,d}=0,\mu_{B}=0)>T_{\mathrm{cep}}(m_{u,d}^{\mathrm{phys}},\mu_{B}^{\mathrm{cep}})\;.$
(20)

For completeness, we need to also discuss an alternative scenario, where the
chiral phase transition in the massless limit is second order all the way to
$T=0$. At least from a lattice perspective, this is not excluded so far, but
crucially depends on whether there is any non-trivial
$m^{c}_{u,d}(\mu)$-dependence in the continuum limit. Moreover, a recent
investigation of the chiral nucleon-meson and chiral quark-meson models finds
the phase transition for $m=0$ at $T=0$ to turn second order, once
fluctuations are included Brandes et al. (2021). In such a scenario there is
no tricritical point and no first-order transition anywhere. Instead, non-
vanishing quark masses remove the entire second-order line and the chiral
transition would be analytic crossover exclusively for physical quark masses.

Figure 12: (Left): Relation of the tentative QCD phase diagram with physical
light quark masses (back plane) to the chiral limit (front plane) according to
Halasz et al. (1998); Hatta and Ikeda (2003). (Right): If the entire chiral
transition line in the massless limit is of second order, the transition at
the physical point is crossover everywhere.

###  5.1 The Crossover at Small Baryon Densities

There are several methods that have been used so far to extract information
about the phase structure at the physical point for small baryon density. All
of them introduce some approximation which can be controlled as long as
$\mu/T\raise 1.29167pt\hbox{$<$\kern-7.5pt\raise-4.73611pt\hbox{$\sim$}}1$:
(i) Reweighting Fodor and Katz (2002), (ii) Taylor expansion in $\mu/T$ Allton
et al. (2002) and (iii) analytic continuation from imaginary chemical
potential de Forcrand and Philipsen (2002); D’Elia and Lombardo (2003). When
the QCD pressure is expressed as a series in baryon chemical potential,

$\frac{p(T,\mu_{B})}{T^{4}}=\frac{p(T,0)}{T^{4}}+\sum_{n=1}^{\infty}\frac{1}{2n!}\chi^{B}_{2n}(T)\left(\frac{\mu_{B}}{T}\right)^{2n}\;,\quad\chi^{B}_{2n}(T)=\frac{\partial^{2n}(\frac{p}{T^{4}})}{\partial(\frac{\mu_{B}}{T})^{2n}}\Big{|}_{\mu_{B}=0}\;,$
(21)

the Taylor coefficients are the baryon number fluctuations evaluated at zero
density, which can also be computed by fitting to untruncated results at
imaginary $\mu_{B}$. This permits full control of the systematics between (ii)
and (iii). These coefficients are presently known up to $2n=8$ on
$N_{\tau}=16$ lattices, Figure 13 (left), and in principle also observable
experimentally. For a review of the equation of state relating to heavy ion
phenomenology, see Karsch et al. (2016); Ratti (2019). Note also, that this
low density regime appears to be accessible by complex Langevin simulations
without recourse to series expansions, albeit not yet for physical quark
masses Attanasio et al. (2020). This offers an additional cross check between
different methods.

Figure 13: (Left): Baryon number fluctuations
$\chi^{B}_{2},\chi^{B}_{4},\chi^{B}_{8}$ from the lattice in comparison with
the CEM model. (Right): $b_{1}$ computed directly from Equation (26) by the WB
collaboration, and reverse engineered using CEM from HotQCD baryon number
fluctuations. From Vovchenko et al. (2019).

An important quantity is the pseudo-critical temperature marking the “phase
boundary” between the chirally broken and restored regimes. Since the chiral
transition at the physical point corresponds to an analytic crossover with a
non-zero order parameter everywhere, there are no truly distinct “phases” and
no unambiguous definition of a transition temperature exists. In general,
definitions based on different observables will give different pseudo-critical
temperatures, even in the thermodynamic limit, contrary to the unique
locations of singularities for true phase transitions. While this is an issue
when comparing with an experimental situation, for theoretical investigations
it is convenient to stick to the observables representing the true order
parameter in the appropriate limit, i.e., the susceptibility of an
appropriately normalised chiral condensate in this case. Following as an
implicitly defined function from the partition function, the pseudo-critical
temperature can be similarly expressed as a power series in chemical
potential,

$\frac{T_{pc}(\mu_{B})}{T_{pc}(0)}=1-\kappa_{2}\left(\frac{\mu_{B}}{T_{pc}(0)}\right)^{2}+\ldots,$
(22)

with $T_{pc}(0)=156.5(1.5)$ MeV Bazavov et al. (2019). Continuum extrapolated
results for the leading coefficient are collected in Table 22, the sub-leading
coefficient $\kappa_{4}$ is compatible with zero at the current accuracy. This
is a remarkable result telling us that up to $\mu_{B}\raise
1.29167pt\hbox{$<$\kern-7.5pt\raise-4.73611pt\hbox{$\sim$}}3T$ the dependence
of thermodynamic quantities on chemical potential is rather weak and can be
accurately described by a truncated leading-order Taylor series in chemical
potential. [H] Summary of continuum-extrapolated values for $\kappa_{2}$ in
Equation (22) . $\kappa_{2}$ Action Ref. 0.0158(13) imag. $\mu$, stout-smeared
staggered Bellwied et al. (2015) 0.0135(20) imag. $\mu$, stout-smeared
staggered Bonati et al. (2019) 0.0145(25) Taylor, stout-smeared staggered
Bonati et al. (2019, 2018) 0.016(5) Taylor, HISQ Bazavov et al. (2019)

We now have the necessary information to obtain a conservative bound on the
location of a possible critical point, which according to Figure 12 sits on
the pseudo-critical line of a strengthening crossover. Using the central value
from Equation (10) for the chiral critical temperature and imposing the model-
independent ordering $T_{\mathrm{cep}}<T_{c}=132$ MeV, the chemical potential
of a critical point must satisfy

$\mu_{B}^{\mathrm{cep}}>3.1\;T_{pc}(0)\approx 485\;\mathrm{MeV}.$ (23)

###  5.2 The Search for a Critical Point

For any power series of a function with a given domain of analyticity in its
complex argument, the radius of convergence gives the distance between the
expansion point and the nearest singularity. This implies that the location
$(T_{c},\mu_{B}^{c})$ of a non-analytic QCD phase transition constitutes an
upper bound on the radius of convergence of the series Equation (21), or that
of any other thermodynamic function. Turning this around one may search for a
critical point: If a finite radius of convergence can be extracted from the
pressure series for real parameter values, it should signal a phase
transition. The simplest estimator is the ratio test of consecutive
coefficients, whose extrapolation yields the radius of convergence,

$r=\lim_{n\rightarrow\infty}r_{2n}\;,\quad
r_{2n}=\left|\frac{2n(2n-1)\chi^{B}_{2n}}{\chi^{B}_{2n+2}}\right|\;.$ (24)

In practice, however, only the first few coefficients are available and an
extrapolation is not feasible. For a compilation of available results, see
Bazavov et al. (2017). Moreover, the ratio estimator is inappropriate for
series with irregular signs and complex singularities, where it fails to
converge.

This can be illustrated by modelling the lattice data in a spirit similar to
the hadron resonance gas descriptions, such that higher coefficients become
available and different scenarios can be tested for compatibility with the
data. As an example, consider the fugacity expansion of baryon number density.
At imaginary chemical potential, this is a Fourier series whose coefficients
can be computed on the lattice without sign problem,

$\displaystyle\frac{n_{B}}{T^{3}}|_{\mu_{B}=i\theta_{B}T}$ $\displaystyle=$
$\displaystyle i\sum_{k}b_{k}(T)\sin(k\mu_{B}/T)\;,$ (25) $\displaystyle\
b_{k}(T)$ $\displaystyle=$
$\displaystyle\frac{2}{\pi}\int_{0}^{\pi}d\theta_{B}\;{\rm
Im}\Big{(}\frac{n_{B}(T,i\theta_{B}T)}{T^{3}}\Big{)}\sin(k\theta_{B})\;.$ (26)

In Vovchenko et al. (2018) a cluster expansion model (CEM) was proposed, which
takes the first two coefficients as input from a lattice calculation Vovchenko
et al. (2017), and expresses all higher coefficients in terms of these,

$b_{k}(T)=\alpha_{k}^{SB}\frac{[b_{2}(T)]^{k-2}}{[b_{1}(T)]^{k-1}}\;,\quad
k=3,4,\ldots$ (27)

The $\alpha_{k}^{SB}$ are $T$-independent and fixed to reproduce the Stefan-
Boltzmann limit. This recursion implies that only two-body interactions are
included, and corresponds to a truncated virial expansion which is valid for
sufficiently dilute systems. The model now predicts all coefficients $b_{k\geq
3}$ and allows for an all-order closed expression,

$\displaystyle\frac{n_{B}(T,\mu_{B})}{T^{3}}=-\frac{2}{27\pi^{2}}\,\frac{\hat{b}_{1}^{2}}{\hat{b}_{2}}\,\left\\{4\pi^{2}\,\left[\textrm{Li}_{1}\left(x_{+}\right)-\textrm{Li}_{1}\left(x_{-}\right)\right]+3\,\left[\textrm{Li}_{3}\left(x_{+}\right)-\textrm{Li}_{3}\left(x_{-}\right)\right]\right\\}\;,$
(28)
$\displaystyle\mbox{with}\quad\hat{b}_{1,2}=\displaystyle\frac{b_{1,2}(T)}{b_{1,2}^{\rm
SB}},\quad
x_{\pm}=\displaystyle-\frac{\hat{b}_{2}}{\hat{b}_{1}}\;e^{\pm\mu_{B}/T},\quad\textrm{Li}_{s}(z)=\displaystyle\sum_{k=1}^{\infty}\frac{z^{k}}{k^{s}}.$

All existing lattice data in the crossover regime are reproduced with
excellent accuracy, as the examples in Figure 13 (left) show. Moreover, Figure
13 (right) compares the coefficient $b_{1}$ computed directly from its
defining equation, Equation (26), with a calculation from a combination of
$\chi^{B}_{i}$, by inverting the CEM. Note that all orders of the fugacity
expansion enter this calculation. The remarkable quantitative agreement is
only possible, if both lattice calculations (using different actions) give
equivalent results and all coefficients of the CEM are sufficiently close to
the true QCD values.

Figure 14: (Left): Comparison of estimators for the radius of convergence in
CEM. (Right): The resulting radius of convergence as a function of $T$
predicts the RW-transition at imaginary $\mu_{B}$. From Vovchenko et al.
(2018).

With all coefficients of the fugacity expansion available, one can study the
radius of convergence of CEM, Figure 14. The ratio estimator fails to
converge, because of the irregular signs of higher order coefficients (it
works for equal or alternating signs). On the other hand, the Mercer-Roberts
estimator,

$r_{n}=\left|\frac{c_{n+1}c_{n-1}-c_{n}^{2}}{c_{n+2}c_{n}-c_{n+1}^{2}}\right|^{1/4}\;,$
(29)

converges and extrapolates to a unique radius of convergence, even for
coefficients $c_{n}$ pertaining to different observables, as must be the case
for a true singularity. A number of improved estimators with faster
convergence properties has been proposed in Giordano and Pásztor (2019).

The red points in Figure 14 (right) show the resulting radius of convergence
of the CEM for different temperatures. At high temperatures a singularity is
predicted at a distance $|\mu_{B}/T|\approx\pi$ and $T>T_{pc}$, which is
quantitatively compatible with the first-order Roberge-Weiss transition in the
imaginary $\mu_{B}$ direction. At lower temperatures there is no Roberge-Weiss
transition and the radius of convergence quickly grows. This is an intriguing
result, given that no information from imaginary chemical potential has been
used as input. It shows that the general approach to detect a non-analytic
phase transition is viable, provided that a good estimator for the radius of
convergence is used and that sufficiently many coefficients are available
(more than 100 in this case!). The correct identification of the Roberge-Weiss
transition by the radius of convergence conversely implies that the CEM of
lattice QCD does not have another phase transition closer to the origin, i.e.,
a possible critical endpoint in the real direction must satisfy

$\mu_{B}^{\mathrm{cep}}>\pi T\;.$ (30)

This is fully consistent with Equation (23) but derived by completely
different methods.

Of course, modelling higher coefficients in terms of lower ones is not unique.
The simplest alternative is a description of the baryon number fluctuations up
to $\chi^{B}_{8}$ by a polynomial model, equally without singularity Fodor et
al. (2019). Another one is provided by a rational function model, which can
account for singularities and can be applied both to QCD or to a chiral model
with a phase transition Almási et al. (2019). Equations of state compatible
with lattice data and including an Ising critical point in predefined
locations have also been constructed Parotto et al. (2020). While the
properties of models are not those of QCD, these analyses altogether do show
that there is no sign of criticality in the real $\mu_{B}$ direction from the
presently available, continuum extrapolated lattice data at zero or imaginary
chemical potential, but at best in their higher order completions. It will
thus be interesting to test higher order coefficients without modelling by
Padé approximants or conformal maps as demonstrated in a Gross-Neveu model
Basar (2021), or by resummation schemes in $(\mu_{B}/T)$ Mondal et al. (2021).

Having to go term by term in an expansion can be avoided, if the radius of
convergence is instead determined by the Lee-Yang zero Yang and Lee (1952);
Lee and Yang (1952) closest to the origin. Using reweighting to real chemical
potential, this was the strategy employed in the first lattice prediction of a
critical point on $N_{\tau}=4$ lattices using unimproved rooted staggered
fermions Fodor and Katz (2002). However, it is now understood that the closest
Lee-Yang zero was caused by a spectral gap between the unrooted taste quartets
of the zeros, after Taylor expanding the reweighting factor, rather than by a
phase transition Giordano et al. (2020). This is related to the general
problem of staggered taste quartets splitting up when they cross a branch cut
of a rooted determinant Golterman et al. (2006). It has then been proposed to
redefine the rooted staggered determinant by a geometric matching procedure
averaging over the quartets, after which it can be represented as an ordinary
polynomial Giordano et al. (2020). A further development Giordano et al.
(2020) concerns the reweighting procedure, which usually is based on sampling
with the phase quenched determinant and reweighting in the phase factor. This
has a well-known overlap problem between the sampled and reweighted ensembles.
To get rid of this, one may neglect the imaginary part of the determinant in
the partition function altogether, which is allowed since it will average to
zero. One can then sample with the real part of the fermion determinant and
reweight in the sign only, which has no overlap problem. Of course, the sign
problem remains and one is still faced with a challenging signal to noise
ratio.

Application of these new methods using the stout-smeared action on coarse
$N_{\tau}=4$ lattices appears to signal a Lee-Yang zero at $\mu_{B}\sim 2.4T$
Giordano et al. (2020), which however is is still far from the continuum.
Simulation results on $N_{\tau}=6$ for the renormalised chiral condensate,

$\langle\bar{\psi}\psi\rangle_{R}(T,\mu)=-\frac{m_{ud}}{f_{\pi}^{4}}\left[\langle\bar{\psi}\psi\rangle_{T,\mu}-\langle\bar{\psi}\psi\rangle_{0,0}\right]\,,$
(31)

are shown in Figure 15, both for real and imaginary chemical potential
Borsanyi et al. (2021). In the left plot, there is no steepening or narrowing
of the chiral crossover yet as the real chemical potential is increased. In
the right plot the chemical potential is varied for fixed temperature $T=140$
MeV. One observes full compatibility of the reweighted real $\mu_{B}$
simulations with the analytic continuation from imaginary $\mu_{B}$
simulations, but with smaller errors. In this case, also, there is no sign of
a non-analyticity, which is again consistent with the results from other
methods.

Figure 15: Renormalised chiral condensate from simulations with stout-smeared
staggered fermions at imaginary (blue) and real (black) chemical potential.
(Left): Temperature scan at various chemical potentials. (Right): $\mu_{B}$
scan at at $T=140$ MeV. Coloured bands result from analytic continuation of
imaginary chemical potential data. From Borsanyi et al. (2021).

##  6 Conclusions

Even though the sign problem of lattice QCD remains unsolved, there has been
considerable progress over the last few years towards phenomenologically
relevant constraints on the QCD phase diagram. This is due to a number of
complementary paths of investigation, each of them refining their methods and
having finer lattices as well as different lattice actions at their disposal.

In the Columbia plot at $\mu_{B}=0$, the region of first-order deconfinement
transitions in the heavy mass corner can be directly simulated, and its
$Z(2)$-critical boundary for $N_{\text{f}}=2$, while not yet continuum
extrapolated, is expected to settle in a region around $m_{PS}\sim 4$ GeV. The
region of first-order chiral phase transitions seen in various lattice
discretisations has been observed to shrink continuously with decreasing
lattice spacing. There is numerical evidence for several $N_{\text{f}}\geq 3$
that their $Z(2)$-critical boundary extrapolates to the chiral limit at finite
lattice spacing, exhibiting tricritical scaling. Such first-order regions are
not connected to the continuum and represent lattice artefacts. If this is
confirmed by further studies, the Columbia plot will look like in Figure 8.

In another recent development, the critical temperature for the chiral
transition in the massless limit has been determined to be in the range
$T\approx 126-140$ MeV, which is significantly lower than the pseudo-critical
temperature at the physical point. According to the expected relationship
between the massless limit and the physical quark masses, this bounds a
possible critical endpoint at finite baryon density to be at $\mu_{B}$
$>$$\sim$ $3T$. An estimate of the radius of convergence, based on a cluster
expansion model consistent with all available zero density baryon number
fluctuations for physical QCD, is fully consistent with this. Novel techniques
to extract the leading Lee-Yang zero of the partition function, either with
improved reweighting methods or by refined series analyses, should be able to
provide additional, independent estimates of the proximity of a phase
transition to the physical point in the near future.

Parts of the work reported here were funded by the Deutsche
Forschungsgemeinschaft (German Research Foundation), grant number
315477589-TRR 211 (Strong-interaction matter under extreme conditions).

The author declares no conflict of interest.

References

## References

  * Gasser and Leutwyler (1984) Gasser, J.; Leutwyler, H.  Chiral Perturbation Theory to One Loop.  Ann. Phys. 1984, 158, 142, doi:10.1016/0003-4916(84)90242-2.
  * Leutwyler (1994) Leutwyler, H.  Principles of chiral perturbation theory.  arXiv 1994, arXiv:hep-ph/9406283.
  * Gattringer and Lang (2010) Gattringer, C.; Lang, C.B.  Quantum Chromodynamics on the Lattice; Springer: Berlin/Heidelberg, Germany, 2010; Volume 788, doi:10.1007/978-3-642-01850-3.
  * Philipsen (2010) Philipsen, O.  Lattice QCD at non-zero temperature and baryon density.  In Les Houches Summer School: Session 93: Modern Perspectives in Lattice QCD: Quantum Field Theory and High Performance Computing; Oxford Univesity Press: Oxford UK, 2010; pp. 273–330.
  * Aarts (2016) Aarts, G.  Introductory lectures on lattice QCD at nonzero baryon number.  J. Phys. Conf. Ser. 2016, 706, 022004, doi:10.1088/1742-6596/706/2/022004.
  * Montvay and Munster (1997) Montvay, I.; Munster, G.  Quantum Fields on a Lattice; Cambridge Monographs on Mathematical Physics; Cambridge University Press: Cambridge, UK , 1997, doi:10.1017/CBO9780511470783.
  * DeGrand and Detar (2006) DeGrand, T.; Detar, C.E.  Lattice Methods for Quantum Chromodynamics; World Scientific, Singapore, 2006.
  * Symanzik (1983) Symanzik, K.  Continuum Limit and Improved Action in Lattice Theories. 1. Principles and phi**4 Theory.  Nucl. Phys. B 1983, 226, 187–204, doi:10.1016/0550-3213(83)90468-6.
  * Nielsen and Ninomiya (1981) Nielsen, H.B.; Ninomiya, M.  No Go Theorem for Regularizing Chiral Fermions.  Phys. Lett. B 1981, 105, 219–223, doi:10.1016/0370-2693(81)91026-1.
  * Yang and Lee (1952) Yang, C.N.; Lee, T.D.  Statistical theory of equations of state and phase transitions. 1. Theory of condensation.  Phys. Rev. 1952, 87, 404–409, doi:10.1103/PhysRev.87.404.
  * Lee and Yang (1952) Lee, T.D.; Yang, C.N.  Statistical theory of equations of state and phase transitions. 2. Lattice gas and Ising model.  Phys. Rev. 1952, 87, 410–419, doi:10.1103/PhysRev.87.410.
  * de Forcrand (2009) de Forcrand, P.  Simulating QCD at finite density.  PoS 2009, LAT2009, 010, doi:10.22323/1.091.0010.
  * Aoki et al. (2006) Aoki, Y.; Endrodi, G.; Fodor, Z.; Katz, S.D.; Szabo, K.K.  The Order of the quantum chromodynamics transition predicted by the standard model of particle physics.  Nature 2006, 443, 675–678, doi:10.1038/nature05120.
  * Brown et al. (1990) Brown, F.R.; Butler, F.P.; Chen, H.; Christ, N.H.; Dong, Z.h.; Schaffer, W.; Unger, L.I.; Vaccarino, A.  On the existence of a phase transition for QCD with three light quarks.  Phys. Rev. Lett. 1990, 65, 2491–2494, doi:10.1103/PhysRevLett.65.2491.
  * McLerran and Svetitsky (1981) McLerran, L.D.; Svetitsky, B.  Quark Liberation at High Temperature: A Monte Carlo Study of SU(2) Gauge Theory.  Phys. Rev. D 1981, 24, 450, doi:10.1103/PhysRevD.24.450.
  * Boyd et al. (1996) Boyd, G.; Engels, J.; Karsch, F.; Laermann, E.; Legeland, C.; Lutgemeier, M.; Petersson, B.  Thermodynamics of SU(3) lattice gauge theory.  Nucl. Phys. B 1996, 469, 419–444, doi:10.1016/0550-3213(96)00170-8.
  * Shirogane et al. (2016) Shirogane, M.; Ejiri, S.; Iwami, R.; Kanaya, K.; Kitazawa, M.  Latent heat at the first order phase transition point of SU(3) gauge theory.  Phys. Rev. D 2016, 94, 014506, doi:10.1103/PhysRevD.94.014506.
  * Saito et al. (2011) Saito, H.; Ejiri, S.; Aoki, S.; Hatsuda, T.; Kanaya, K.; Maezawa, Y.; Ohno, H.; Umeda, T.  Phase structure of finite temperature QCD in the heavy quark region.  Phys. Rev. D 2011, 84, 054502, doi:black10.1103/PhysRevD.85.079902.
  * Ejiri et al. (2020) Ejiri, S.; Itagaki, S.; Iwami, R.; Kanaya, K.; Kitazawa, M.; Kiyohara, A.; Shirogane, M.; Umeda, T.  End point of the first-order phase transition of QCD in the heavy quark region by reweighting from quenched QCD.  Phys. Rev. D 2020, 101, 054505, doi:10.1103/PhysRevD.101.054505.
  * Cuteri et al. (2021) Cuteri, F.; Philipsen, O.; Schön, A.; Sciarra, A.  Deconfinement critical point of lattice QCD with $N_{f}$=2 Wilson fermions.  Phys. Rev. D 2021, 103, 014513, doi:10.1103/PhysRevD.103.014513.
  * Kiyohara et al. (2021) Kiyohara, A.; Kitazawa, M.; Ejiri, S.; Kanaya, K.  Finite-size scaling around the critical point in the heavy quark region of QCD. _arXiv_ 2021, arXiv:2108.00118.
  * Fischer et al. (2015) Fischer, C.S.; Luecker, J.; Pawlowski, J.M.  Phase structure of QCD for heavy quarks.  Phys. Rev. D 2015, 91, 014024, doi:10.1103/PhysRevD.91.014024.
  * Fromm et al. (2012) Fromm, M.; Langelage, J.; Lottini, S.; Philipsen, O.  The QCD deconfinement transition for heavy quarks and all baryon chemical potentials.  J. High Energy Phys. 2012, 2012, 42, doi:10.1007/JHEP01(2012)042.
  * Fromm et al. (2013) Fromm, M.; Langelage, J.; Lottini, S.; Neuman, M.; Philipsen, O.  Onset Transition to Cold Nuclear Matter from Lattice QCD with Heavy Quarks.  Phys. Rev. Lett. 2013, 110, 122001, doi:10.1103/PhysRevLett.110.122001.
  * Pisarski and Wilczek (1984) Pisarski, R.D.; Wilczek, F.  Remarks on the Chiral Phase Transition in Chromodynamics.  Phys. Rev. 1984, D29, 338–341, doi:10.1103/PhysRevD.29.338.
  * Gausterer and Sanielevici (1988) Gausterer, H.; Sanielevici, S.  Can the Chiral Transition in QCD Be Described by a Linear $\sigma$ Model in Three-dimensions?  Phys. Lett. B 1988, 209, 533–537, doi:10.1016/0370-2693(88)91188-4.
  * Butti et al. (2003) Butti, A.; Pelissetto, A.; Vicari, E.  On the nature of the finite temperature transition in QCD.  J. High Energy Phys. 2003, 2003, 029, doi:10.1088/1126-6708/2003/08/029.
  * Pelissetto and Vicari (2013) Pelissetto, A.; Vicari, E.  Relevance of the axial anomaly at the finite-temperature chiral transition in QCD.  Phys. Rev. D 2013, 88, 105018, doi:10.1103/PhysRevD.88.105018.
  * Braun et al. (2020a) Braun, J.; Fu, W.j.; Pawlowski, J.M.; Rennecke, F.; Rosenblüh, D.; Yin, S.  Chiral susceptibility in ( 2+1 )-flavor QCD.  Phys. Rev. D 2020, 102, 056010, doi:10.1103/PhysRevD.102.056010.
  * Braun et al. (2020b) Braun, J.; Leonhardt, M.; Pawlowski, J.M.; Rosenblüh, D.  Chiral and effective $U(1)_{\rm A}$ symmetry restoration in QCD. _arXiv_ 2020, arXiv:2012.06231.
  * Ding et al. (2019) Ding, H.T.; others.  Chiral Phase Transition Temperature in ( 2+1 )-Flavor QCD.  Phys. Rev. Lett. 2019, 123, 062002, doi:10.1103/PhysRevLett.123.062002.
  * Kaczmarek et al. (2020) Kaczmarek, O.; Karsch, F.; Lahiri, A.; Schmidt, C.  Universal scaling properties of QCD close to the chiral limit. _arXiv_ 2020, arXiv:2010.15593.
  * Iwasaki et al. (1996) Iwasaki, Y.; Kanaya, K.; Kaya, S.; Sakai, S.; Yoshie, T.  Finite temperature transitions in lattice QCD with Wilson quarks: Chiral transitions and the influence of the strange quark.  Phys. Rev. D 1996, 54, 7010–7031, doi:10.1103/PhysRevD.54.7010.
  * Bonati et al. (2014) Bonati, C.; de Forcrand, P.; D’Elia, M.; Philipsen, O.; Sanfilippo, F.  Chiral phase transition in two-flavor QCD from an imaginary chemical potential.  Phys. Rev. 2014, D90, 074030, doi:10.1103/PhysRevD.90.074030.
  * Cuteri et al. (2018) Cuteri, F.; Philipsen, O.; Sciarra, A.  QCD chiral phase transition from noninteger numbers of flavors.  Phys. Rev. D 2018, 97, 114511, doi:10.1103/PhysRevD.97.114511.
  * Philipsen and Pinke (2016) Philipsen, O.; Pinke, C.  The $N_{f}=2$ QCD chiral phase transition with Wilson fermions at zero and imaginary chemical potential.  Phys. Rev. 2016, D93, 114507, doi:10.1103/PhysRevD.93.114507.
  * Kotov et al. (2021) Kotov, A.Y.; Lombardo, M.P.; Trunin, A.  QCD transition at the physical point, and its scaling window from twisted mass Wilson fermions. _arXiv_ 2021, arXiv:2105.09842.
  * Nakamura et al. (2019) Nakamura, Y.; Kuramashi, Y.; Ohno, H.; Takeda, S.  Critical endpoint in the continuum limit and critical endline at $N_{\rm T}=6$ of the finite temperature phase transition of QCD with clover fermions.  PoS 2019, LATTICE2019, 053, doi:10.22323/1.363.0053.
  * Clarke et al. (2021) Clarke, D.A.; Kaczmarek, O.; Karsch, F.; Lahiri, A.; Sarkar, M.  Sensitivity of the Polyakov loop and related observables to chiral symmetry restoration.  Phys. Rev. D 2021, 103, L011501, doi:10.1103/PhysRevD.103.L011501.
  * de Forcrand et al. (2007) de Forcrand, P.; Kim, S.; Philipsen, O.  A QCD chiral critical point at small chemical potential: Is it there or not?  PoS 2007, LAT2007, 178.
  * Varnhorst (2015) Varnhorst, L.  The $N_{f}$=3 critical endpoint with smeared staggered quarks.  PoS 2015, LATTICE2014, 193, doi:10.22323/1.214.0193.
  * Karsch et al. (2001) Karsch, F.; Laermann, E.; Schmidt, C.  The Chiral critical point in three-flavor QCD.  Phys. Lett. 2001, B520, 41–49, doi:10.1016/S0370-2693(01)01114-5.
  * Karsch et al. (2004) Karsch, F.; Allton, C.R.; Ejiri, S.; Hands, S.J.; Kaczmarek, O.; Laermann, E.; Schmidt, C.  Where is the chiral critical point in three flavor QCD?  Nucl. Phys. B Proc. Suppl. 2004, 129, 614–616, doi:10.1016/S0920-5632(03)02659-8.
  * Bazavov et al. (2017) Bazavov, A.; Ding, H.T.; Hegde, P.; Karsch, F.; Laermann, E.; Mukherjee, S.; Petreczky, P.; Schmidt, C.  Chiral phase structure of three flavor QCD at vanishing baryon number density.  Phys. Rev. D 2017, 95, 074505, doi:10.1103/PhysRevD.95.074505.
  * Jin et al. (2015) Jin, X.Y.; Kuramashi, Y.; Nakamura, Y.; Takeda, S.; Ukawa, A.  Critical endpoint of the finite temperature phase transition for three flavor QCD.  Phys. Rev. D 2015, 91, 014508, doi:10.1103/PhysRevD.91.014508.
  * Jin et al. (2017) Jin, X.Y.; Kuramashi, Y.; Nakamura, Y.; Takeda, S.; Ukawa, A.  Critical point phase transition for finite temperature 3-flavor QCD with non-perturbatively O($a$) improved Wilson fermions at $N_{\rm t}=10$.  Phys. Rev. 2017, D96, 034523, doi:10.1103/PhysRevD.96.034523.
  * Kuramashi et al. (2020) Kuramashi, Y.; Nakamura, Y.; Ohno, H.; Takeda, S.  Nature of the phase transition for finite temperature $N_{\rm f}=3$ QCD with nonperturbatively O($a$) improved Wilson fermions at $N_{\rm t}=12$.  Phys. Rev. D 2020, 101, 054509, doi:10.1103/PhysRevD.101.054509.
  * de Forcrand and D’Elia (2017) de Forcrand, P.; D’Elia, M.  Continuum limit and universality of the Columbia plot.  PoS 2017, LATTICE2016, 081.
  * Ohno et al. (2018) Ohno, H.; Kuramashi, Y.; Nakamura, Y.; Takeda, S.  Continuum extrapolation of the critical endpoint in 4-flavor QCD with Wilson-Clover fermions.  PoS 2018, LATTICE2018, 174, doi:10.22323/1.334.0174.
  * Rajagopal (1995) Rajagopal, K.  The Chiral phase transition in QCD: Critical phenomena and long wavelength pion oscillations. _World Sci._ 1995, doi:10.1142/9789812830661_0009.
  * Lawrie and Sarbach (1984) Lawrie, I.; Sarbach, S.  Theory of tricritical points. In Phase Transitions and Critical Phenomena; Domb, C.; Lebowitz, J., Eds.; Academic Press, Cambridge MA, USA, 1984; Volume 9, p. 1.
  * de Forcrand and Philipsen (2007) de Forcrand, P.; Philipsen, O.  The Chiral critical line of N(f) = 2+1 QCD at zero and non-zero baryon density.  J. High Energy Phys. 2007, 0701, 077, doi:10.1088/1126-6708/2007/01/077.
  * Cuteri et al. (2021) Cuteri, F.; Philipsen, O.; Sciarra, A.  On the order of the QCD chiral phase transition for different numbers of quark flavours. _arXiv_ 2021, arXiv:2107.12739.
  * Aoki (1984) Aoki, S.  New Phase Structure for Lattice QCD with Wilson Fermions.  Phys. Rev. D 1984, 30, 2653, doi:10.1103/PhysRevD.30.2653.
  * Aoki (1986) Aoki, S.  A Solution to the U(1) Problem on a Lattice.  Phys. Rev. Lett. 1986, 57, 3136, doi:10.1103/PhysRevLett.57.3136.
  * Sharpe and Singleton (1998) Sharpe, S.R.; Singleton, Jr, R.L.  Spontaneous flavor and parity breaking with Wilson fermions.  Phys. Rev. D 1998, 58, 074501, doi:10.1103/PhysRevD.58.074501.
  * Farchioni et al. (2005) Farchioni, F.; Frezzotti, R.; Jansen, K.; Montvay, I.; Rossi, G.C.; Scholz, E.; Shindler, A.; Ukita, N.; Urbach, C.; Wetzorke, I.  Twisted mass quarks and the phase structure of lattice QCD.  Eur. Phys. J. C 2005, 39, 421–433, doi:10.1140/epjc/s2004-02078-9.
  * Burger et al. (2013) Burger, F.; Ilgenfritz, E.M.; Kirchner, M.; Lombardo, M.P.; Müller-Preussker, M.; Philipsen, O.; Urbach, C.; Zeidlewicz, L.  Thermal QCD transition with two flavors of twisted mass fermions.  Phys. Rev. 2013, D87, 074508, doi:10.1103/PhysRevD.87.074508.
  * Ali Khan et al. (2000) Ali Khan, A.; others.  Phase structure and critical temperature of two flavor QCD with renormalization group improved gauge action and clover improved Wilson quark action.  Phys. Rev. D 2000, 63, 034502, doi:10.1103/PhysRevD.63.034502.
  * Tetradis and Litim (1996) Tetradis, N.; Litim, D.F.  Analytical solutions of exact renormalization group equations.  Nucl. Phys. B 1996, 464, 492–511, doi:10.1016/0550-3213(95)00642-7.
  * Del Debbio and Keegan (2011) Del Debbio, L.; Keegan, L.  RG flows in 3D scalar field theory.  PoS 2011, LATTICE2011, 061, doi:10.22323/1.139.0061.
  * Alford et al. (1999) Alford, M.G.; Kapustin, A.; Wilczek, F.  Imaginary chemical potential and finite fermion density on the lattice.  Phys. Rev. D 1999, 59, 054502, doi:10.1103/PhysRevD.59.054502.
  * de Forcrand and Philipsen (2002) de Forcrand, P.; Philipsen, O.  The QCD phase diagram for small densities from imaginary chemical potential.  Nucl. Phys. B 2002, 642, 290–306, doi:10.1016/S0550-3213(02)00626-0.
  * D’Elia and Lombardo (2003) D’Elia, M.; Lombardo, M.P.  Finite density QCD via imaginary chemical potential.  Phys. Rev. D 2003, 67, 014505, doi:10.1103/PhysRevD.67.014505.
  * Roberge and Weiss (1986) Roberge, A.; Weiss, N.  Gauge Theories With Imaginary Chemical Potential and the Phases of QCD.  Nucl. Phys. B 1986, 275, 734–745, doi:10.1016/0550-3213(86)90582-1.
  * de Forcrand and Philipsen (2010) de Forcrand, P.; Philipsen, O.  Constraining the QCD phase diagram by tricritical lines at imaginary chemical potential.  Phys. Rev. Lett. 2010, 105, 152001, doi:10.1103/PhysRevLett.105.152001.
  * Bonati et al. (2011) Bonati, C.; Cossu, G.; D’Elia, M.; Sanfilippo, F.  The Roberge-Weiss endpoint in $N_{f}$ = 2 QCD.  Phys. Rev. D 2011, 83, 054505, doi:10.1103/PhysRevD.83.054505.
  * Philipsen and Pinke (2014) Philipsen, O.; Pinke, C.  Nature of the Roberge-Weiss transition in $N_{f}=2$ QCD with Wilson fermions.  Phys. Rev. D 2014, 89, 094504, doi:10.1103/PhysRevD.89.094504.
  * de Forcrand and Philipsen (2008) de Forcrand, P.; Philipsen, O.  The Chiral critical point of N(f) = 3 QCD at finite density to the order (mu/T)**4.  J. High Energy Phys. 2008, 11, 012, doi:10.1088/1126-6708/2008/11/012.
  * Jin et al. (2015) Jin, X.Y.; Kuramashi, Y.; Nakamura, Y.; Takeda, S.; Ukawa, A.  Curvature of the critical line on the plane of quark chemical potential and pseudoscalar meson mass for three-flavor QCD.  Phys. Rev. 2015, D92, 114511, doi:10.1103/PhysRevD.92.114511.
  * Philipsen and Sciarra (2020) Philipsen, O.; Sciarra, A.  Finite Size and Cut-Off Effects on the Roberge-Weiss Transition in $N_{\text{f}}=2$ QCD with Staggered Fermions.  Phys. Rev. D 2020, 101, 014502, doi:10.1103/PhysRevD.101.014502.
  * Czaban et al. (2016) Czaban, C.; Cuteri, F.; Philipsen, O.; Pinke, C.; Sciarra, A.  Roberge-Weiss transition in $N_{\text{f}}=2$ QCD with Wilson fermions and $N_{\tau}=6$.  Phys. Rev. 2016, D93, 054507, doi:10.1103/PhysRevD.93.054507.
  * Bonati et al. (2019) Bonati, C.; Calore, E.; D’Elia, M.; Mesiti, M.; Negro, F.; Sanfilippo, F.; Schifano, S.F.; Silvi, G.; Tripiccione, R.  Roberge-Weiss endpoint and chiral symmetry restoration in $N_{f}=2+1$ QCD.  Phys. Rev. D 2019, 99, 014502, doi:10.1103/PhysRevD.99.014502.
  * Goswami et al. (2019) Goswami, J.; Karsch, F.; Lahiri, A.; Neumann, M.; Schmidt, C.  Critical end points in (2+1)-flavor QCD with imaginary chemical potential.  PoS 2019, CORFU2018, 162, doi:10.22323/1.347.0162.
  * Halasz et al. (1998) Halasz, A.M.; Jackson, A.D.; Shrock, R.E.; Stephanov, M.A.; Verbaarschot, J.J.M.  On the phase diagram of QCD.  Phys. Rev. D 1998, 58, 096007, doi:10.1103/PhysRevD.58.096007.
  * Hatta and Ikeda (2003) Hatta, Y.; Ikeda, T.  Universality, the QCD critical / tricritical point and the quark number susceptibility.  Phys. Rev. D 2003, 67, 014028, doi:10.1103/PhysRevD.67.014028.
  * Stephanov et al. (1998) Stephanov, M.A.; Rajagopal, K.; Shuryak, E.V.  Signatures of the tricritical point in QCD.  Phys. Rev. Lett. 1998, 81, 4816–4819, doi:10.1103/PhysRevLett.81.4816.
  * Brandes et al. (2021) Brandes, L.; Kaiser, N.; Weise, W.  Fluctuations and phases in baryonic matter.  Eur. Phys. J. A 2021, 57, 243, doi:10.1140/epja/s10050-021-00528-2.
  * Fodor and Katz (2002) Fodor, Z.; Katz, S.D.  A New method to study lattice QCD at finite temperature and chemical potential.  Phys. Lett. B 2002, 534, 87–92, doi:10.1016/S0370-2693(02)01583-6.
  * Allton et al. (2002) Allton, C.R.; Ejiri, S.; Hands, S.J.; Kaczmarek, O.; Karsch, F.; Laermann, E.; Schmidt, C.; Scorzato, L.  The QCD thermal phase transition in the presence of a small chemical potential.  Phys. Rev. D 2002, 66, 074507, doi:10.1103/PhysRevD.66.074507.
  * Karsch et al. (2016) Karsch, F.; others.  Conserved Charge Fluctuations from Lattice QCD and the Beam Energy Scan.  Nucl. Phys. A 2016, 956, 352–355, doi:10.1016/j.nuclphysa.2016.01.008.
  * Ratti (2019) Ratti, C.  QCD at non-zero density and phenomenology.  PoS 2019, LATTICE2018, 004, doi:10.22323/1.334.0004.
  * Attanasio et al. (2020) Attanasio, F.; Jäger, B.; Ziegler, F.P.G.  Complex Langevin simulations and the QCD phase diagram: Recent developments.  Eur. Phys. J. A 2020, 56, 251, doi:10.1140/epja/s10050-020-00256-z.
  * Bellwied et al. (2015) Bellwied, R.; Borsanyi, S.; Fodor, Z.; Günther, J.; Katz, S.D.; Ratti, C.; Szabo, K.K.  The QCD phase diagram from analytic continuation.  Phys. Lett. B 2015, 751, 559–564, doi:10.1016/j.physletb.2015.11.011.
  * Bonati et al. (2019) Bonati, C.; D’Elia, M.; Negro, F.; Sanfilippo, F.; Zambello, K.  The curvature of the chiral pseudocritical line from LQCD: Analytic continuation and Taylor expansion compared.  Nucl. Phys. A 2019, 982, 198–200, doi:10.1016/j.nuclphysa.2018.08.026.
  * Bonati et al. (2018) Bonati, C.; D’Elia, M.; Negro, F.; Sanfilippo, F.; Zambello, K.  Curvature of the pseudocritical line in QCD: Taylor expansion matches analytic continuation.  Phys. Rev. D 2018, 98, 054510, doi:10.1103/PhysRevD.98.054510.
  * Bazavov et al. (2019) Bazavov, A.; others.  Chiral crossover in QCD at zero and non-zero chemical potentials.  Phys. Lett. B 2019, 795, 15–21, doi:10.1016/j.physletb.2019.05.013.
  * Vovchenko et al. (2019) Vovchenko, V.; Steinheimer, J.; Philipsen, O.; Pasztor, A.; Fodor, Z.; Katz, S.D.; Stoecker, H.  Lattice-based QCD equation of state at finite baryon density: Cluster Expansion Model.  Nucl. Phys. A 2019, 982, 859–862, doi:10.1016/j.nuclphysa.2018.10.068.
  * Bazavov et al. (2017) Bazavov, A.; others.  The QCD Equation of State to $\mathcal{O}(\mu_{B}^{6})$ from Lattice QCD.  Phys. Rev. D 2017, 95, 054504, doi:10.1103/PhysRevD.95.054504.
  * Vovchenko et al. (2018) Vovchenko, V.; Steinheimer, J.; Philipsen, O.; Stoecker, H.  Cluster Expansion Model for QCD Baryon Number Fluctuations: No Phase Transition at $\mu_{B}/T<\pi$.  Phys. Rev. D 2018, 97, 114030, doi:10.1103/PhysRevD.97.114030.
  * Vovchenko et al. (2017) Vovchenko, V.; Pasztor, A.; Fodor, Z.; Katz, S.D.; Stoecker, H.  Repulsive baryonic interactions and lattice QCD observables at imaginary chemical potential.  Phys. Lett. B 2017, 775, 71–78, doi:10.1016/j.physletb.2017.10.042.
  * Giordano and Pásztor (2019) Giordano, M.; Pásztor, A.  Reliable estimation of the radius of convergence in finite density QCD.  Phys. Rev. D 2019, 99, 114510, doi:10.1103/PhysRevD.99.114510.
  * Fodor et al. (2019) Fodor, Z.; Giordano, M.; Günther, J.N.; Kapás, K.; Katz, S.D.; Pásztor, A.; Portillo, I.; Ratti, C.; Sexty, D.; Szabó, K.K.  Trying to constrain the location of the QCD critical endpoint with lattice simulations.  Nucl. Phys. A 2019, 982, 843–846, doi:10.1016/j.nuclphysa.2018.12.015.
  * Almási et al. (2019) Almási, G.A.; Friman, B.; Morita, K.; Redlich, K.  Fourier coefficients of the net baryon number density and their scaling properties near a phase transition.  Phys. Lett. B 2019, 793, 19–25, doi:10.1016/j.physletb.2019.04.023.
  * Parotto et al. (2020) Parotto, P.; Bluhm, M.; Mroczek, D.; Nahrgang, M.; Noronha-Hostler, J.; Rajagopal, K.; Ratti, C.; Schäfer, T.; Stephanov, M.  QCD equation of state matched to lattice data and exhibiting a critical point singularity.  Phys. Rev. C 2020, 101, 034901, doi:10.1103/PhysRevC.101.034901.
  * Basar (2021) Basar, G.  Universality, Lee-Yang singularities and series expansions. _arXiv_ 2021, arXiv:2105.08080.
  * Mondal et al. (2021) Mondal, S.; Mukherjee, S.; Hegde, P.  Lattice QCD Equation of State for Nonvanishing Chemical Potential by Resumming Taylor Expansion. _arXiv_ 2021, arXiv:2106.03165.
  * Fodor and Katz (2002) Fodor, Z.; Katz, S.D.  Lattice determination of the critical point of QCD at finite T and mu.  J. High Energy Phys. 2002, 03, 014, doi:10.1088/1126-6708/2002/03/014.
  * Giordano et al. (2020) Giordano, M.; Kapas, K.; Katz, S.D.; Nogradi, D.; Pasztor, A.  Radius of convergence in lattice QCD at finite $\mu_{B}$ with rooted staggered fermions.  Phys. Rev. D 2020, 101, 074511, doi:10.1103/PhysRevD.101.074511.
  * Golterman et al. (2006) Golterman, M.; Shamir, Y.; Svetitsky, B.  Breakdown of staggered fermions at nonzero chemical potential.  Phys. Rev. D 2006, 74, 071501, doi:10.1103/PhysRevD.74.071501.
  * Giordano et al. (2020) Giordano, M.; Kapas, K.; Katz, S.D.; Nogradi, D.; Pasztor, A.  New approach to lattice QCD at finite density; results for the critical end point on coarse lattices.  J. High Energy Phys. 2020, 05, 088, doi:10.1007/JHEP05(2020)088.
  * Borsanyi et al. (2021) Borsanyi, S.; Fodor, Z.; Giordano, M.; Katz, S.D.; Nogradi, D.; Pasztor, A.; Wong, C.H.  Lattice simulations of the QCD chiral transition at real baryon density. _arXiv_ 2021, arXiv:2108.09213.