# Phase structure of 2+1-flavour QCD and the magnetic equation of state

Fei Gao  Institut für Theoretische Physik, Universität Heidelberg,
Philosophenweg 16, 69120 Heidelberg, Germany  Jan M. Pawlowski  Institut für
Theoretische Physik, Universität Heidelberg, Philosophenweg 16, 69120
Heidelberg, Germany  ExtreMe Matter Institute EMMI, GSI, Planckstr. 1, 64291
Darmstadt, Germany

###### Abstract

We determine the chiral phase structure of $2+1$-flavour QCD in dependence of
temperature and the light flavour quark mass with Dyson-Schwinger equations.
Specifically, we compute the renormalised chiral condensate and its
susceptibility. The latter is used to determine the (pseudo)critical
temperature for general light current quark masses. In the chiral limit we
obtain a critical temperature of about 141 MeV. This result is in quantitative
agreement with recent functional renormalisation group results in QCD, and is
compatible with the respective lattice results. We also compute the order
parameter potential of the light chiral condensate and map out the regime in
the phase diagram which exhibits quasi-massless modes, and discuss the
respective chiral dynamics.

###### pacs:

11.30.Rd, 12.38.Aw, 05.10.Cc, 12.38.Mh, 12.38.Gc

##  I Introduction

The dynamics of QCD matter in heavy ion collisions (HIC) depends crucially on
the chiral properties of QCD at finite temperature and density. For physical
quark masses, QCD exhibits a smooth chiral crossover behaviour from the
perturbative high temperature region of gluons and nearly massless light
quarks to the hadronic phase with a sizeable spontaneous chiral symmetry
breaking, for reviews see Luo and Xu (2017); Adamczyk  _et al._ (2017);
Andronic  _et al._ (2018); Stephanov (2006); Andersen  _et al._ (2016);
Shuryak (2017); Pawlowski (2014); Roberts and Schmidt (2000); Fischer (2019);
Yin (2018).

Naturally, the small current quark mass of the light up and down quarks begs
the question, whether QCD under the extreme conditions of a HIC is close to
the chiral limit. A positive answer allows and supports the phenomenological
access to the dynamics in HICs via systematic chiral expansions, similar to
the extremely successful chiral perturbation theory ($\chi$PT) in the vacuum.
The latter is based on an expansion in the pion mass that works in the regime
with quasi-massless pions. Consequently for the answer to the above
phenomenologically important question we have to map out the regime with
quasi-massless modes in the phase diagram of QCD.

Moreover, in the chiral limit QCD is expected to exhibit a second order phase
transition, and the value $T_{c0}$ of the critical temperature, the respective
universality class as well as the size of the critical regime may lead to
constraints on the location of the potential critical end point (CEP) at
finite density. This has led to detailed studies of the current quark mass
dependence of the chiral crossover, the magnetic equation of state of QCD, for
recent works with functional QCD and lattice QCD, see Ding  _et al._ (2019);
Braun  _et al._ (2020a); Kotov  _et al._ (2020). While not being fully
conclusive, functional studies in low energy effective theories of QCD
suggest, that the critical regime is rather small and is restricted to pion
masses lower than $m_{\pi}\sim 1-10$ MeV, see e.g. Braun and Klein (2008,
2009); Braun  _et al._ (2011a), for a recent review Klein (2017). Given, that
the respective low energy effective theories encode the full chiral dynamics
of QCD, it is unlikely that the inclusion of gluonic fluctuations lead to an
increase of the scaling window, and this is corroborated by a recent study in
functional renormalisation group (fRG) study in first principle QCD Braun  _et
al._ (2020a).

In the present work we access the magnetic equation of state within a
generalised functional approach, that combines Dyson-Schwinger equations (DSE)
and the functional renormalisation group (fRG). This approach was set-up and
used in Gao and Pawlowski (2020, 2021); Gao  _et al._ (2021a) for the phase
structure of QCD as well as quantitative computations of QCD correlation
functions in the vacuum. It extends and utilises previous computations in
first principles QCD with functional approaches, for fRG works see e.g. Braun
(2009); Braun  _et al._ (2011b); Fister and Pawlowski (2011); Mitter  _et al._
(2015); Braun  _et al._ (2016); Rennecke (2015); Fu  _et al._ (2016); Rennecke
and Schaefer (2017); Cyrol  _et al._ (2016, 2018a, 2018b); Fu  _et al._
(2020a, b); Leonhardt  _et al._ (2020); Braun  _et al._ (2020b), for DSE works
see e.g. Roberts and Schmidt (2000); Qin  _et al._ (2011); Fischer  _et al._
(2011); Fischer and Mueller (2011); Fischer  _et al._ (2014a, b); Eichmann
_et al._ (2016a); Gao  _et al._ (2016); Gao and Liu (2016); Fischer (2019);
Gunkel  _et al._ (2019); Isserstedt  _et al._ (2019); Reinosa  _et al._ (2015,
2017); Maelger  _et al._ (2018a, b, 2020); Aguilar  _et al._ (2017, 2018a,
2018b). For related lattice studies see e.g. Bazavov  _et al._ (2012);
Borsanyi  _et al._ (2013, 2014); Bonati  _et al._ (2015); Bellwied  _et al._
(2015); Bazavov  _et al._ (2017a, b); Bonati  _et al._ (2018); Borsanyi  _et
al._ (2018); Bazavov  _et al._ (2019); Guenther  _et al._ (2018); Ding  _et
al._ (2019).

Specifically, we compute the quark and gluon propagators and the quark-gluon
vertex at finite temperature. The quark propagator is then utilised to compute
the renormalised light chiral condensate, whose thermal susceptibility defines
the chiral crossover temperature. With these observables we discuss the
magnetic equation of state in comparison to other functional and lattice
results, and deduce the respective chiral transition temperature in dependence
of the light current quark mass. We also compute the order parameter potential
of the light chiral condensate, and discuss the regime in the phase diagram of
QCD which exhibits quasi-massless modes.

In Section II we briefly review the functional approach used here. The
definition of the chiral condensate and the determination of the current quark
mass parameters at the physical point are discussed in Section III. There, we
also evaluate the systematic error of the set-up. In Section IV we compute the
magnetic equation of state as well as the mass dependence of the transition
temperature. These results are compared to other functional QCD results as
well as lattice results, also leading to an estimate for the critical
temperature in the chiral limit with a combined systematic error estimate. In
Section V we derive a formula for the order parameter of the light chiral
condensate, and compute it together with simple polynomial fits. Moreover, we
also parameterise the current quark mass dependence of the light chiral
condensate, and use these results in combination to estimate the validity
regime of chiral expansion schemes. In Section VI we summarise our findings.

##  II fRG-assisted Dyson-Schwinger equations

In this Section we review the functional QCD approach used here and put
forward in Gao and Pawlowski (2020, 2021); Gao  _et al._ (2021a). In this
approach the Dyson-Schwinger equations (DSEs) at finite temperature and
density have been expanded about QCD within another parameter set (e.g.
flavour, temperature, density). In Gao and Pawlowski (2020, 2021), the
correlation functions of two-flavour QCD in the vacuum have been used as
input, obtained in the functional renormalisation group (fRG) approach, Cyrol
_et al._ (2018b). In the present work we also utilise 2+1 flavour data from
Gao and Pawlowski (2021) and the recent 2+1 flavour vacuum QCD precision data
from Gao  _et al._ (2021a). The latter data are in particular used as
benchmark tests.

Specifically we use 2 and 2+1 flavour vacuum data of the gluon propagator and
the quark-gluon vertex. This input enables us to compute the quark propagator
in the vacuum from its gap equation, see Figure 1. The quantitative accuracy
of this result is an important self-consistency check of the approach. Then
the gap equation is solved at finite temperature for different current quark
masses. This requires the gluon propagator and quark-gluon vertex at finite
temperature, and we expand the respective DSEs about two-flavour QCD, hence
only solving for the quark mass and temperature dependence, see Figure 2 and
Figure 3.

In the following we briefly recapitulate our setup, more details can be found
in Gao and Pawlowski (2020, 2021); Gao  _et al._ (2021a).

###  II.1 Quark gap equation and chiral phase transition

The quark gap equation relates the inverse quark propagator
$\Gamma^{(2)}_{q\bar{q}}$ to its classical counter part $S^{(2)}_{q\bar{q}}$,
to the quark and gluon propagators, and the classical and full quark-gluon
vertex, see Figure 1. In the vacuum we write

$\displaystyle\Gamma^{(2)}_{q\bar{q}}(p)=Z_{q}(p^{2})\left[\text{i}\,/{p}+M_{q}(p)\right]\,,$
(1)

where

$\displaystyle\Gamma^{(n)}_{\Phi_{i_{1}}\cdots\Phi_{i_{n}}}(p_{1},...,p_{n})=\frac{\delta\Gamma[\Phi]}{\delta\Phi_{i_{1}}(p_{1})\cdots\delta\Phi_{i_{n}}(p_{n})}\,,$
(2)

denotes the 1PI correlation functions of QCD, with
$\phi=(A_{\mu},c,\bar{c},q,\bar{q})$, for more details see e.g. Cyrol  _et
al._ (2018a); Fu  _et al._ (2020b); Gao and Pawlowski (2020, 2021); Gao  _et
al._ (2021a).

At finite temperature the rest frame singles out the temporal direction and we
write,

$\displaystyle\Gamma^{(2)}_{q\bar{q}}(p)=Z^{{\mathrel{\raisebox{1.2pt}{\scriptsize$\parallel$}}}}_{q}(p)\Bigr{[}\text{i}\,\gamma_{0}\,p_{0}+M_{q}(p)\Bigr{]}+Z_{q}(p)\,\text{i}\,\vec{\gamma}\,\vec{p}\,,$
(3)

where $p_{0}$ are now thermal Matsubara frequencies $\omega_{n}$,

$\displaystyle\\!\\!A_{\mu},c,\bar{c}:\ \omega_{n}=2\pi Tn\,,\quad q,\bar{q}:\
\omega_{n}=$ $\displaystyle\,2\pi T\left(n+\frac{1}{2}\right)\,.$ (4)

In eq. 3, $Z^{{\mathrel{\raisebox{1.2pt}{\scriptsize$\parallel$}}}}_{q}$ is
the wave function renormalisation for the mode parallel to the rest frame, and
$Z_{q}=Z^{\bot}_{q}$, is the wave function renormalisation for the transverse
modes perpendicular to the rest frame. With the parametrisation eq. 3,
$M_{q}(p)$ is defined as the pole mass. The transverse modes carry more weight
in the DSE loop integrals, and we shall use the approximation
$Z^{{\mathrel{\raisebox{1.2pt}{\scriptsize$\parallel$}}}}_{q}/Z_{q}\approx 1$
in the vertices. In this approximation the difference between the longitudinal
color-electric and transverse color-magnetic dressings is ignored. In Section
III.3 the respective systematic error is discussed in detail.

Within this approximation the quark DSE at finite temperature reads,

$\displaystyle\Gamma^{(2)}_{q\bar{q}}(\tilde{p})-S^{(2)}_{q\bar{q}}(\tilde{p})=$
$\displaystyle\,\int\sum\frac{dq_{0}}{2\pi}\\!\int\frac{d^{3}{q}}{(2\pi)^{3}}\;\Biggl{[}G_{AA}{}_{\mu\nu}^{ab}(q+p)\quad$
$\displaystyle\times\frac{\lambda^{a}}{2}{(-ig\gamma_{\mu})}G_{q\bar{q}}({q})\left[\Gamma^{(3)}_{q\bar{q}A}\right]^{b}_{\nu}(q,-p)\Biggr{]}\,.$
(5)

In Equation 5, $S^{(2)}_{q\bar{q}}$ is the inverse of classical quark
propagator, the full quark propagator
$G_{q\bar{q}}=(1/\Gamma^{(2)})_{q\bar{q}}$ reads

$\displaystyle
G_{q\bar{q}}(q)=-\frac{Z^{{\mathrel{\raisebox{1.2pt}{\scriptsize$\parallel$}}}}_{q}(p)\Bigr{[}\text{i}\,\gamma_{0}\,p_{0}-M_{q}(p)\Bigr{]}+Z_{q}(p)\,\text{i}\,\vec{\gamma}\,\vec{p}}{\left[Z^{{\mathrel{\raisebox{1.2pt}{\scriptsize$\parallel$}}}}_{q}(p)\right]^{2}\Bigl{[}p_{0}^{2}+M_{q}^{2}(p)\Bigr{]}+Z^{2}_{q}(p)\,\vec{p}^{2}}\,.$
(6)

Finally, the gluon propagator is defined by $G_{AA}=(1/\Gamma^{(2)})_{AA}$,
and is detailed below, see eq. 8 and eq. 8. All momenta are counted as
incoming.

Figure 1: Quark gap equation. Lines with a blob are full propagators, that
without is the classical quark propagator. The Vertex with blob is the full
quark-gluon vertex, that without is the classical one.

Equation 5 can be solved with the knowledge of the gluon propagator $G_{AA}$
and the quark-gluon vertex $\Gamma_{q\bar{q}A}^{(3)}$. In Section II.2 and
Section II.3 we detail how these correlation functions at finite temperature
are computed on the basis of the vacuum two flavour fRG-data from Cyrol  _et
al._ (2018a).

###  II.2 DSE for the gluon propagator

The DSE for the inverse gluon propagator for general flavours $N_{f}$ at
finite temperature is expanded about that in the vacuum for two-flavour QCD.
Schematically this reads,

$\displaystyle\left.\Gamma^{(2)}_{AA}(p)\right|_{T,N_{f}}=\left.\Gamma^{(2)}_{AA}(p)\right|_{0,2}+\Delta\Gamma^{(2)}_{AA}(p)\,,$
(7)

and leaves us with a DSE for the difference $\Delta\Gamma^{(2)}_{AA}$ between
the full inverse gluon propagator and that of vacuum $2$-flavour QCD, depicted
in Figure 2.

At finite temperature, the gluon two-point function has color-electric and
color-magnetic components,

$\displaystyle\Gamma^{(2)}_{AA}(p)=p^{2}\left[Z^{M}_{A}(p)\Pi^{M}_{\mu\nu}(p)+Z^{E}_{A}(p)\Pi^{E}_{\mu\nu}(p)\right]+\frac{p_{\mu}p_{\nu}}{\xi}\,,$
(8)

with the color-magnetic dressing $Z^{M}_{A}$ and the color-electric one,
$Z^{E}_{A}$. In the present work we use the Landau gauge, $\xi\to 0$. The
projection operators onto the color-electric and color-magnetic directions in
eq. 8 read,

$\displaystyle\Pi^{M}_{\mu\nu}(p)=(1-\delta_{0\mu})(1-\delta_{0\nu})\left(\delta_{\mu\nu}-\frac{p_{\mu}p_{\nu}}{\vec{p}^{2}}\right)\,,$
$\displaystyle\Pi^{E}_{\mu\nu}(p)=\delta_{\mu\nu}-\frac{p_{\mu}p_{\nu}}{{p}^{2}}-\Pi^{M}_{\mu\nu}(p)\,.$
(9)

As discussed in detail in Gao and Pawlowski (2020, 2021), such difference DSEs
are very stable and converge very quickly in an iterative procedure about the
initial value $\Delta\Gamma^{(n)}=0$. This can be traced back to the
relatively small size of the thermal correction of the gluon propagator and
quark-gluon vertex for the temperatures of interest, see e.g. Fu  _et al._
(2020b); Braun  _et al._ (2020a); Fischer (2019); Gunkel  _et al._ (2019); Gao
and Pawlowski (2020, 2021). The small size of the corrections also support a
further approximation of the difference DSE for the gluon propagator with

$\displaystyle
T\sum_{\omega_{n}}\textrm{loop}_{T}(q,p)-\int\frac{d\omega}{2\pi}\textrm{loop}_{\textrm{vac}}(q,p)$
$\displaystyle=$
$\displaystyle\left[T\sum_{\omega_{n}}\textrm{loop}_{\textrm{vac}}(q,p)-\int\frac{d\omega}{2\pi}\textrm{loop}_{\textrm{vac}}(q,p)\right]$
$\displaystyle+T\sum_{\omega_{n}}\Biggl{[}\textrm{loop}_{T}(q,p)-\textrm{loop}_{\textrm{vac}}(q,p)\Biggr{]}$
$\displaystyle\approx$
$\displaystyle\left[T\sum_{\omega_{n}}\textrm{loop}_{\textrm{vac}}(q,p)-\int\frac{d\omega}{2\pi}\textrm{loop}_{\textrm{vac}}(q,p)\right]\,,$
(10)

with $q_{0}=\omega_{n}$ at $T\neq 0$ and $q_{0}=\omega$ at $T=0$, and
$\textrm{loop}(q,p)$ stands for the loops in the second line of Figure 2. The
third line vanishes for $\Delta\Gamma^{(2)}_{AA}=0$, and is negligible for
small $\Delta\Gamma^{(2)}_{AA}$. Accordingly we have dropped it, but have we
have monitored its irrelevance in our explicit computation. The results for
the gluon propagator in this approximation agree well with respective finite
temperature results, see Gao and Pawlowski (2020); Gao  _et al._ (2021a); Fu
_et al._ (2020b), for more details we refer to these works.

Figure 2: Gluon DSE for the difference $\Delta\Gamma_{AA}^{(2)}$ between the
full gluon propagator and the vacuum $2+1$-flavour gluon propagator. Lines and
vertices with black blobs are full $N_{f}=2+1$ propagators and vertices at
finite temperature. Lines and vertices with grey blobs are full vacuum
propagators and vertices for $N_{f}=2$. The square bracket contains the
temperature fluctuations and does not require renormalisation. We have dropped
the part of the classical gluon propagator that carries the renormalisation of
the strange loop.

Moreover, as in Gao and Pawlowski (2020, 2021) we drop the difference between
color-magnetic and color-electric dressings: only the color-magnetic dressing
$Z^{E}_{A}(p)$ is computed, and the color-electric one is approximated as
$Z^{E}_{A}(p)\approx Z^{M}_{A}(p)$. This is the same approximation also used
for the quark propagator, and the respective systematic error or rather is
smallness is discussed in Section III.3.

###  II.3 Quark-gluon vertex at finite $T$

The DSE for the quark-gluon vertex for $2+1$-flavour QCD at finite temperature
is expanded about its two- or 2+1 flavour counter part in the vacuum, computed
in Cyrol  _et al._ (2018a) (2 flavour) and Fu  _et al._ (2020b); Gao and
Pawlowski (2020, 2021); Gao  _et al._ (2021a) (2+1 flavour), see also Williams
(2015); Williams  _et al._ (2016); Aguilar  _et al._ (2017, 2018b). We have
checked, that the results of expansions about the two and 2+1 flavour vacuum
results agree quantitatively. The expansion about 2+1 flavour vacuum QCD
reads,

$\displaystyle\left.\Gamma^{(3)}_{q\bar{q}A}(p_{1},p_{2})\right|_{T,N_{f}}\\!\\!=\left.\Gamma^{(3)}_{q\bar{q}A}(p_{1},p_{2})\right|_{0,N_{f}}\\!\\!+\Delta\Gamma^{(3)}_{q\bar{q}A}(p_{1},p_{2})\,.$
(11)

for the light quarks, $q=l$, and the strange quark, $q=s$.

At finite temperature the vertex has color-electric and magnetic components.
We identify them within our $O(4)$-symmetric approximation, as the splitting
is weak and we only keep the vacuum tensor structures. A complete basis in the
vacuum with twelve basis elements is given by the eight transverse and four
longitudinal projections of the Lorentz tensors

$\displaystyle\begin{array}[]{lcl}\left[{\cal
T}^{(1)}_{q\bar{q}A}\right]^{\mu}(p,q)=-i\gamma^{\mu}\,,&&\left[{\cal
T}^{(5)}_{q\bar{q}A}\right]^{\mu}(p,q)=i{k\hbox
to0.0pt{\hss/}}_{+}k_{-}^{\mu}\,,\\\\[8.61108pt] \left[{\cal
T}^{(2)}_{q\bar{q}A}\right]^{\mu}(p,q)=k_{-}^{\mu}\,,&&\left[{\cal
T}^{(6)}_{q\bar{q}A}\right]^{\mu}(p,q)=ik\hbox
to0.0pt{\hss/}_{-}k_{-}^{\mu}\,,\\\\[8.61108pt] \left[{\cal
T}^{(3)}_{q\bar{q}A}\right]^{\mu}(p,q)={k\hbox
to0.0pt{\hss/}_{-}}\gamma^{\mu}\,,&&\left[{\cal
T}^{(7)}_{q\bar{q}A}\right]^{\mu}(p,q)=\frac{i}{2}[p\hbox
to0.0pt{\hss/},q\hbox to0.0pt{\hss/}]\gamma^{\mu}\,,\\\\[8.61108pt]
\left[{\cal T}^{(4)}_{q\bar{q}A}\right]^{\mu}(p,q)=k\hbox
to0.0pt{\hss/}_{+}\gamma^{\mu}\,,&&\left[{\cal
T}^{(8)}_{q\bar{q}A}\right]^{\mu}(p,q)=-\frac{1}{2}[p\hbox
to0.0pt{\hss/},q\hbox to0.0pt{\hss/}]k_{-}^{\mu}\,,\end{array}$ (16)

multiplied by the Gell-Mann matrices $\lambda^{a}$, see e.g. Cyrol  _et al._
(2018a); Gao and Pawlowski (2020). The transverse projection operator is
denoted by $\Pi^{\bot}_{\mu\nu}(k)=\,\delta_{\mu\nu}-k_{\mu}k_{\nu}/k^{2}$,
and the longitudinal one is given by $\Pi^{\parallel}=\mathbbm{1}-\Pi^{\bot}$.
In the Landau gauge only the transverse part of the vertex enters the
(transverse) DSEs, and hence suffices to access the full dynamics of QCD. With
eq. 16 the full transverse quark-gluon vertex is given by,

$\displaystyle\left[\Gamma^{(3)}_{q\bar{q}A}\\!\right]^{a}_{\mu}(p,q)=\frac{\lambda^{a}}{2}\sum_{i=1}^{8}\lambda^{(i)}_{q\bar{q}A}(p,q)\Pi^{\bot}_{\mu\nu}(k_{+})\left[{\cal
T}^{(i)}_{q\bar{q}A}\\!\right]_{\nu}\\!(p,q)\,.$ (17)

Figure 3: Quark-gluon DSE for the difference $\Delta\Gamma_{q\bar{q}A}^{(3)}$
between the full quark-gluon vertex (with black blob) and the vacuum
$2$-flavour quark-gluon vertex (with grey blob). The temperature fluctuations
in the square bracket do not require renormalisation, as in the gluon DSE
Figure 2.

Equation 11 leaves us with the task of solving the Dyson-Schwinger equation
for the difference. In the present work we take into account all tensor
structures in the diagrams in the DSE for the quark-gluon vertex. We derive
the DSE by considering the functional DSE for $\delta\Gamma/\delta A_{\mu}$
and then taking $q,\bar{q}$-derivatives. Within this hierarchy we have only
consider the diagrams, that have a perturbative one-loop counterpart, and the
resulting DSE is depicted in Figure 3. The diagrams dropped can be partly
understood or taken into account as vertex dressings of the bare vertices in
Figure 3.

In our opinion, the most prominent omission in the present approximation is
the one-loop diagram with a four-quark vertex that carries (off-shell) two-
quark resonances such as mesons and diquarks. Note that this omission only
concerns the thermal part of the vertex, in the fRG computation underlying the
input data in two-flavour vacuum QCD, Cyrol  _et al._ (2018a), the Fierz-
complete four-quark interaction vertex has been taken into account. In any
case, the omission of the thermal correction of the dominant scalar-
pseudoscalar channel ($\sigma$-pion) potentially leads to an underestimation
of the infrared dynamics of QCD and in a linear estimate it enhances chiral
symmetry breaking and hence the critical temperatures. Moreover, these
contributions are relevant for the emergence of critical scaling including the
values of the critical exponents in the chiral limit, see in particular
Fischer and Mueller (2011). Their inclusion is subject of ongoing work.

##  III Physical point and error estimates

In Section III.1 we define our order parameter for the chiral phase
transition, the renormalised chiral condensate. We use the physical point in
the vacuum for setting our scales, see Section III.2, and discuss estimates of
the systematic error of the magnetic equation of state and deduced observables
such as the chiral transition temperature in Section III.3.

###  III.1 Chiral transition temperature and the renormalised chiral
condensate

The chiral transition temperature is determined by the peak of the
susceptibility of the renormalised light chiral condensate $\Delta_{l,R}$,
with $l=u,d$ in the present isospin-symmetric approximation, defined by the
thermal part of the light chiral condensate $\Delta_{l}$,

$\displaystyle\Delta_{l,R}=\frac{1}{{\cal
N}_{R}}\Bigl{[}\Delta_{l}(T)-\Delta_{l}(0,0)\Bigr{]}\,.$ (18a) with
$\displaystyle\Delta_{l}\simeq-\frac{1}{2}\,T\sum_{n\in\mathbb{Z}}\int\frac{d^{3}q}{(2\pi)^{3}}{\text{tr}}\,G_{l\bar{l}}(q)\,,$
(18b)

where the factor 1/2 in eq. 18b cancels that from the sum over the light
flavours, and $G_{l\bar{l}}$ is the light quark propagator, see eq. 6. We
remark that the definition of the chiral condensate in eq. 18b lacks the
further factor $m_{l}$ used in Braun  _et al._ (2020a). While such a
definition ensures RG-invariance of the condensate, the present definition is
that used also in chiral perturbation theory. Moreover, eq. 18b has to be
renormalised, for a respective detailed discussion within the present
approximation see Gao  _et al._ (2021a). For the present work the details of
the renormalisation are not relevant: it drops out in eq. 18a due to the
subtraction, which renders $\Delta_{l,R}$ finite. The factor ${\cal N}_{R}$ in
eq. 18a is a convenient normalization which leaves $\Delta_{l,R}$
dimensionless, and we choose ${\cal N}_{R}=m^{4}_{\pi}$.

###  III.2 Physical point

In Gao  _et al._ (2021a) it has been shown, that the current approximation
leads to a chiral condensate in quantitative agreement with recent lattice
estimates, see e.g. the recent FLAG-review Aoki  _et al._ (2020). We expect
that this accuracy propagates to the renormalised chiral condensate, built on
the same approximation scheme. Note however, that the difference definition
eliminates the explicit dependence of $\Delta_{l,R}$ on the vacuum condensate.

For the estimate of the pion mass we apply the one-loop formula from chiral
perturbation theory ($\chi$PT) Gasser and Leutwyler (1984) (valid in
$2$-flavour QCD),

$\displaystyle m^{2}_{\pi}=$
$\displaystyle\,M^{2}\left(1-\frac{x}{2}\,\mathrm{Log}\frac{\Lambda_{3}^{2}}{M^{2}}+{\cal
O}(x^{2})\right)\,,$ $\displaystyle f_{\pi}=$
$\displaystyle\,f_{\pi,\chi}\left(1+x\,\mathrm{Log}\frac{\Lambda_{4}^{2}}{M^{2}}+{\cal
O}(x^{2})\right)\,,$ (19a) with the expansion parameter $x$, that depends on
the mass parameter $M$. The latter is the tree-level value of the pion mass,
given by the Gell-Mann–Oakes–Renner (GMOR) relation, see e.g. Bender  _et al._
(1998); Gao and Liu (2018), for reviews see Fischer (2006); Bashir  _et al._
(2012); Eichmann  _et al._ (2016b). We have, $\displaystyle
M^{2}=\,\frac{2\,\Delta_{l,\chi}}{f^{2}_{\pi,\chi}}\,m_{l}\,,\qquad
x=\frac{M^{2}}{(4\pi f_{\pi,\chi})^{2}}\,.$ (19b) Equation 19a also depends on
the low energy constants
$\displaystyle\Lambda_{3}=640(19)\,\textrm{MeV}\,,\quad\Lambda_{4}=1030(21)\,\textrm{MeV}\,,$
(19c) that can be computed in a Bethe-Salpeter–DSE approach. Here we simply
use the recent FLAG estimates in Aoki  _et al._ (2020), for a discussion of
the dependence on the RG-scheme see Gao  _et al._ (2021a). There, the
$\Lambda_{n}$ follow from
$\bar{l}_{n}=\log{\Lambda_{n}}/{m_{\pi,\textrm{phys}}}$ and
$\bar{l}_{3}=3.07(64)$, and $\bar{l}_{4}=4.02(45)$ for $2+1$ flavour QCD. For
sufficiently light quark masses $m_{l}$ considered here the ${\cal
O}(m_{l}^{3})$ is known to be small both on the lattice, Aoki  _et al._
(2020), and in quantitative Bethe-Salpeter computations, see e.g. Eichmann
_et al._ (2016b); Chen  _et al._ (2018). For the decay constant, we utilise
the FLAG results in Aoki  _et al._ (2020) for the pion decay constants in the
chiral limit, $\displaystyle f_{\pi,\chi}=86.7\,\textrm{MeV}\,,$ (19d)

in the isospin symmetric limit. The vacuum quark condensate $\Delta_{l,\chi}$
can be extracted from the UV behaviour of the corresponding constituent quark
mass, see Gao  _et al._ (2021a). In the vacuum, the current approximation
reduces to the the same quantitatively reliable approximation used in Gao  _et
al._ (2021a). However, in the present work we expand the system of equations
about the fRG vacuum solution for two-flavour QCD, while in Gao  _et al._
(2021a) the DSE for the quark-gluon vertex was solved. Accordingly, the
present benchmark results in the vacuum are in quantitative agreement with
that in Gao  _et al._ (2021a), but not identical. Note that this remarkable
agreement is yet another self-consistency check of functional methods. In the
chiral limit we find,

$\displaystyle\Delta_{l,\chi}(\mu)=\,(308.9(8)\,\textrm{MeV})^{3}\,,$ (20)

in quantitative agreement with the lattice estimate
$\Delta_{l,\chi}(\mu)=\,315(6)\,\textrm{MeV})^{3}$ and the DSE results in Gao
_et al._ (2021a). The above relations allow us to fix the physical point: the
current quark masses $m_{q}(\mu)$ with an RG scale of $\mu=40$ GeV are fixed
such that we have in the isospin symmetric limit with $m_{l}=m_{u}=m_{d}$,

$\displaystyle
m_{\pi,\textrm{phys}}=138\,\textrm{MeV}\,,\qquad\frac{m_{s,\textrm{phys}}}{m_{l,\textrm{phys}}}=27\,.$
(21a) This leads us to $\displaystyle
m_{l,\textrm{phys}}(\mu)=2.49\,\textrm{MeV}\,,\qquad
m_{s,\textrm{phys}}(\mu)=66.7\,\textrm{MeV}\,,$ (21b) with the prediction
$\displaystyle f_{\pi,\textrm{phys}}=92.4\,\textrm{MeV}\,,$ (21c)

in line with the physical value and the lattice estimate in Aoki  _et al._
(2020). We emphasise that this perfect agreement originates in the
quantitative agreement of the chiral condensate in the chiral limit,
$\Delta_{l,\chi}(\mu)$ in eq. 20, from the present DSE computation with the
respective lattice estimate. As discussed in Gao  _et al._ (2021a),
$\Delta_{l,\chi}(\mu)$ is the benchmark observable in the vacuum which shows
the quantitative reliability of the present approximation of the DSE and fRG
in Cyrol  _et al._ (2018a). Additionally, the correct value for
$f_{\pi,\textrm{phys}}$ is a test of the 2-flavour approximations used in eq.
19a.

###  III.3 Systematic error estimate

The present functional approach, or rather the approximation level of the
systematic vertex expansion used here, has passed successfully numerous
benchmark comparisons, a systematic comparison can be found in Gao  _et al._
(2021a), a discussion of the systematic error estimate is done in Section
III.3. Apart from the impressive quantitative agreement of benchmark
observables such as the chiral condensate, results of different functional
approaches that constitute different resummation schemes also agree
quantitatively. In conclusion we consider the approximation level used here as
sufficiently converged in the vacuum.

In turn, at finite temperature we use additional approximations. In the
following we discuss the ensuing systematic error of the present approximation
level. The most relevant sources for the systematic error are,

  * (i)

No feedback of the order parameter potential: feeding back the order parameter
potential into the diagrams encodes the critical dynamics beyond mean field,
and we expect deviations from the full result in the scaling regime for
$m_{l}\to 0$.

  * (ii)

No thermal splits in the vertices, most prominently in the quark-gluon vertex:
we identify electric and magnetic dressings and compute this uniform dressing
by magnetic projections.

Let us first discuss the impact of (i): The quantitative access to the
critical dynamics for $m_{l}\to 0$, and hence $m_{\pi}\to 0$, requires the
inclusion of a dynamical order parameter potential. Heuristically speaking the
latter takes into account the multi-quark scatterings of the resonant pion
channel that become relevant in the limit of infinite correlation length. This
approximation is well-tested in the vacuum, and the neglection of the chiral
dynamics encoded in the higher multi-quark scattering vertices is negligible
there. In the quantitative fRG-studies Mitter  _et al._ (2015); Cyrol  _et
al._ (2018a), the higher-order scatterings have been taken into account
systematically via dynamical hadronisation in terms of the scattering of the
respective resonant interaction channels, for a recent review see Dupuis  _et
al._ (2021). This allows a quantitative access to the critical chiral
dynamics, and the effects on vertices and propagators are very small: The
respective changes are well within the systematic error estimates for the
given approximations. Moreover, the results for quark and gluon propagators as
well as the quark-gluon vertex with the full order parameter potential are in
quantitative agreement with that of the DSE-approach used here, see Gao  _et
al._ (2021a). In the context of the present work this is most impressively
confirmed by the chiral condensate in the chiral limit, $\Delta_{l,\chi}$ in
eq. 20, which is in quantitative agreement with other benchmark computations.
In conclusion, the critical dynamics does not affect the results in the
vacuum, and the scaling regime for $m_{l}\to 0$ is very small.

At finite temperature these findings strongly suggest that the higher order
multi-quark scatterings only play a r$\hat{\textrm{o}}$le in the critical
scaling regime for temperatures $T$ close to the chiral phase transition
temperature $T_{c0}$ in the chiral limit and small light current quark masses
$m_{l}\to 0$. These higher scattering processes have been taken into account
in Braun  _et al._ (2020a) for pion masses $m_{\pi}\geq 30$ MeV without any
sign of criticality. Taking also into account the small difference between
mean field scaling and full $O(4)$ or $O(2)$ scaling, we estimate that the
systematic error of neglecting the higher order multi-quark scatterings is
negligible also for $m_{\pi}<30$ MeV. Note that this estimate applies to the
absolute value of observables such as $T_{c}$ as a function of the light quark
mass. As already emphasised at the end of Section II.3, the approximation (i)
has to be improved for the access to critical scaling, and in particular for a
computation of the critical exponents, see Fischer and Mueller (2011). This is
subject of ongoing work.

In summary, we are led to a systematic error estimate for (i), which is
subleading for the observables studied here, and is dominated by other
approximation effects.

We proceed with the evaluation of (ii), the missing thermal splits in
propagators and vertices. This affects the contribution of the zeroth
Matsubara frequency $\omega_{0}$, that is either $\omega_{0}=0$ for the gluons
and ghost, and $\omega_{0}=\pi T$ for the quarks. There, the spatial momentum
dependence is changed for $\vec{q}{\,}^{2}\lesssim(2\pi T)^{2}$. In turn, the
thermal effects for higher Matsubara frequencies are negligible. Within
diagrams, this thermal change of correlation functions is suppressed with the
measure factor $\vec{q}{\,}^{2}$ of the spatial momentum integration. This
combination of the decay of thermal contributions for
$\vec{q}{\,}^{2}\to\infty$ and its phase space suppression for
$\vec{q}{\,}^{2}\to 0$, suggests the feasibility of $O(4)$-symmetric
approximations. Here, $O(4)$ refers to the Euclidean space-time symmetry in
the vacuum. In the present work we implement this idea with using projections
on the magnetic (spatial) dressings of propagators and vertices. Then, the
electric dressing is identified with the magnetic one. The above line of
reasoning has been checked within different theories within the functional
renormalisation group approach, Dupuis  _et al._ (2021), including Yang-Mills
theory at finite temperature, Cyrol  _et al._ (2018b). Finally, we have
extended the approximation used in the present work with the thermal splits
for vertices and propagators, which will be discussed in detail in a
forthcoming work. The respective results corroborate the arguments above and
lead to negligible modifications of the observables including the chiral
condensate. Moreover, the crossover temperature (for physical quark masses)
only increases by less than 3 MeV, Gao  _et al._ (2021b).

Finally, we discuss the systematic error estimate for such an approximation.
To begin with, the above suppression of thermal contributions for
$\vec{q}{\,}^{2}\lesssim(2\pi T)^{2}$ is undone for correlation functions,
whose functional relations(fRG or DSEs) have sizeable infrared contributions
for $\vec{q}{\,}^{2}/m^{2}_{\textrm{gap}}\to 0$. Here, $m^{2}_{\textrm{gap}}$
is the characteristic mass scale of QCD: below this mass scale, momentum
fluctuations do not contribute significantly to the physics. This entails that
$m^{2}_{\textrm{gap}}$ is triggered by the gluon mass gap as well as he
constituent quark mass. Note however, that such a survey also has to take into
account effective degrees of freedom, and in particular the pion. It is here,
where (i) comes into play. Again, from detailed studies in the vacuum we know
that the characteristic scale for fluctuations for the quark mass function and
consequently chiral condensate is $m_{\textrm{gap}}\approx 1$ GeV. This
structural analysis is well confirmed within fRG-studies, where momentum
fluctuations are resolved iteratively momentum-shell by momentum-shell for
$\vec{q}{\,}^{2}\approx k^{2}$, the latter being the infrared cutoff scale:
the infrared decay of momentum fluctuations is clearly visible. The above
argument also entails, that a phenomenological infrared enhancement of
vertices, as often used in functional studies, has to be taken with a grain of
salt as it may over-enhance thermal (and chemical potential) effects.

Apart from these structural arguments the present setup offers several simple
self-consistency checks: Firstly, instead of using the magnetic projection for
the computation of the $O(4)$-symmetric dressing we can either average over
all Euclidean directions or we even take the electric projection. We have
checked that this only leads to variations of the chiral transition
temperature of a few MeV, indeed safely below $5$ MeV. Secondly, the slope or
curvature of the chiral transition temperature $T_{c}(\mu_{B})$ at vanishing
baryon chemical potential, $\mu_{B}=0$, is a direct check of the correct
infrared dynamics at finite temperature and density: despite the differences,
i.e. the lack of direct density contributions to gluons and further
correlations with vanishing baryon number, thermal and density effects share
the same tendencies. In particular, increasing the infrared vertex strength of
the quark-gluon vertex increases the size of thermal and chemical potential
effects. Within the present approximation level of the DSE approach to QCD
that underlies the current study, it is quantitatively agreeing with the
lattice prediction, Gao and Pawlowski (2021, 2020). This also holds true for
the related fRG study in Fu  _et al._ (2020b), which underlies the
investigation of the magnetic equation of state in Braun  _et al._ (2020a).

In summary, the Euclidean $O(4)$-symmetric approximation used here is sound,
and the systematic error estimate is achieved by applying an artificial
infrared enhancement to the quark-gluon vertex that leads to a variation of
the chiral transition temperature $T_{c}$ at physical quark masses of $5$ MeV.
The respective systematic error band is then shown in our depiction of the
results for the chiral transition temperature, see Figure 5 and Figure 6.

##  IV Magnetic equation of state

A first main result is given by the chiral transition temperature
$T_{c}(m_{l})$ in dependence of the ratio $H=m_{l}/m_{s}$ of current quark
masses, that can be derived from the magnetic equation of state: the
transition temperature is obtained from the peak of the susceptibility of the
renormalised light chiral condensate defined in Section III.1. The results are
discussed in Section IV.1, including a comparison with the respective lattice
results, Ding  _et al._ (2019); Kaczmarek  _et al._ (2020a, b); Karsch (2019);
Kotov  _et al._ (2021), as well as the functional QCD study (fRG) in Braun
_et al._ (2020a). For the comparison with the latter results we also present
$T_{c}(m_{\pi})$. This requires $m_{\pi}(m_{l})$, for which we use eq. 19a,
valid at one-loop up to $O(x^{3})$. The latter higher order terms are known to
be small for the pion masses $m_{\pi}\leq 140$ MeV considered in Braun  _et
al._ (2020a). In Section IV.2 we discuss constraints on the equation of states
as well as the merits of future combined studies, utilising both, functional
results as well as lattice simulations.

###  IV.1 Transition temperature and transition line

The chiral transition temperature is computed from the peak of the
renormalised light quark susceptibility $\chi^{(l,R)}_{M}$, given by,

$\displaystyle\chi_{M}^{\
}=\chi^{(l,R)}_{M}=-\frac{\partial\Delta_{l,R}}{\partial m_{l}}\,.$ (22)

With eq. 22 we define the chiral transition temperature as the peak position
of $\chi^{(l,R)}_{M}$,

$\displaystyle
T_{c}=T^{(l,R)}_{c}=T_{\textrm{peak}}\,,\quad\quad\left.\frac{\partial\chi^{(l,R)}_{M}}{\partial
T}\right|_{T=T_{\textrm{peak}}}=0\,.$ (23)

In Braun  _et al._ (2020a), also the magnetic susceptibilities for the light
chiral condensate $\chi^{(l)}_{M}$ and the reduced chiral condensate
$\chi^{(l,s)}_{M}$ have been considered. There, it has been shown, that the
respective transition temperatures $T_{c}^{(i)}$ with
$(i)=(l)\,,\,(l,R)\,,\,(l,s)$ agree well in the crossover regime. Evidently,
$T_{c}^{(l)}=T_{c}^{(l,R)}$. The second approximate relation,
$T_{c}^{(l,s)}\approx T_{c}^{(l,R)}$ follows from the fact, that the
transition temperatures only vary within the width of the crossover in the
first place, and the thermal dependence of the strange quark condensate is far
reduced. Moreover, all transition temperatures agree trivially in the second
order regime.

Figure 4: Temperature dependence of the renormalised magnetic susceptibility
$\chi^{(l,R)}_{M}$, eq. 22, for various pion masses. The results are compared
to the respective fRG-results Braun  _et al._ (2020a).

We have benchmarked our scales in the vacuum with the physical point with a
isospin-symmetric pion mass, $m_{\pi,\textrm{phys}}=138$ MeV, see eq. 21. The
chiral transition temperature at the physical point now follows as

$\displaystyle T_{c}=154.7^{+5}_{-5}\,\textrm{MeV}\,,$ (24)

with the systematic error estimate (ii) discussed in Section III.3. It encodes
the systematic error of the strength of the quark-gluon coupling. The latter
has been varied such that $T_{c}(m_{\pi,\textrm{phys}})$ for physical pion
masses with $\pm$5 MeV.

In Figure 4 we depict the chiral susceptibilities for $m_{\pi}=30,55,80,140$
MeV in comparison to the fRG results in Braun  _et al._ (2020a). The
respective transition temperatures are depicted in Figure 5 as a function of
$H=m_{l}/m_{s}$. One clearly sees the $\sqrt{H}$-behaviour. Note that
different definitions of the transition temperature lead to minor changes of
the location of phase transition in the crossover regime, but they all agree
in the second order regime in the chiral limit. We also depict the pion mass,
obtained in the approximation eq. 19a, as a function of current quark mass in
Figure 5.

In the present approximation to the $m_{l}$-dependence of the pion mass, eq.
19a, evidently $m_{\pi}\propto\sqrt{H}$ holds close to the chiral limit. In
turn, a linear dependence $m_{\pi}\propto H$ emerges for asymptotically large
light current quark masses. Together with the approximate, but quantitative,
$\sqrt{H}$-behaviour of the chiral transition temperature, this leads us to
the linear dependence $T_{c}\propto m_{\pi}$ close to the chiral limit. The
linear dependence is clearly visible in Figure 6, and holds true even up to
$H=1$, that is $m_{l}=m_{s}$.

In the scaling regime of QCD close to the chiral limit, the chiral phase
transition temperature should show critical scaling with $T_{c}\sim
H^{1/(\beta\delta)}$. Here, the critical exponents $\beta,\delta$ determine
the universality class of QCD in the chiral limit: $m_{s}$ fixed, $H\to 0$.

A detailed discussion in view of the universality class of QCD in the chiral
limit of the functional QCD results in the present work and in Braun  _et al._
(2020a), as well as the lattice results of Ding  _et al._ (2019); Kotov  _et
al._ (2021) is given in Section V. For the time being, we just remark that the
approximation to functional QCD with the fRG in Braun  _et al._ (2020a)
incorporates criticality in terms of the dynamics of the order parameter and
the order parameter potential. This has been tested both in generic low energy
effective theories and in QCD, see Dupuis  _et al._ (2021). Consequently, if
critical scaling is present in QCD, it is well-captured within the
approximation to QCD used in Braun  _et al._ (2020a), leading to non-trivial
$\beta,\delta$’s. Note also that the extraction of non-trivial scaling is
technically not challenging in the fRG approach. In short, non-trivial scaling
including sub-leading contributions can hardly be overlooked in fRG
computations.

Figure 5: Chiral transition temperature for $N_{f}=2+1$ flavours as a function
of $H=m_{l}/m_{s}$ in comparison to the fRG-results from Braun  _et al._
(2020a). The crossover temperature $T_{c}$ is defined by eq. 23. The inset
shows $T_{c}$ as a function of $\sqrt{H}$. Both curves are well approximated
by $\sqrt{H}$ for $H\lesssim 0.25$. The blue error band is a systematic error
estimate on the DSE results: it is obtained by varying the overall strength of
the quark-gluon coupling such that $T_{c}(m_{\pi})$ changes in a range of
$T_{c}=155\pm 5$ MeV. Figure 6: Chiral transition temperature for $N_{f}=2+1$
flavours as a function of the pion mass in comparison to the fRG-results from
Braun  _et al._ (2020a) . $T_{c}\propto m_{\pi}$ is a very good approximation
for $H\lesssim 0.25$. The inset shows the pion mass as a function of
$H=m_{l}/m_{s}$ within $\chi$PT, see eq. 19. The blue error band is a
systematic error estimate on the DSE results: it is obtained by varying the
overall strength of the quark-gluon coupling such that $T_{c}(m_{\pi})$
changes in a range of $T_{c}=155\pm 5$ MeV.

In Braun  _et al._ (2020a) it has been shown, that critical scaling is absent
(in the approximation used there) at least for pion masses $m_{\pi}\gtrsim 30$
MeV. Hence, the results in Braun  _et al._ (2020a) do not suggest a large
scaling regime. This corroborate previous detailed scaling analyses including
the study of finite volume scaling in advanced low energy effective theories
of QCD within the past two decades, see in particular Braun  _et al._ (2011a)
and the review Klein (2017).

In conclusion, the above mentioned intricacy of the determination of the
universality class is of subleading importance for the extraction of the
chiral transition temperature in the chiral limit. Owing to the small critical
regime and the proximity of the critical scaling to the mean-field one, we
arrive at a critical temperature in the chiral limit of

$\displaystyle T_{c0}=141.3^{+3.0}_{-3.3}\,\textrm{MeV}\,.$ (25)

The error estimate in eq. 25 is derived from the systematic error estimate
(ii), leading to an error of $\pm$5 MeV for the physical point, see eq. 24.
The critical temperature in the chiral limit, eq. 25 agrees quantitatively
with the fRG-estimate in Braun  _et al._ (2020a), while it is consistent with
the lattice estimates in Ding  _et al._ (2019) with $T_{c0}=132^{+3}_{-6}$
MeV, and Kotov  _et al._ (2021) with $T_{c0}=134^{+6}_{-4}$ MeV. We emphasise
that the coincidence of the chiral transition temperature in the chiral limit
as well as for all light current quark masses investigated with the fRG
results is rather non-trivial in view of the different approximations
involved: while the present work includes, in contradistinction to Braun  _et
al._ (2020a), the full quark-gluon vertex for all temperatures, the fRG-
computation in Braun  _et al._ (2020a) includes multi-scattering vertices of
the scalar-pseudoscalar four-quark channel ($\sigma$ and pion). The latter
allow for a quantitative access to criticality and the size of the scaling
regime, which is missing in the present work. In combination, a consistent
picture emerges: the scaling regime is rather small, and the largest
systematic error left within functional approaches concerns the slope of the
chiral transition temperature $T_{c}(H)$.

|  $D(\pi)$ |  $m_{\pi}$ [MeV]
---|---|---
|  |  0 |  $30$ |  $40$ |  $55$ |  $70$ |  $80$ |  $100$ |  $110$ |  $120$ |  $140$ |  $210$ |  $370$
$T_{c}$[MeV] | fRG-DSE |  0.101 |  141.3 |  $144.2$ |  $145.3$ |  $146.5$ |  $148.3$ |  $149.1$ |  $151.3$ |  $152.1$ |  $153.2$ |  $155.4$ |  $161.4$ |  $177.4$
fRG |  0.10 |  – |  $145.3$ |  $146.4$ |  $148.0$ |  $149.6$ |  $150.5$ |  $152.7$ |  $153.6$ |  $154.8$ |  $156.3$ |  – |  –
hotQCD: $N_{\tau}$ = 8 |  0.106 |  – |  – |  – |  $150.9(4)$ |  – |  $153.9(3)$ |  – |  $157.9(3)$ |  – |  $161.0(1)$ |  – |  –
hotQCD: $N_{\tau}$ = 12 |  – |  – |  – |  – |  – |  – |  $149.7(3)$ |  – |  $155.6(6)$ |  – |  $158.2(5)$ |  – |  –
hotQCD: Cont. Ex. |  0.189 |  132 |  – |  – |  – |  – |  $145.6(4)$ |  – |  $151.1(6)$ |  – |  $157(2)$ |  – |  –
KLT: fixed scale |  0.143 |  – |  – |  – |  – |  – |  – |  – |  – |  – |  $157.8(7)$ |  $172(3)$ |  $197(2)$
Table 1: crossover temperatures $T_{c}$, defined as the peak position of the
reduced susceptibility, eq. 23 with $\chi_{M}^{(l,s)}$, for various pion
masses from the present computation (fRG-DSE) in comparison with results from
Braun  _et al._ (2020a) (fRG), and recent lattice QCD studies: Ding  _et al._
(2019) (hotQCD) with $N_{\tau}=8,12$ and Kotov  _et al._ (2021) (KLT) with a
fixed scaled approach. These results are also depicted in Figure 7.

###  IV.2 Combined error estimate for the chiral transition temperature

This situation calls for a combined functional-lattice study: the chiral limit
(and the continuum limit) are difficult to access within lattice computations,
while functional computations are well-controlled in the chiral limit.
However, lattice computations provide benchmark results for current quark
masses, where the continuum extrapolation can safely be done. Given the simple
dependence of $T_{c}$ on the current quark mass, as predicted by functional
approaches, a specifically important benchmark result is that of the slope of
$T_{c}(H)$ in this regime. This goes beyond the scope of the present work,
here we only present a preliminary discussion on the basis of the existing
functional and lattice data.

To begin with, we list the data of the current computation in comparison with
the fRG-results, Braun  _et al._ (2020a), and lattice results, Ding  _et al._
(2019); Kotov  _et al._ (2021) in Table 1. Given the long regime with linear
dependence of $T_{c}(m_{\pi})$ on the pion mass (or a square root dependence
on the light current quark mass) discussed in the last section, we consider

$\displaystyle D(m_{\pi})=\frac{T_{c}(m_{\pi})-T_{c0}}{T_{c0}}\,,$ (26)

the weighted difference of the transition temperatures in the chiral limit and
at the physical point. In eq. 26 the errors stems from the systematic error
estimate related to (ii) discussed in Section III.3. For the sake of a direct
comparison to the fRG and lattice results we should use the transition
temperature for $m_{\pi}=140$ MeV. In the present work, the respective
critical temperature is given by

$\displaystyle
T_{c}(m_{\pi}=140\,\textrm{MeV})=155.4^{+5.1}_{-5.5}\textrm{MeV}\,.$ (27)

see also Table 1. The error estimate in eq. 27 is derived from the systematic
error estimate (ii), leading to an error of $\pm$5 MeV for the physical point,
see eq. 24. With eq. 27 and eq. 25 we get

$\displaystyle D(m_{\pi}=140\,\textrm{MeV})=0.101^{+0.013}_{-0.018}\,,$ (28)

where the error estimate derives from that of $T_{c0}$ in eq. 25 and eq. 27.
Note also that the ratio grows with increasing coupling strength of the quark-
gluon vertex.

Figure 7: Linear extrapolation of $T_{c}$ as a function of $M_{\pi}$ with the
lattice data from Ding  _et al._ (2019) (hotQCD) and Kotov  _et al._ (2021)
(Kotov et al) and the present functional results and that of Braun  _et al._
(2020a) (fRG: Braun et al). The red area indicates extrapolations of the
$N_{\tau}=12$ results with or without the data point at $m_{\pi}=80$ MeV.

The result eq. 28, obtained in the present DSE approach has to be compared
with the respective predictions from the fRG approach in Braun  _et al._
(2020a) and lattice predictions Ding  _et al._ (2019); Kotov  _et al._ (2021)
in Table 1. While the functional results are in quantitative agreement, the
lattice results and functional results are still compatible, if one takes into
account an enlarged systematic error for the continuum and infinite volume
extrapolation on the lattice as that indicated in Table 1. In our opinion,
such an enlarged systematic error follows straightforwardly from the very
thorough studies of the finite volume extrapolation in low energy effective
theories in Braun  _et al._ (2006, 2011a); Klein (2017). Interestingly, these
works indicate a rather strong flattening of the ratio $D(m_{\pi})$, eq. 26,
in the infinite volume limit. Moreover, whether or not such a strong effect is
also present in the continuum limit is difficult to assess, leading to a
conservative extension of the systematic error in comparison to that shown in
in Table 1.

However, as discussed before, the large systematic error can be crunched down
in the combined analysis with functional approaches and lattice simulations:
The situations suggests a precise lattice determination of the slope within a
current quark mass regime, in which the continuum limit can be taken safely.
Moreover, this regime should be within the regime of validity of the square
root dependence $T_{c}(m_{l})\propto\sqrt{H}$. For these current quark masses
chiral expansions such as chiral perturbation theory hold, which is discussed
in the next section.

##  V Quasi-massless modes and chiral expansions

The results for the magnetic equation of state and the deduced chiral phase
structure presented in the last section and depicted in Figure 5 and Figure 6
await an explanation of the underlying dynamics. In particular we expect the
occurrence of quasi-massless modes as well as (partial) chiral symmetry (e.g.
$O(4)$ or $O(2)$) close to the chiral transition line for sufficiently low
current quark masses. The presence of the latter is relevant for the
reliability of chiral models used for heavy ion phenomenology, and in
particular for chiral transport, e.g. Bluhm  _et al._ (2019, 2020); Grossi
_et al._ (2021); Florio  _et al._ (2021). Trivially, massless modes and chiral
symmetry is present in the critical scaling regime close the the chiral limit.
This has spurred many studies of the size of the critical region, and we
discuss our present findings in Section V.1 in comparison to that in the
literature. In short, the present work is compatible with the findings in
previous functional studies that hint at a very small critical region.
However, we emphasise that the validity of these chiral transport models or
chiral fluid dynamics is not restricted to the critical region, they simply
require quasi-massless modes and (partial) chiral symmetry. These properties
may be present even far outside the critical region, and are tied to the
chiral dynamics.

The details of the chiral dynamics are conveniently accessed with the order
parameter potential, $V_{\textrm{eff}}(\Delta_{l})$. Higher order terms
corresponded to multi-quark scattering processes in the scalar interaction
channel, see in particular Fu  _et al._ (2020b). In Section V.2 we derive the
order parameter potential in the DSE approach for the first time, and discuss
the relation of these processes for the scaling of the chiral condensate as
well as that of the transition temperature with the light current quark mass.
In Section V.3 the potential different scaling regimes are then derived from
the quark gap equation as well as the magnetic susceptibility.

Finally, in Section V.4 we use our results on the temperature and current
quark mass dependence of the chiral condensate as well as that on the
effective order parameter potential to estimate the size of the chiral regime
with quasi-massless modes, and the sub-regimes with different chiral dynamics.

###  V.1 Critical scaling and the size of the critical region

The size of the regime with critical scaling and the respective universality
class have been discussed at length in the literature, see e.g. Braun  _et
al._ (2020a); Ding  _et al._ (2019); Kotov  _et al._ (2021) . In particular it
has been argued in Braun  _et al._ (2020a), that the scaling window for
$m_{\pi}\to 0$ in 2+1 flavour QCD is very small and critical scaling is not
visible for $m_{\pi}\gtrsim 30$ MeV. This finding in functional QCD is
supported by respective findings in low energy effective theories, see e.g.
Braun and Klein (2008, 2009); Braun  _et al._ (2011a) and the review Klein
(2017). In the latter works it is also shown that apparent critical scaling
can be found in regimes far away from the scaling window: These regimes do not
show critical scaling, but it may require a high accuracy of the numerical
data for extracting its absence.

While this intricate situation has only been demonstrated explicitly for low
energy effective theories, the respective theories carry the dynamics of the
critical modes in QCD. In terms of scaling, QCD should be seen as a driven
model where the driving terms are provided by the glue dynamics. The latter
dynamics does not trigger chiral criticality as it carries no chiral symmetry.
Typically, the scaling windows of driven models with a non-critical driving
force shrink in comparison to the model in the absence of driving forces. In
conclusion the above findings strongly suggest a very small scaling window in
QCD for pion masses $m_{\pi}\ll m_{\pi,\textrm{phys}}$.

The present DSE study cannot add much to this specific question as the current
approximation does not include the full dynamics of the chiral critical modes
in contradistinction to the fRG study in Braun  _et al._ (2020a). However, as
shown in Section IV.1 and Section IV.2, the present results agree
quantitatively with that of Braun  _et al._ (2020a), which corroborates the
irrelevance of critical scaling for the magnetic EoS, the transition
temperature $T_{c}(H)$ and in particular for $T_{c0}$. Explicitly, we consider
the finite volume expansion in the 3d $O(4)$ universality class as used e.g.
in Ding  _et al._ (2019),

$\displaystyle
T_{p}(H,L)=T_{c0}\left(1+\frac{z_{X}(z_{L})}{z_{0}}H^{\frac{1}{\beta\delta}}\right)+c_{X}H^{1-\frac{1}{\delta}+\frac{1}{\beta\delta}},$
(29)

where $L$ is the size of the system, and $L\rightarrow\infty$ in the continuum
limit. Equation 29 also holds for other universality classes, and it requires
the computation of the two critical exponents $\beta$ and $\delta$. The former
one can be derived from the temperature-dependence of the chiral condensate in
the chiral limit in the critical region with $T\to T_{c0}$ from below,

$\displaystyle\Delta_{l,\chi}(T)\propto\left(\frac{T_{c0}-T}{T_{c0}}\right)^{\beta}\,,\quad\beta=\frac{1}{2}\nu(d-2+\eta_{\Delta})\,,$
(30)

with the spacial dimension $d$. The critical exponent $\delta$ describes the
scaling with the light current quark mass in the critical region with $H\to 0$
and $T=T_{c0}$,

$\displaystyle\Delta_{l}(T_{c0},H)\propto
H^{\frac{1}{\delta}}\,,\qquad\delta=\frac{2+d-\eta_{\Delta}}{d-2+\eta_{\Delta}}\,.$
(31)

In Section V.3 we will derive (non-critical) scaling laws for the chiral
condensates in order to determine the size of the regime with massless modes
rather than that with critical scaling. In any case, for the present purpose
of critical scaling the asymptotic cases are,

  1. (1)

Asymptotically large critical region: Then all relevant pion masses are
covered by the scaling regime, and we can globally use eq. 29. In case of 3d
$O(4)$ scaling, the scaling exponents are given by $\beta=0.379$,
$\delta=4.820$ and $1/(\beta\delta)=0.547$. For comparison, the 3d $O(2)$
scaling exponents are given by $\beta=0.349$, $\delta=4.780$ and
$1/(\beta\delta)=0.599$, the 3d $Z(2)$ scaling exponents of a potential CEP
are given by $\beta=0.326$, $\delta=4.805$ and $1/(\beta\delta)=0.638$.

  2. (2)

Asymptotically small critical regime with 3d $O(4)$ scaling: This implies that
the order parameter potential has a polynomial expansion for most parameter
values. We shall see that

$\displaystyle T_{p}(H,L)=T_{c0}\left(1+a(L)H^{\frac{1}{2}}+b(L)H\right)\,,$
(32)

fits the data well even in the chiral limit $m_{l}=0$ with $m_{\pi}=0$. The
$H$-scaling of $T_{p}$ is linked to that of $\Delta_{l}$ in the presence of
quasi-massless modes, which is discussed in detail in the next sections,
Section V.2 and Section V.3.

Now we use our data to extrapolate to $T_{c0}$ assuming either case (1) with
3d $O(4)$ scaling or case (2). We found that both extrapolations lead to
$T_{c0}\approx 141$ MeV. Similar results are achieved for $O(2)$ or $Z(2)$
critical exponents. Indeed, our data are matched best with eq. 32, and the
lattice data are also well compatible with this trivial relation that follows
from a simple chiral expansion. We also note that a possible volume dependence
of $c_{X}$ and further volume dependences largely add to a systematic error
estimate, as discussed at the end of Section IV.2, see also Figure 7. This
leads us to a conservative combined estimate for $T_{c0}$,

$\displaystyle 132\,\textrm{MeV}\lesssim T_{c0}\lesssim 141\,\textrm{MeV}\,,$
(33)

where the range is solely due to the combined systematic error estimate based
on functional results from the present DSE approach and from the fRG
computation in Braun  _et al._ (2020a), as well as lattice results from Ding
_et al._ (2019); Kotov  _et al._ (2021). While the two approaches also allow
for smaller (lattice QCD) or larger (functional QCD) temperatures within the
respective error estimates, in combination these temperatures are disfavoured.

In summary this situation calls for a combined study with functional and
lattice approaches for crunching down the relatively large systematic error in
eq. 33, exploiting the respective strengths. This task is left to future work.

###  V.2 Order parameter potential

All scalings, and in particular the critical scaling discussed in the last
Section, can be conveniently extracted from the temperature-dependent order
parameter potential $V_{\textrm{eff}}(\Delta_{l},m_{l})$ at fixed $m_{s}$. It
also offers a more direct access to the related physics: The effective
potential of the order parameter $\Delta_{l}$ can be derived from the
effective action $\Gamma[A,c,\bar{c},q,\bar{q};\Delta_{l},m_{l}]$ at fixed
$m_{s}$, where we have introduced a current for $\Delta_{l}$ and applied a
respective Legendre transform. In the fRG approach this is done with dynamical
hadronisation, more details in the present context can be found in Fu  _et
al._ (2020b); Braun  _et al._ (2020a); Dupuis  _et al._ (2021). Then, the
order parameter potential is given by a convex function
$V_{\textrm{eff}}(\Delta_{l},m_{l})=\Gamma[0;\Delta_{l},m_{l}]/$(volume).
Moreover, it can be shown that the only term with the light current quark mass
is linear in $m_{l}$. For a recent detailed discussion of this fact see Fu
_et al._ (2020b). We have

$\displaystyle
V_{\textrm{eff}}(\Delta_{l},m_{l})=V_{\chi}(\Delta_{l})-c_{\Delta}\,m_{l}\,\Delta_{l}\,,$
(34)

with a dimensionless constant $c_{\Delta}$. While such a split very much
resembles chiral expansions, it is important to note that eq. 34 is valid for
all masses. The temperature and current-quark mass dependent expectation value
of the condensate is determined by the solution of the equation of motion
(EoM) $\Delta_{\textrm{EoM}}$,

$\displaystyle\left.\frac{\partial
V_{\chi}(\Delta_{l})}{\partial\Delta_{l}}\right|_{\Delta_{l}=\Delta_{l,\textrm{EoM}}(T,m_{l})}=c_{\Delta}\,m_{l}\,.$
(35)

Hence, eq. 34 carries the thermal and current mass dependence of the chiral
condensate, and can be used to extract the transition temperature. The
potential $\bar{V}_{\chi}$ gives us easy access to the scaling behaviour of
the pion mass. On the EoM for $\Delta_{l}$ we have,

$\displaystyle m_{\pi}^{2}\propto\frac{\partial
V_{\chi}(\Delta_{l})}{\partial\Delta_{l}^{2}}=\frac{c_{\Delta}m_{l}}{2\Delta_{l,\textrm{EoM}}(T,m_{l})}\,,$
(36)

which also allows us to recover $\chi$PT in the vacuum, for a detailed
discussion see e.g. Fu  _et al._ (2020b). The effective potential also carries
the (quark) pressure with
$V_{\chi}(\Delta_{\textrm{EoM}})-c_{\Delta}m_{l}\,\Delta_{\textrm{EoM}}$,
which has been used in Isserstedt  _et al._ (2021) for a determination of the
EoS. The temperature-independent normalisation $c_{\Delta}$ can be
conveniently fixed at large temperatures, where the pressure resumes the
Stefan-Boltzmann limit. However, in the present work we are not interested in
the latter, and for the sake of simplicity of the following considerations we
normalise eq. 34 with $c_{\Delta}$,

$\displaystyle\bar{V}_{\textrm{eff}}(\Delta_{l},m_{l})=\bar{V}_{\chi}(\Delta_{l})-m_{l}\,\Delta_{l}\,,$
(37)

with $\bar{V}_{\chi}(\Delta_{l})=V_{\chi}(\Delta_{l})/c_{\Delta}$ is the
potential in the chiral limit, divided by the dimensionless constant
$c_{\Delta}$. Formulated in terms of $\bar{V}_{\chi}$, the EoM in eq. 35 turns
into

$\displaystyle\frac{\partial\bar{V}_{\chi}(\Delta_{l})}{\partial\Delta_{l}}=m_{l}(\Delta_{l})\,.$
(38)

Equation 38 has already been adapted to our present results: we have computed
$\Delta_{l}(T,m_{l})$ for a given $H=m_{l}$, which provides us with pairs of
$(\Delta_{l}(T),m_{l}(T))$ for all temperatures. Accordingly, the chiral
effective potential $V_{\chi}(\Delta_{l},T)$ for
$\Delta_{l}(T)>\Delta_{l,\chi}(T)$ can be represented as

Figure 8: Chiral effective potential $V_{\chi}(\Delta_{l})$ as a function of
the chiral condensate $\Delta_{l}$ and temperature $T$. The boundary of the
flat regime is the solution of the EoM in the chiral limit.
$V_{\chi}(\Delta_{l})$ is fit well with the polynomial potential eq. 40
(broken regime) and eq. 42 (symmetric regime), the respective couplings are
depicted in Figure 9.

$\displaystyle\bar{V}_{\chi}(\Delta_{l},T)=\int_{\Delta_{\chi}(T)}^{\Delta_{l}}d\Delta\,m_{l}(\Delta,T)\,.$
(39a) For $\Delta_{l}<\Delta_{i,\chi}$ we have
$\bar{V}_{\textrm{eff}}(\Delta_{l},T)=\bar{V}_{\textrm{eff}}(\Delta_{l,\chi},T)$
due to convexity of $V_{\chi}$. The $\Delta$-integral in eq. 39a can be turned
into an $H$-integral with the help of the chiral susceptibility eq. 22,
$\displaystyle\bar{V}_{\chi}(\Delta_{l},T)=-\frac{1}{2}\int_{0}^{m_{l}^{2}(\Delta_{l})}dm^{2}_{l}\,\chi^{(l)}(m_{l},T)\,,$
(39b)

for $\Delta_{l}>0$. In the presence of a first order phase transition, eq. 39
has to be slightly modified, but a similar expression still holds true. We
also note that the potential is even in $\Delta_{l}\to-\Delta_{l}$ which
reflects the symmetry under $m_{l}\to-m_{l}$.

This leaves us with the following task: For a given temperature $T$ we compute
numerically the $m_{l}$\- or $H$-dependence of the light chiral condensate
$\Delta_{l}$, leading to a fit of $\Delta_{l}(T,H)$. The integral
representation eq. 39 of the chiral effective potential requires
$m_{l}(T,\Delta_{l})$ or $H(T,\Delta_{l})$, which is extracted by numerically
inverting this relation. The result is depicted in Figure 8. The flat region
for $\Delta_{l}\leq\Delta_{l,\chi}(T)$ shrinks with temperature and vanishes
for $T\geq T_{c0}$. Moreover, we clearly see the flattening of the potential
in the vicinity of $T_{c}$ for $\Delta_{l}>\Delta_{l,\chi}(T)$, indicating the
emergence of quasi-massless modes.

Figure 9: Temperature dependence of the mass and coupling parameters in the
polynomial potentials eq. 40 with eq. 41 (broken regime) and eq. 42 (symmetric
regime). The couplings have been normalised by their values at $T=0$. At
vanishing temperature, $T=0$, we have $m_{\Delta}^{2}=0.24$ GeV2,
$\lambda_{4}=3.3\times 10^{3}$, $\lambda_{6}=3\times 10^{6}$ GeV-2.

We also can provide simple polynomial fits for the effective potential shown
in Figure 8: In the broken regime we consider a polynomial potential,

$\displaystyle
V_{\chi}(\Delta_{l})=\frac{\lambda_{4}}{8}(\Delta_{l}^{2}-\Delta_{l,\chi}^{2})^{2}+\frac{\lambda_{6}}{24}(\Delta_{l}^{2}-\Delta_{l,\chi}^{2})^{3}\,,$
(40)

where light chiral condensate in the chiral limit, $\Delta_{l,\chi}\geq 0$,
minimises the potential. Equation 40 is the potential of the scalar
($\sigma$-) mode, and the respective ’mass’ in the broken regime is

$\displaystyle
m_{\Delta_{l}}^{2}=\partial_{\Delta_{l}}^{2}V_{\chi}(\Delta_{l}=\Delta_{l,\chi})=\lambda_{4}\,\Delta_{l,\chi}^{2}\,.$
(41)

Hence, eq. 40 provides a three-parameter fit to the potential in the broken
regime in an expansion about the minimum $\Delta_{l,\chi}>0$. The parameters
are $\Delta_{l,\chi},\lambda_{4},\Lambda_{6}$, the mass at the minimum is a
derived quantity. However, with eq. 41 we can also use $m_{\Delta_{l}}^{2}$ as
parameter instead of $\Delta_{l,\chi}^{2}=m_{\Delta_{l}}^{2}/\lambda_{4}$,
leaving us with the three parameters
$m_{\Delta_{l}}^{2},\lambda_{4},\Lambda_{6}$.

Figure 10: Order parameter potential $V_{\chi}(\Delta_{l})$ (red squares) for
the chiral phase transition temperature in the chiral limit, $T_{c0}$. The
potential is shown for the regime eq. 44 with $\Delta_{\textrm{ch}}=0.026$
GeV3. We also display the fits computed from the regime eq. 43 with
$\Delta_{\textrm{max}}=0.017$ GeV3 (_dashed curve_), as well as that computed
from the regime $\Delta_{l}\in[\Delta_{\textrm{max}},\Delta_{\textrm{ch}}]$
(_solid curve_). Outside the regime eq. 44 the latter fit fails.

In the symmetric regime we have the fixed expansion point $\Delta_{l,\chi}=0$,
and the expansion point parameter is in any case traded for the explicit mass
parameter. Again this leaves us with the three parameters,
$m_{\Delta_{l}}^{2},\lambda_{4},\Lambda_{6}$, and the the potential is given
by,

$\displaystyle
V_{\chi}(\Delta_{l})=\frac{m_{\Delta_{l}}^{2}}{2}\Delta_{l}^{2}+\frac{\lambda_{4}}{8}\Delta_{l}^{4}+\frac{\lambda_{6}}{24}\Delta_{l}^{6}\,.$
(42)

Both, in the broken and in the symmetric regime we have a Gaußian mass term,
and and a standard $\Delta_{l}^{4}$ interaction term with the interaction
strength $\lambda_{4}$. The $\Delta_{l}^{6}$-term carries higher order
interactions of the scalar interaction channel of quark–anti-quark
scatterings. The parameters are now obtained in a fit regime about the
respective minimum $\Delta_{l,\chi}$, which is one of the parameters in the
broken regime. We use a $\chi^{2}$-fit in a fit regime

$\displaystyle\Delta_{l}\in[\Delta_{l,\chi}(T),\Delta_{\textrm{max}}(T)]\,.$
(43)

In eq. 43, $\Delta_{\textrm{max}}$ is the maximal $\Delta_{l}$, for which the
fit parameters as well as the $\chi^{2}$ is stable. In the vicinity of
$T_{c0}$ the respective $\Delta_{\textrm{max}}$ is shrinking, and we also use
additional polynomial fits in the regime
$\in[\Delta_{\textrm{max}},\Delta_{\textrm{ch}}]$ for an estimate of the
regime dominated by the $\Delta_{l}^{4}$ term. It is suggestive to speculate
that chiral expansions work in the combination of these regimes, that is the
regime

$\displaystyle{}\Delta_{l}\in[\Delta_{l,\chi}(T),\Delta_{\textrm{ch}}(T)]\,.$
(44)

Note however, that the regime eq. 44 does not sustain a global fit as close to
$T_{c0}$ the potential is dominated by a higher order polynomial.

The temperature dependence of the three parameters, normalised by their value
at vanishing temperature, is depicted in Figure 9. The steep drop of both,
$m_{\Delta_{l}}^{2}(T)/m_{\Delta_{l}}^{2}(0)$ and
$\lambda_{4}(T)/\lambda_{4}(0)$ for $T\to T_{c0}$ from below signals the
approach to the chiral phase transition. The drop of $\lambda_{4}$ entails
that the $\lambda_{6}$-term dominates, leading to an $H^{1/5}$ scaling, as
discussed below. Note however, that this does not necessarily entail critical
scaling in this regime. In Figure 10 we show the full potential at the chiral
transition temperature in the chiral limit, $T_{c0}$, (red squares) together
with the monomial fit with $\Delta_{l}^{6}$ in the regime eq. 43 with
$\Delta_{\textrm{max}}(T_{c0})=0.017$ GeV3 as well as a polynomial fit in the
regime $\Delta_{l}\in[\Delta_{\textrm{max}},\Delta_{\textrm{ch}}]$ with
$\Delta_{\textrm{ch}}=0.026$ GeV3. For chiral condensates larger than
$\Delta_{\textrm{ch}}$ we approach the regime of asymptotically large current
quark masses.

We proceed with discussing different regimes in the $T,H$-plane with different
chiral or non-chiral dynamics. This dynamics is carried by the interaction
terms in the effective potential $V_{\chi}(\Delta_{l})$. We note that for
(fractional) monomial scaling of $\Delta_{l}$ with $H$ the respective term in
the effective action can be deduced analytically. This gives us an analytic
understanding of the $H$-scaling of $\Delta_{L}$ and the underlying QCD
dynamics. For example, monomial potentials lead to specific monomial
$H$-scalings,

$\displaystyle
V_{n}(\Delta_{l})\approx\frac{\lambda_{n}}{n!}(\Delta_{l}^{2}-\Delta^{2}_{l,\chi})^{\frac{n}{2}}\,,\quad\to\quad\Delta_{l}(H)\propto
H^{\frac{1}{n-1}}\,,$ (45)

for large $H$. Equation 45 provides the link of the $H$-dependence of the
condensate for different temperature with the dynamics of the chiral
condensate as $\Delta_{l}^{n}$ terms in the effective potential simply entail
the multi-scattering of scalar quark-antiquark channels. This relation is
directly seen within dynamical hadronisation, see e.g. Fu  _et al._ (2020b);
Braun  _et al._ (2020a).

The inversion property $\Delta^{n}\to H^{1/(n-1)}$ displayed in eq. 45 entail
that higher order terms in the potential lead to terms that dominate an
expansion of the chiral condensate for small $H$. Consequently one can hope
for accurate global fits $\Delta_{l}(H)$ in the potential validity regime eq.
44 of chiral expansions. This property will be used in section V.4.

Moreover, we emphasise that specifically for $n=4,6$ the scaling in eq. 45 is
the 3d mean field scaling, $\delta=3$, and 3d dimensional scaling, $\delta=5$,
respectively. The latter is derived from eq. 31 with a vanishing anomalous
dimension, $\eta_{\Delta}=0$. Given the small anomalous dimensions for 3d
$O(4)$, 3d $O(2)$ (and 3d $Z_{2}$) universality classes, the dimensional
scaling is very close to the full universal scalings, and in particular to the
$O(4)$-scaling. In the present DSE approach a non-vanishing anomalous
dimension $\eta_{\Delta}$ requires a full resolution of the (scalar-
pseudoscalar) four-quark interaction, typically included via pion
contributions. We have not taken into account (the thermal part) of the
resonant scalar-pseudoscalar four-quark interaction leading to
$\eta_{\Delta}=0$. Note however, that this has been done in Braun  _et al._
(2020a) for $m_{\pi}\geq 30$ MeV with no sign of criticality. Still, as
mentioned before, the present work does not add much to the intricate question
of the size of the critical regime. Note also that in the vacuum we expect an
$H^{1/3}$-scaling for large, but not asymptotically large H, which we call
trivial scaling as it is far away from criticality (at asymptotically large
$H$ we expect an $H^{3}$ scaling fixed by dimensionality). We emphasise that
this makes it extremely difficult to disentangle trivial scaling regimes from
critical scaling regimes. This entails, that both, the proof of presence as
well as that of the absence of critical scaling requires data with exceedingly
small statistical error. This has been studied in detail in Braun and Klein
(2008, 2009); Klein (2017), advocating both, a small critical regime and the
necessity for exceedingly small statistical error.

Indeed, the worst case in this respect is a trivial scaling regime close to
$T_{c0}$, that exhibits an $H^{1/5}$-scaling, stemming from a
$\Delta_{l}^{6}$-potential. As the quadratic term vanishes at $T_{c0}$,
$\lambda_{2}(T\to T_{c0})\to 0$, the presence of this scenario is linked to
the fate of $\lambda_{4}(T\to T_{c0})$. We shall see that also
$\lambda_{4}(T\to T_{c0})\to 0$ and we may be left with a trivial
$H^{1/5}$-scaling that turns into a critical one with
$H^{1/(\delta\beta)}$-scaling in the critical regime. How to distinguish this
scenario from the one with a large regime with critical scaling, a more
refined analysis is required. This analysis is deferred to a forthcoming
publication, fQCD collaboration (2021).

While this intricacy makes it difficult to disentangle critical from trivial
scaling, it is very good news for chiral expansions and extrapolations. In
summary, the present analysis gives full access to the phenomenologically
relevant validity regime of chiral expansions with quasi-massless modes. It
also corroborates the point of view, that for phenomenological applications
and experimental signatures in heavy ion collisions the search for critical
scaling and its signatures is rather a red herring than a property to look
for.

###  V.3 Quasi-massless modes

With these preparations we now can tackle the phenomenologically important
question after the regime with quasi-massless modes close to the transition
line, already mentioned in the introduction of Section V. Chiral transport
models or chiral fluid dynamics both rely on the existence of quasi-massless
modes such as the pion or the sigma mode close to the potential critical end
point in the QCD phase structure. Certainly critical modes are quasi-massless
in the critical regime, but this property is more generally connected to the
validity regime of chiral expansions.

The most prominent chiral expansion is chiral perturbation theory ($\chi$PT),
an expansion about QCD in the chiral limit in the vacuum. It is a low energy
effective theory of QCD, and the expansion coefficients or low energy
constants are either determined from first principles QCD computation in the
chiral limit or they are fixed by experimental data. We have already used
$\chi$PT results for the determination of the pion mass in the present
approach, see eq. 19a and below. $\chi$PT is impressively successful in the
vacuum as well as small enough temperatures and densities. In $\chi$PT, the
chiral condensate has a chiral expansion similar to that in eq. 19a for the
pion mass and pion decay constant. We parametrise,

$\displaystyle\Delta_{l}(T,H)=$
$\displaystyle\,\Delta_{l,\chi}(T)\Bigr{[}1+\left(c_{1}+c_{\textrm{log}}\,\log
H\right)\,H+O(H^{2})\Bigr{]}\,,$ (46a) with temperature-dependent coefficients
$c_{\textrm{log}}(T)$, $c_{1}(T)$. The estimate eq. 46a can only hold for
temperatures below $T_{c0}$, so at most $T<T_{c0}$.

Note also that eq. 46a can be reformulated as a chiral expansion such as eq.
19a with the expansion parameter eq. 19b. The logarithm in eq. 46a originates
in massless pion loops in the vacuum within the expansion about $H=0$ with
$m_{\pi}=0$. There are higher order logarithms well-known in chiral
perturbation theory, which we discard for the present estimate. Moreover, at
sufficiently large temperature we do not expect a logarithmic scaling any
more, as the infrared singularities at finite temperature are 3d ones and have
a rational scaling. Indeed, we shall see that eq. 46a does not work for
$T\gtrsim 50$ MeV, see Figure 13.

Finally, for asymptotically large quark masses the current quark mass
dominates the constituent quark mass and hence the condensate. There we expect
a trivial linear $m_{l}$-scaling of the renormalised condensate. Whether this
transition already takes place for $H\leq 1$ is a dynamical question, which
also depends on the size of the strange quark mass.

In summary, we infer from this discussion, that an $H$-expansion of the chiral
condensate for small $H$ and small temperature at least has to incorporate the
logarithm and polynomial terms known from $\chi$PT. In turn, for larger $H$
and temperatures the logarithm is absent.

Hence, for sufficiently large temperatures we use a rational expansion without
logarithm, as the latter is absent in a three-dimensional chiral expansion,

$\displaystyle\Delta^{(\chi\textrm{EP})}_{l}(T,H)\approx$
$\displaystyle\,\Delta_{l,\chi}(0)\Bigl{(}c_{0}+c_{r_{1}}H^{r_{1}}+c_{r_{2}}H^{r_{2}}+c_{1}H\Bigr{)}\,,$
(46b)

for the extension of the $\chi$PT regime to a chiral expansion regime
($\chi$EP) at larger temperatures and small current quark masses. Here,
$H^{r_{1}}$ is the leading scaling for $H\to 0$, and hence $r_{1}<r_{2}$. As
in the vacuum, we expect a linear $H$-scaling of the chiral condensate for
sufficiently large current quark masses. Then, the explicit chiral symmetry
breaking via $H$ dominates over the spontaneous one. This linear regime is
expected to set-in at smaller current quark masses for larger temperature, as
the chiral condensate in the chiral limit, which measures the amount of
spontaneous chiral symmetry breaking, reduces with the temperature. Hence, eq.
46b may also allow for good global fits, in contradistinction to $\chi$PT,
where the logarithmic small $H$-term would also dominate the large $H$-regime.

We proceed by discussing the different $H$-scalings we expect to see from
general diagrammatic arguments as well as within the present approximation.

To begin with, for constant vertices the gap equation allow for a simple
extraction of mean field scaling for $\Delta_{l}\to 0$. We shall use that for
the purpose of a leading order scaling analysis that the wave functions
$Z_{q}(p)$ and $Z_{q}^{\parallel}(p)$ of the quark propagator in eq. 6 only
show a mild momentum dependence, see e.g. the recent DSE review Fischer
(2019), and Gao and Pawlowski (2021, 2020); Gao  _et al._ (2021a), where the
DSE framework used in the present work is put forward. Hence, for the present
scaling analysis the quark propagator is well approximated by

$\displaystyle\frac{M_{q}-\text{i}\,/{p}}{p^{2}+M_{q}^{2}}\,,$ (47)

where $M_{q}=M_{q}(p=0)$ is the infrared value of the constituent quark mass.
Within this approximation the quark gap equation is an integral equation for
$M_{q}$, and its resolution with a gapped gluon propagator readily leads to

$\displaystyle M_{q}(H)\propto H^{\frac{1}{3}}+O(H)\,,$ (48)

in the limit $H\to 0$. With $\Delta_{l}\propto M_{l}$ this entails the mean
field scaling for the chiral condensate. It is well-known that the behaviour
eq. 48 is already seen with the crude approximation
$G_{AA}(p)\approx\delta(p^{2}-m_{\textrm{gap}}^{2})$ with $H\ll
m_{\textrm{gap}}$. Note however, that this scaling follows as trivial scaling
from a polynomial order parameter potential $V_{\chi}$ with a (Gaußian) mass
term proportional to $\Delta_{l}^{2}$ and a $\Delta_{l}^{4}$-term (Ginzburg-
Landau potential).

Figure 11: Leading diagram in the chiral limit, $m_{\pi}\to 0$ in the
diagrammatic representation of the magnetic susceptibility. The left diagram
is the (resummed) fish diagram with two dressed four quark vertices, the right
diagram shows the dominant scalar–pseudo-scalar $(\sigma,\vec{\pi})$-channel
in the chiral limit. The pions and the $\sigma$ mode are depicted by a quark
double line, indicating its $l\bar{l}$ content.

For the full scaling the DSE requires a non-trivial quark-gluon vertex with a
full back coupling of fluctuations, including the critical ones we expect the
full $H$-scaling to be present in the gap equation. The scaling analysis of
this case is conveniently done with the chiral susceptibility: As a current
quark mass derivative of the chiral condensate it is an integrated four-quark
correlation function. It has a diagrammatic representation in terms of its DSE
of fRG equation, which already diagrammatically includes part of the back
reaction required for the full scaling. Its chiral limit, $H\to 0$, is
governed by the most divergent diagram for $H\to 0$. In the present DSE
approach this analysis is best done after introducing effective degrees of
freedom for the resonant scalar-pseudo-scalar interaction four-quark channels,
as done in the fRG approach with dynamical hadronisation, in the present
context see Fu  _et al._ (2020b); Braun  _et al._ (2020a). Alternatively one
may employ the respective resummations. In the first case this leads us to a
quark-pion box diagram, in the latter case to a fish diagram with two dressed
four-quark-vertices, see Figure 11.

For the sake of a scaling analysis the pseudoscalar resonant four-quark
interaction channels (or pion propagators) may be approximated by
$1/(p^{2}+m_{\pi}^{2})$ away from the critical scaling regime, and the light
and strange quark propagators may still be approximated by eq. 47.

Naturally, the analysis of both diagrams leads to the same chiral scaling, and
in the vacuum it leads to the well-known results in chiral perturbation
theory. At sufficiently large temperature the scaling analysis leads to

$\displaystyle\frac{\partial\Delta_{l}}{\partial
H}\propto\frac{M_{l}}{m_{\pi}(m_{\pi}+M_{l})^{3}}\,.$ (49)

Equation 49 originates from the 3d spatial momentum integral of the zeroth
Matsubara mode.

In a regime with $\Delta_{l}(H\to 0)\to 0$, we can substitute the derivative
w.r.t. $H$ in eq. 49 by $1/H$. From the integral representation of the
(renormalised) chiral condensate, see e.g. Gao and Pawlowski (2021, 2020); Gao
_et al._ (2021a) we can safely assume that $M_{q}\propto\Delta_{l}$, which
also follows straightforwardly within dynamical hadronisation, see Fu  _et
al._ (2020b). We also use the relation between chiral condensate and pion
mass, eq. 36. Note that the latter already implies in the limit $H\to 0$, that
$\Delta_{l}\propto H^{r_{1}}$ with $r_{1}<1$. This follows from
$m^{2}_{\pi}\propto H/\Delta_{l}\to 0$.

With the relations discussed above, eq. 49 turns into,

$\displaystyle\Delta_{l}\propto
H^{\frac{1}{2}}\,\frac{\Delta_{l}^{3}}{\left(c\,H^{1/2}+\Delta_{l}^{3/2}\right)^{3}}\,,$
(50)

for $\Delta_{l}\to 0$. The constant $c$ in eq. 50 is irrelevant for the
following discussion of the leading scaling behaviour.

We now proceed with a detailed analysis of the consequences of eq. 49 for the
$H$-scaling of the chiral condensate for small $H\to 0$:

$M_{q,\chi}>0$: This condition is tantamount to large spontaneous symmetry
breaking, and hence also a large chiral condensate. This situation holds true
for temperatures smaller and sufficiently far from the chiral transition
temperature. There, the constituent quark mass $M_{q}$ as well as the chiral
condensate $\Delta_{l}$ have a sizeable non-vanishing part for $H\to 0$. Then,
eq. 49 reduces to

$\displaystyle\frac{\partial\Delta_{l}}{\partial
H}\propto\frac{1}{m_{\pi}}\,.$ (51)

Moreover, the pion mass is proportional to $H^{1/2}$ as in the vacuum, see eq.
19 and eq. 36. This leads us to

$\displaystyle\Delta_{l}\propto c_{0}+c_{\frac{1}{2}}H^{\frac{1}{2}}\,,$ (52)

that is $r_{1}=1/2$. This scaling is the analogue of the $H\log H$ scaling in
the vacuum, the difference coming from the dimensional reduction 4d$\to$3d
present in the Matsubara zero mode. Clearly, in this regime quasi-massless
modes are present, as it solely originates in the small pion mass in
comparison to other mass scales. We close the discussion of this regime with
the remark, that in this limit we also expect a chiral scaling of the vertices
in the diagrams for the chiral susceptibility with powers or inverse powers of
$m_{\pi}^{2}$. This originates in the infrared singularities similar to those
leading to the $1/m_{\pi}$ scaling in eq. 51. This may change the scaling eq.
52. Indeed, the latter naively relates to a $\Delta_{l}^{3}$ term in
$V_{\textrm{eff}}$, which is at odds with the symmetry of the effective
potential under $\Delta_{l}\to-\Delta_{l}$.

$M_{q,\chi}\approx 0$: Sufficiently close to the chiral phase transition, both
the constituent quark mass $M_{q}$ as well as the chiral condensate melt away.
For $r_{1}\leq 1/3$, which potentially includes the mean field scaling, we
have $c\,H^{\frac{1}{2}}+\Delta_{l}^{2}\propto\Delta_{l}^{2}$ in the limit
$H\to 0$. Then, eq. 50 leads to

$\displaystyle\Delta_{l}\propto H^{\frac{1}{5}}\,.$ (53)

Equation 53 comprises the consistent dimensional scaling for small spontaneous
chiral symmetry breaking scales at finite temperature. In terms of the order
parameter potential it relates to a monomial $\Delta_{l}^{6}$ potential, see
eq. 45, and hence higher order scatterings of the scalar-pseudoscalar channel.
Note that this scaling is also the critical one in approximations with
vanishing anomalous dimension (of the pion), $\eta=0$, see eq. 31.

In summary, the scalings eq. 48, eq. 52 and eq. 53 indicate the presence of
quasi-massless modes. Finally, in the critical region, $r_{1}=1/\delta$ and
$r_{2}$ is the subleading scaling coefficient The critical regime has already
been discussed in detail in Section V.1 and we refer the reader to this
Section. Naturally, quasi-massless modes are present in the critical regime.

This concludes our analysis of the potential scalings of the chiral condensate
with the light current quark mass: In regimes with the scalings eq. 48, eq.
52, eq. 53 and the critical scaling eq. 31 quasi-massless modes are present,
admitting chiral expansion schemes. This supports the construction of
phenomenological chiral models in this regime.

###  V.4 Validity regime of chiral expansions

In summary we are led to a global fit for the chiral condensate
$\Delta_{l}(T,H)$ in the validity range of chiral expansions, excluding the
$\chi$PT region. In this regime Equation 46b with $r_{1}=1/5$ and $r_{2}=1/3$
is valid, as already argued in Section V.3. This reads,

$\displaystyle\Delta_{l}(T,H)\approx$
$\displaystyle\,\Delta_{l,\chi}(0)\Bigl{(}c_{0}+c_{\frac{1}{5}}H^{\frac{1}{5}}+c_{\frac{1}{3}}H^{\frac{1}{3}}+c_{1}H\Bigr{)}\,,$
(54)

and the temperature-dependent coefficients $c_{i}(T)$ with $i=0,1/3,1/5,1$ are
depicted in Figure 12. The different terms follow from a polynomial potential:
the linear $H$ term from $\Delta_{l}^{2}$, the $H^{1/3}$ term from
$\Delta_{l}^{4}$, and the $H^{1/5}$ term from $\Delta_{l}^{6}$, in combination
this leads us to the potential eq. 40 already discussed before.

The highest order term $\Delta_{l}^{6}$ in this potential relates to a
stronger chiral dynamics including multi-scatterings. This terms dominates for
small spontaneous symmetry breaking close to the chiral transition temperature
$T_{c0}$, and can be inferred from the scaling analysis of the susceptibility,
see also the discussion about Equation 53. In conclusion, as already discussed
below eq. 45, we expect the higher order terms in the potential to be relevant
for small $\Delta_{l}$ and hence small $H$. Therefore, global fits to
$\Delta_{l}(T,H)$ may work well.

The above analytic analysis is confirmed within a numerical leading order
analysis for smaller quark masses and close to the chiral transition
temperature. The scaling coefficient $r_{1}$ as well as $\beta$ have been
determined from our numerical data as

$\displaystyle\beta={0.400}({9})\,,\qquad r_{1}={0.200}({6})\,,$ (55)

with $\beta$ being extracted from $c_{0}(T)$, as
$\Delta_{l,\chi}(T)=\Delta_{l,\chi}(0)\,c_{0}$ for $T\leq T_{c0}$. In this
regime the coefficient $c_{1}(T)$ of the linear $H$ term decreases towards the
chiral transition temperature $T_{c0}$, where it vanishes, see Figure 12. We
emphasise that the scaling eq. 55 in this regime with $r=0.2$ is close to, but
distinct from the critical scaling, $r_{O(4)}\approx 0.207$ for $O(4)$-scaling
and $r_{O(2)}\approx 0.209$ for $O(2)$-scaling. This can be inferred from the
respective $\chi^{2}$. We rush to add, that within the present approximation
we do not expect full critical scaling in the critical regime, as discussed
before. However, in the corresponding fRG analysis in Braun  _et al._ (2020a)
even a polynomial scaling with $r=1/5$ can be distinguished from the critical
scaling with either $r_{O(4)}$ or $r_{O(2)}$. For the meson masses considered
there ($m_{\pi}\gtrsim 30$ MeV, that is $m_{l}\gtrsim 0.1$ MeV), no critical
scaling but polynomial scaling has been found.

Figure 12: Temperature dependence of the coefficients
$c_{0},c_{1},c_{\frac{1}{3}},c_{\frac{1}{5}}$ in the global fit eq. 54 for
$\Delta_{l}(T,H)$.

The coefficients in Equation 54 are now obtained from $H$-fits in regimes
$0\leq H\leq H_{\textrm{max}}$ for $\Delta_{l}(T,H)$ for a given temperature
$T$. Here, $H_{\textrm{max}}$ is the maximal $H$ for which the coefficients as
well as the $\chi^{2}$ error is stable. In turn, the complementary regime (for
$H\leq 1$) with $H_{\textrm{max}}<H\leq 1$ has a linear $H$ scaling and does
not admit a chiral expansion.

Figure 13: Validity range of chiral expansions in the current quark mass and
temperature plane: yellow (vacuum $\chi$PT ), blue (dominating
$H^{1/3}$-scaling), green (dominating $H^{1/5}$-scaling). For the definition
of dominating scaling see Section V.4. The grey area indicates the linear $H$
scaling for large $m_{l}$ (no chiral expansion). The green bullet indicates
the intersection point $(T\,,\,H)\approx(152\,{\rm MeV}\,,\,0.023)$ with
$m_{\pi}\approx 110$ MeV of the green regime with the chiral transition line.
The blue bullet indicates the intersection point $(T\,,\,H)\approx(166\,{\rm
MeV}\,,\,0.14)$ with $m_{\pi}\approx 270$ MeV of the blue regime with the
chiral transition line. We have also include an estimate of the small critical
scaling regime (red area, discussed above eq. 29). The vertical line indicates
QCD with physical current quark masses, that is $H=m_{l}/m_{s}=1/27$.

The fit Equation 54 with the temperature dependent coefficients depicted in
Figure 12 allows us to provide a more detailed depiction of all the different
regimes. For the present qualitative analysis we define sub-regimes in the
validity regime of chiral expansions with the relative dominance of the
respective $c_{r_{i}}H^{r_{i}}$ term over the other term:
$c_{r_{i}}H^{r_{i}}\geq c_{r_{j}}H^{r_{j}}$ for all $j=1/3,1/5,1$: The
respective regimes are indicated as (green ($H^{1/5}$), blue) ($H^{1/3}$),
yellow (vacuum $\chi$PT) in Figure 13. The regime with large light current
quark masses is indicated with grey.

Specifically, we find a large regime with dominant $H^{1/3}$-scaling. As
discussed before, this behaviour originates from an order parameter potential
with a $\Delta_{l}^{4}$ interaction term. A subleading Gaußian
$\Delta_{l}^{2}$-term leads to an additional linear term in $H$. The absence
of higher order terms with $\Delta_{l}^{2n}$ with $n\geq 3$ in the full
effective potential $V_{\textrm{eff}}$ indicates a weak chiral dynamics.
Hence, while such a scaling resembles (critical) mean field scaling, we rather
interpret this as trivial chiral scaling, related to a weak chiral dynamics.
Evidently, the latter works in favour of the convergence of the expansion.

Most importantly, in our opinion, the presence or absence of critical scaling
in a small critical regime is not relevant for the phenomenological
application of chiral transport models or chiral fluid dynamics, however, the
presence of quasi-massless modes is: clearly, in the regime with
$H^{1/5}$-scaling (green in Figure 13) quasi-massless modes are present
($\lambda_{2}\approx 0$) and dominate the dynamics. A detailed analysis of the
respective order parameter potential at $T_{c0}$ has already been provided in
Section V.2, see in particular Figure 10 and the related discussion. In
particular, in this regime the lowest order scattering of the scalar
interaction channel is absent and the dynamics is only carried by the
$\Delta^{6}$-term. This suggests to study the impact of this dynamics on QCD
transport.

##  VI Summary

In this work we have studied the magnetic equation of $2$\- and $2+1$-flavour
QCD, mostly concentrating on the latter physical case. This has been done in a
generalised functional first principles approach to QCD, set-up in Gao and
Pawlowski (2020, 2021); Gao  _et al._ (2021a). There it has been used for the
phase structure of QCD as well as precision computations in the vacuum.
Importantly, in this approach no phenomenological parameters has to be
included.

We have computed the quark condensate and chiral susceptibility for general
light current quark masses, while keeping the strange quark mass fixed in
$2+1$-flavour QCD. The magnetic equation of state and the light current quark
mass dependence of the chiral transition temperature is in quantitative
agreement with that from the recent functional renormalisation group study
Braun  _et al._ (2020a), see in particular Figure 4, Figure 5 and Figure 6.
Moreover, the critical temperature in the chiral limit is given by
$T_{c0}=141$ MeV. We also discussed a combined systematic error estimate from
the present functional results and the present work, as well as lattice
results from Ding  _et al._ (2019); Kotov  _et al._ (2021), leading to a range
$132\,\textrm{MeV}\lesssim T_{c0}\lesssim 141\,\textrm{MeV}$, see eq. 33. We
stress that this estimate could be much reduced in a combined study with
functional and lattice approaches, exploiting the respective strengths.

So far, functional computations suggest a very small scaling window with
critical scaling from low energy effective theories, Braun and Klein (2008,
2009); Braun  _et al._ (2011a); Klein (2017), as well as functional QCD, Braun
_et al._ (2020a), see Section V. The present approximation does not
incorporate the full back coupling of the chiral dynamics, and hence cannot
add to this intricate question. However, it is an interesting observation,
that the magnetic equation of state from the present DSE computation without
the inclusion of the potential critical dynamics agrees quantitatively with
that in Braun  _et al._ (2020a), where the potential critical dynamics is
taken care of for $m_{\pi}\gtrsim 30$ MeV. This hints at the quantitative
irrelevance of the critical dynamics for a large range of pion masses.

More importantly, we have argued that for phenomenological applications to
heavy ion physics in the vicinity of the chiral phase transitions the presence
or absence of a large critical regime is not important. Instead, it is the
presence of quasi-massless modes which allows for the use of chiral transport
models or chiral fluid dynamics, e.g. Bluhm  _et al._ (2019, 2020); Grossi
_et al._ (2021); Florio  _et al._ (2021). These transport models require QCD
input, most prominently in terms of the full dispersion, see Bluhm  _et al._
(2019) for first steps in this direction, and the full interaction. In this
paper we have contributed to the computation of the latter in computing the
full chiral order parameter potential for the first time within the DSE
approach, see Section V.2, and in particular Figure 8 and Figure 9. This
computation utilises the light current quark mass dependence of the chiral
condensate, whose knowledge also allows us to estimated the size of the regime
with quasi-massless modes in Section V, see in particular Figure 13 and the
respective discussion. Importantly, we find a large regime with quasi-massless
modes. In particular we can identify the validity regime of vacuum $\chi$PT
(red regime in Figure 13) and a regime about the phase transition temperature
in the chiral limit, $T_{c0}$ (green regime in Figure 13, the respective
potential is shown in Figure 10) which allow for a direct expansion about the
chiral limit. Finally, we can identify a large regime with a Ginsburg Landau
type order parameter potential (blue regime in Figure 13), which also should
allow for chiral expansions. The asymptotic regime with large current quark
masses in indicated with grey regime in Figure 13.

The present results provide highly welcome support for the current use of
chiral transport models and chiral fluid dynamics in heavy ion collisions at
small densities. Currently, we refine the present qualitative analysis, as
well as extending it to the high density regime including the potential CEP.
We hope to report on the respective results in the near future.

Acknowledgements
We thank J. Braun, G. Eichmann, C. F. Fischer, W.-j. Fu, A. Y. Kotov, M. P.
Lombardo, J. Papavassiliou, F .Rennecke, B.-J. Schaefer, C. Schmidt, A. Trunin
and N. Wink for discussions. We also thank the other members of the fQCD
collaboration Braun  _et al._ (2021) for discussions and collaboration on
related subjects. F. Gao is supported by the Alexander von Humboldt
foundation. This work is supported by EMMI and the BMBF grant 05P18VHFCA. It
is part of and supported by the DFG Collaborative Research Centre SFB 1225
(ISOQUANT) and the DFG under Germany’s Excellence Strategy EXC - 2181/1 -
390900948 (the Heidelberg Excellence Cluster STRUCTURES).

## References

  * Luo and Xu (2017) X. Luo and N. Xu, Nucl. Sci. Tech. 28, 112 (2017), arXiv:1701.02105 [nucl-ex] .
  * Adamczyk  _et al._ (2017) L. Adamczyk _et al._ (STAR), Phys. Rev. C96, 044904 (2017), arXiv:1701.07065 [nucl-ex] .
  * Andronic  _et al._ (2018) A. Andronic, P. Braun-Munzinger, K. Redlich, and J. Stachel, Nature 561, 321 (2018), arXiv:1710.09425 [nucl-th] .
  * Stephanov (2006) M. A. Stephanov,  _Proceedings, 24th International Symposium on Lattice Field Theory (Lattice 2006): Tucson, USA, July 23-28, 2006_ , PoS LAT2006, 024 (2006), arXiv:hep-lat/0701002 [hep-lat] .
  * Andersen  _et al._ (2016) J. O. Andersen, W. R. Naylor, and A. Tranberg, Rev. Mod. Phys. 88, 025001 (2016), arXiv:1411.7176 [hep-ph] .
  * Shuryak (2017) E. Shuryak, Rev. Mod. Phys. 89, 035001 (2017), arXiv:1412.8393 [hep-ph] .
  * Pawlowski (2014) J. M. Pawlowski, Nucl. Phys. A931, 113 (2014).
  * Roberts and Schmidt (2000) C. D. Roberts and S. M. Schmidt, Prog. Part. Nucl. Phys. 45, S1 (2000), arXiv:nucl-th/0005064 [nucl-th] .
  * Fischer (2019) C. S. Fischer, Prog. Part. Nucl. Phys. 105, 1 (2019), arXiv:1810.12938 [hep-ph] .
  * Yin (2018) Y. Yin, (2018), arXiv:1811.06519 [nucl-th] .
  * Ding  _et al._ (2019) H. T. Ding _et al._ , Phys. Rev. Lett. 123, 062002 (2019), arXiv:1903.04801 [hep-lat] .
  * Braun  _et al._ (2020a) J. Braun, W.-j. Fu, J. M. Pawlowski, F. Rennecke, D. Rosenblüh, and S. Yin, Phys. Rev. D 102, 056010 (2020a), arXiv:2003.13112 [hep-ph] .
  * Kotov  _et al._ (2020) A. Y. Kotov, M. P. Lombardo, and A. M. Trunin, Eur. Phys. J. A 56, 203 (2020), arXiv:2004.07122 [hep-lat] .
  * Braun and Klein (2008) J. Braun and B. Klein, Phys. Rev. D 77, 096008 (2008), arXiv:0712.3574 [hep-th] .
  * Braun and Klein (2009) J. Braun and B. Klein, Eur. Phys. J. C 63, 443 (2009), arXiv:0810.0857 [hep-ph] .
  * Braun  _et al._ (2011a) J. Braun, B. Klein, and P. Piasecki, Eur. Phys. J. C 71, 1576 (2011a), arXiv:1008.2155 [hep-ph] .
  * Klein (2017) B. Klein, Phys. Rept. 707-708, 1 (2017), arXiv:1710.05357 [hep-ph] .
  * Gao and Pawlowski (2020) F. Gao and J. M. Pawlowski, Phys. Rev. D 102, 034027 (2020), arXiv:2002.07500 [hep-ph] .
  * Gao and Pawlowski (2021) F. Gao and J. M. Pawlowski, Phys. Lett. B 820, 136584 (2021), arXiv:2010.13705 [hep-ph] .
  * Gao  _et al._ (2021a) F. Gao, J. Papavassiliou, and J. M. Pawlowski, Phys. Rev. D 103, 094013 (2021a), arXiv:2102.13053 [hep-ph] .
  * Braun (2009) J. Braun, Eur. Phys. J. C64, 459 (2009), arXiv:0810.1727 [hep-ph] .
  * Braun  _et al._ (2011b) J. Braun, L. M. Haas, F. Marhauser, and J. M. Pawlowski, Phys. Rev. Lett. 106, 022002 (2011b), arXiv:0908.0008 [hep-ph] .
  * Fister and Pawlowski (2011) L. Fister and J. M. Pawlowski, (2011), arXiv:1112.5440 [hep-ph] .
  * Mitter  _et al._ (2015) M. Mitter, J. M. Pawlowski, and N. Strodthoff, Phys. Rev. D91, 054035 (2015), arXiv:1411.7978 [hep-ph] .
  * Braun  _et al._ (2016) J. Braun, L. Fister, J. M. Pawlowski, and F. Rennecke, Phys. Rev. D94, 034016 (2016), arXiv:1412.1045 [hep-ph] .
  * Rennecke (2015) F. Rennecke, Phys. Rev. D92, 076012 (2015), arXiv:1504.03585 [hep-ph] .
  * Fu  _et al._ (2016) W.-j. Fu, J. M. Pawlowski, F. Rennecke, and B.-J. Schaefer, Phys. Rev. D94, 116020 (2016), arXiv:1608.04302 [hep-ph] .
  * Rennecke and Schaefer (2017) F. Rennecke and B.-J. Schaefer, Phys. Rev. D96, 016009 (2017), arXiv:1610.08748 [hep-ph] .
  * Cyrol  _et al._ (2016) A. K. Cyrol, L. Fister, M. Mitter, J. M. Pawlowski, and N. Strodthoff, Phys. Rev. D94, 054005 (2016), arXiv:1605.01856 [hep-ph] .
  * Cyrol  _et al._ (2018a) A. K. Cyrol, M. Mitter, J. M. Pawlowski, and N. Strodthoff, Phys. Rev. D97, 054006 (2018a), arXiv:1706.06326 [hep-ph] .
  * Cyrol  _et al._ (2018b) A. K. Cyrol, M. Mitter, J. M. Pawlowski, and N. Strodthoff, Phys. Rev. D97, 054015 (2018b), arXiv:1708.03482 [hep-ph] .
  * Fu  _et al._ (2020a) W.-j. Fu, J. M. Pawlowski, and F. Rennecke, SciPost Phys. Core 2, 002 (2020a), arXiv:1808.00410 [hep-ph] .
  * Fu  _et al._ (2020b) W.-j. Fu, J. M. Pawlowski, and F. Rennecke, Phys. Rev. D 101, 054032 (2020b), arXiv:1909.02991 [hep-ph] .
  * Leonhardt  _et al._ (2020) M. Leonhardt, M. Pospiech, B. Schallmo, J. Braun, C. Drischler, K. Hebeler, and A. Schwenk, Phys. Rev. Lett. 125, 142502 (2020), arXiv:1907.05814 [nucl-th] .
  * Braun  _et al._ (2020b) J. Braun, M. Leonhardt, and M. Pospiech, Phys. Rev. D 101, 036004 (2020b), arXiv:1909.06298 [hep-ph] .
  * Qin  _et al._ (2011) S.-x. Qin, L. Chang, H. Chen, Y.-x. Liu, and C. D. Roberts, Phys. Rev. Lett. 106, 172301 (2011), arXiv:1011.2876 [nucl-th] .
  * Fischer  _et al._ (2011) C. S. Fischer, J. Luecker, and J. A. Mueller, Phys. Lett. B702, 438 (2011), arXiv:1104.1564 [hep-ph] .
  * Fischer and Mueller (2011) C. S. Fischer and J. A. Mueller, Phys. Rev. D84, 054013 (2011), arXiv:1106.2700 [hep-ph] .
  * Fischer  _et al._ (2014a) C. S. Fischer, L. Fister, J. Luecker, and J. M. Pawlowski, Phys. Lett. B732, 273 (2014a), arXiv:1306.6022 [hep-ph] .
  * Fischer  _et al._ (2014b) C. S. Fischer, J. Luecker, and C. A. Welzbacher, Phys. Rev. D90, 034022 (2014b), arXiv:1405.4762 [hep-ph] .
  * Eichmann  _et al._ (2016a) G. Eichmann, C. S. Fischer, and C. A. Welzbacher, Phys. Rev. D93, 034013 (2016a), arXiv:1509.02082 [hep-ph] .
  * Gao  _et al._ (2016) F. Gao, J. Chen, Y.-X. Liu, S.-X. Qin, C. D. Roberts, and S. M. Schmidt, Phys. Rev. D93, 094019 (2016), arXiv:1507.00875 [nucl-th] .
  * Gao and Liu (2016) F. Gao and Y.-x. Liu, Phys. Rev. D94, 076009 (2016), arXiv:1607.01675 [hep-ph] .
  * Gunkel  _et al._ (2019) P. J. Gunkel, C. S. Fischer, and P. Isserstedt, Eur. Phys. J. A55, 169 (2019), arXiv:1907.08110 [hep-ph] .
  * Isserstedt  _et al._ (2019) P. Isserstedt, M. Buballa, C. S. Fischer, and P. J. Gunkel, Phys. Rev. D100, 074011 (2019), arXiv:1906.11644 [hep-ph] .
  * Reinosa  _et al._ (2015) U. Reinosa, J. Serreau, and M. Tissier, Phys. Rev. D92, 025021 (2015), arXiv:1504.02916 [hep-th] .
  * Reinosa  _et al._ (2017) U. Reinosa, J. Serreau, M. Tissier, and A. Tresmontant, Phys. Rev. D95, 045014 (2017), arXiv:1606.08012 [hep-th] .
  * Maelger  _et al._ (2018a) J. Maelger, U. Reinosa, and J. Serreau, Phys. Rev. D97, 074027 (2018a), arXiv:1710.01930 [hep-ph] .
  * Maelger  _et al._ (2018b) J. Maelger, U. Reinosa, and J. Serreau, Phys. Rev. D98, 094020 (2018b), arXiv:1805.10015 [hep-th] .
  * Maelger  _et al._ (2020) J. Maelger, U. Reinosa, and J. Serreau, Phys. Rev. D 101, 014028 (2020), arXiv:1903.04184 [hep-th] .
  * Aguilar  _et al._ (2017) A. C. Aguilar, J. C. Cardona, M. N. Ferreira, and J. Papavassiliou, Phys. Rev. D 96, 014029 (2017), arXiv:1610.06158 [hep-ph] .
  * Aguilar  _et al._ (2018a) A. C. Aguilar, D. Binosi, C. T. Figueiredo, and J. Papavassiliou, Eur. Phys. J. C78, 181 (2018a), arXiv:1712.06926 [hep-ph] .
  * Aguilar  _et al._ (2018b) A. C. Aguilar, J. C. Cardona, M. N. Ferreira, and J. Papavassiliou, Phys. Rev. D 98, 014002 (2018b), arXiv:1804.04229 [hep-ph] .
  * Bazavov  _et al._ (2012) A. Bazavov _et al._ , Phys. Rev. Lett. 109, 192302 (2012), arXiv:1208.1220 [hep-lat] .
  * Borsanyi  _et al._ (2013) S. Borsanyi, Z. Fodor, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo, Phys. Rev. Lett. 111, 062005 (2013), arXiv:1305.5161 [hep-lat] .
  * Borsanyi  _et al._ (2014) S. Borsanyi, Z. Fodor, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo, Phys. Rev. Lett. 113, 052301 (2014), arXiv:1403.4576 [hep-lat] .
  * Bonati  _et al._ (2015) C. Bonati, M. D’Elia, M. Mariti, M. Mesiti, F. Negro, and F. Sanfilippo, Phys. Rev. D92, 054503 (2015), arXiv:1507.03571 [hep-lat] .
  * Bellwied  _et al._ (2015) R. Bellwied, S. Borsanyi, Z. Fodor, J. Guenther, S. D. Katz, C. Ratti, and K. K. Szabo, Phys. Lett. B751, 559 (2015), arXiv:1507.07510 [hep-lat] .
  * Bazavov  _et al._ (2017a) A. Bazavov _et al._ , Phys. Rev. D95, 054504 (2017a), arXiv:1701.04325 [hep-lat] .
  * Bazavov  _et al._ (2017b) A. Bazavov _et al._ (HotQCD), Phys. Rev. D96, 074510 (2017b), arXiv:1708.04897 [hep-lat] .
  * Bonati  _et al._ (2018) C. Bonati, M. D’Elia, F. Negro, F. Sanfilippo, and K. Zambello, Phys. Rev. D98, 054510 (2018), arXiv:1805.02960 [hep-lat] .
  * Borsanyi  _et al._ (2018) S. Borsanyi, Z. Fodor, J. N. Guenther, S. K. Katz, K. K. Szabo, A. Pasztor, I. Portillo, and C. Ratti, JHEP 10, 205 (2018), arXiv:1805.04445 [hep-lat] .
  * Bazavov  _et al._ (2019) A. Bazavov _et al._ (HotQCD), Phys. Lett. B795, 15 (2019), arXiv:1812.08235 [hep-lat] .
  * Guenther  _et al._ (2018) J. N. Guenther, S. Borsanyi, Z. Fodor, S. K. Katz, K. K. Szab, A. Pasztor, I. Portillo, and C. Ratti,  _Proceedings, 34th Winter Workshop on Nuclear Dynamics (WWND 2018): Guadeloupe, French West Indies, March 25-31, 2018_ , J. Phys. Conf. Ser. 1070, 012002 (2018).
  * Williams (2015) R. Williams, Eur. Phys. J. A51, 57 (2015), arXiv:1404.2545 [hep-ph] .
  * Williams  _et al._ (2016) R. Williams, C. S. Fischer, and W. Heupel, Phys. Rev. D93, 034026 (2016), arXiv:1512.00455 [hep-ph] .
  * Aoki  _et al._ (2020) S. Aoki _et al._ (Flavour Lattice Averaging Group), Eur. Phys. J. C 80, 113 (2020), arXiv:1902.08191 [hep-lat] .
  * Gasser and Leutwyler (1984) J. Gasser and H. Leutwyler, Annals Phys. 158, 142 (1984).
  * Bender  _et al._ (1998) A. Bender, G. I. Poulis, C. D. Roberts, S. M. Schmidt, and A. W. Thomas, Phys. Lett. B431, 263 (1998), arXiv:nucl-th/9710069 [nucl-th] .
  * Gao and Liu (2018) F. Gao and Y.-x. Liu, Phys. Rev. D97, 056011 (2018), arXiv:1702.01420 [hep-ph] .
  * Fischer (2006) C. S. Fischer, J. Phys. G32, R253 (2006), arXiv:hep-ph/0605173 [hep-ph] .
  * Bashir  _et al._ (2012) A. Bashir, L. Chang, I. C. Cloet, B. El-Bennich, Y.-X. Liu, C. D. Roberts, and P. C. Tandy, Commun. Theor. Phys. 58, 79 (2012), arXiv:1201.3366 [nucl-th] .
  * Eichmann  _et al._ (2016b) G. Eichmann, H. Sanchis-Alepuz, R. Williams, R. Alkofer, and C. S. Fischer, Prog. Part. Nucl. Phys. 91, 1 (2016b), arXiv:1606.09602 [hep-ph] .
  * Chen  _et al._ (2018) M. Chen, M. Ding, L. Chang, and C. D. Roberts, Phys. Rev. D 98, 091505 (2018), arXiv:1808.09461 [nucl-th] .
  * Dupuis  _et al._ (2021) N. Dupuis, L. Canet, A. Eichhorn, W. Metzner, J. M. Pawlowski, M. Tissier, and N. Wschebor, Phys. Rept. 910, 1 (2021), arXiv:2006.04853 [cond-mat.stat-mech] .
  * Gao  _et al._ (2021b) F. Gao, J. M. Pawlowski, and C. Schneider, in preparation (2021b).
  * Kaczmarek  _et al._ (2020a) O. Kaczmarek, F. Karsch, A. Lahiri, and C. Schmidt, in  _Criticality in QCD and the Hadron Resonance Gas_ (2020) arXiv:2010.15593 [hep-lat] .
  * Kaczmarek  _et al._ (2020b) O. Kaczmarek, F. Karsch, A. Lahiri, L. Mazur, and C. Schmidt (2020) arXiv:2003.07920 [hep-lat] .
  * Karsch (2019) F. Karsch, PoS CORFU2018, 163 (2019), arXiv:1905.03936 [hep-lat] .
  * Kotov  _et al._ (2021) A. Y. Kotov, M. P. Lombardo, and A. Trunin, Phys. Lett. B 823, 136749 (2021), arXiv:2105.09842 [hep-lat] .
  * Braun  _et al._ (2006) J. Braun, B. Klein, H. J. Pirner, and A. H. Rezaeian, Phys. Rev. D 73, 074010 (2006), arXiv:hep-ph/0512274 .
  * Bluhm  _et al._ (2019) M. Bluhm, Y. Jiang, M. Nahrgang, J. M. Pawlowski, F. Rennecke, and N. Wink,  _Proceedings, 27th International Conference on Ultrarelativistic Nucleus-Nucleus Collisions (Quark Matter 2018): Venice, Italy, May 14-19, 2018_ , Nucl. Phys. A982, 871 (2019), arXiv:1808.01377 [hep-ph] .
  * Bluhm  _et al._ (2020) M. Bluhm, M. Nahrgang, and J. M. Pawlowski, (2020), arXiv:2004.08608 [nucl-th] .
  * Grossi  _et al._ (2021) E. Grossi, A. Soloviev, D. Teaney, and F. Yan, Phys. Rev. D 104, 034025 (2021), arXiv:2101.10847 [nucl-th] .
  * Florio  _et al._ (2021) A. Florio, E. Grossi, A. Soloviev, and D. Teaney, (2021), arXiv:2111.03640 [hep-lat] .
  * Isserstedt  _et al._ (2021) P. Isserstedt, C. S. Fischer, and T. Steinert, Phys. Rev. D 103, 054012 (2021), arXiv:2012.04991 [hep-ph] .
  * fQCD collaboration (2021) fQCD collaboration, in preparation (2021).
  * Braun  _et al._ (2021) J. Braun, Y.-r. Chen, W.-j. Fu, F. Gao, F. Ihssen, A. Geissel, J. Horak, C. Huang, J. M. Pawlowski, F. Rennecke, F. Sattler, B. Schallmo, C. Schneider, Y.-y. Tan, S. Töpfel, R. Wen, J. Wessely, N. Wink, and S. Yin, (2021).