# Hot QCD matter around the chiral crossover:
a lattice study with $\mathcal{O}(a)$-improved Wilson fermions

Ardit Krasniqi  PRISMA+ Cluster of Excellence & Institut für Kernphysik,
Johannes Gutenberg-Universität Mainz, D-55099 Mainz, Germany  Marco Cè
Dipartimento di Fisica, Universitá di Milano-Bicocca and INFN, Sezione di
Milano-Bicocca,
Piazza della Scienza 3, 20126 Milano, Italy  Renwick J. Hudspith  GSI
Helmholtzzentrum für Schwerionenforschung,
64291 Darmstadt, Germany  Harvey B. Meyer  PRISMA+ Cluster of Excellence &
Institut für Kernphysik, Johannes Gutenberg-Universität Mainz, D-55099 Mainz,
Germany  Helmholtz Institut Mainz, Johannes Gutenberg-Universität Mainz,
D-55099 Mainz, Germany

###### Abstract

Using lattice QCD simulations with $\mathcal{O}(a)$-improved Wilson quarks and
physical up, down and strange quark masses, we investigate the properties of
thermal QCD matter at the temperatures $T=\\{128,154,192\\}$ MeV with a fixed
lattice spacing $a=0.064$ fm and volume $V=(6.1\,\text{fm})^{3}$. We find that
the pion quasiparticle, defined as the low-energy pole in the two-point
function of the axial charge, becomes lighter as the temperature increases and
give an argument based on hydrodynamics as to why the pole becomes purely
diffusive above the chiral crossover. We study the thermal modification of the
isovector vector spectral function using the Backus-Gilbert method, finding an
enhancement at low energies and a depletion at energies around 1 GeV. The
analogous study of the axial-vector channel reveals a larger enhancement at
energies below 1 GeV, and we show that these findings are consistent with
rigorous spectral sum rules. The difference between vector and axial-vector
correlators, an order parameter for chiral symmetry, turns out to be overall
suppressed by more than an order of magnitude at the crossover.

††preprint: MITP-24-055

##  I Introduction

Quark matter under extreme conditions of temperature and density is
interesting both from the experimental and theoretical point of view. In the
early universe, on a time scale of microseconds, the strongly-interacting
constituents (quarks and gluons) were in a hot and dense phase called the
Quark-Gluon Plasma (QGP). Heavy ion colliders, like the Large Hadron Collider
(LHC) at CERN and the Relativistic Heavy Ion Collider (RHIC) at BNL enable
similar conditions to be reached in the lab. As a result of its expansion, the
universe gradually cooled down, undergoing a transition to a hadronic phase in
which we now find ourselves. Note that at physical quark masses the transition
is actually a crossover Aoki:2006we  characterized by a pseudocritical
temperature $T_{pc}=156.5(1.5)$ MeV HotQCD:2018 ; see also the earlier results
of Ref. Borsanyi:2010bp .

In the limit of massless quarks, the QCD Lagrangian has a global
$SU(N_{f})_{L}\times SU(N_{f})_{R}$ symmetry corresponding to independent
rotations of the left- and right-handed components of the Dirac fields. This
symmetry is spontaneously broken to $SU(N_{f})_{V}$ at low temperatures and
becomes restored in the high-temperature phase – chiral symmetry restoration.
A non-vanishing value of the scalar density $\langle\bar{\psi}\psi\rangle(T)$
in the chiral limit characterizes the low-temperature phase ($0\leq T\leq
T_{c}\approx 132\,\text{MeV}$ HotQCD:2019 ). By contrast,
$\langle\bar{\psi}\psi\rangle(T)=0$ for $T>T_{c}$, indicating that chiral
symmetry is restored. In other words, the quark condensate
$\langle\bar{\psi}\psi\rangle$ is an order parameter for chiral symmetry
breaking.

The most immediate consequence of a broken chiral symmetry is the existence of
light pions, as implied by Goldstone’s theorem. Thus, starting from the QCD
vacuum, increasing the temperature initially leads to a dilute gas of pions.
As the temperature is further increased, other hadron species also begin to
contribute. At the same time, one expects the excitations of the medium to be
quasiparticles with somewhat modified properties as compared to the standard
hadrons, which are excitations of the vacuum. A natural starting point in the
investigation of the medium’s quasiparticles is to examine the properties of
the pion in the thermal environment Shuryak:1990ie ; Goity:1989gs . The pion
mass and decay constant have been studied to one loop in a thermal chiral
perturbation theory (ChPT) approach Gasser:1987ah . Additionally, the energy
density, the pressure and the quark condensate have been investigated up to
$\mathcal{O}(p^{8})$ in a chiral expansion below the phase transition
Gerber:1988tt . In Ref. Schenk:1993ru  the shift in the pion pole was
calculated as a function of temperature up to second order in the density.
Toublan Toublan:1997rr  calculated also the pion decay constant within thermal
ChPT to two loops and additionally examined the validity of the Gell-
Mann–Oakes–Renner relation at finite temperature. It is not _a priori_ clear
how far up in the temperature this expansion is applicable, since the
partition function is certainly no longer dominated by the pions for $T\gtrsim
100$ MeV Shuryak:1990ie ; Gerber:1988tt . However, the Goldstone-boson nature
of pions guarantees the existence of a divergent static correlation length for
vanishing quark masses Pisarski:1996mt .

In standard thermal ChPT, the quark mass as well as the temperature are
treated as small parameters, resulting in an expansion around $m_{q}=0$ and
$T=0$. In Refs. Son:2001ff ; Son:2002ci , however, Son and Stephanov
investigated perturbations only around $m_{q}=0$, keeping the temperature $T$
fixed to any value in the chirally broken phase. Although an explicit relation
of parameters like the quark condensate and pion decay constant to their $T=0$
counterparts is no longer possible in this framework, the validity of their
results is extended to a regime where neither ChPT nor perturbative QCD is
usable. Since lattice simulations rely on the imaginary-time formalism,
extracting real-time observables such as ‘pole masses’ out of lattice
quantities (e.g. Euclidean correlators) is in general a challenging task.
Nevertheless, at small $u,d$ quark masses the real part of the dispersion
relation of soft pions can be obtained from static Euclidean correlators, with
corrections that vanish in the chiral limit. What happens to the pion pole as
the system crosses over into the chirally restored phase is an aspect that we
investigate both numerically and with guidance from hydrodynamics.

More generally, changes in the spectrum of QCD excitations are an important
aspect of the phenomenology of heavy ion collisions. In particular, the
strongly interacting matter can produce lepton pairs via timelike photons, and
the invariant-mass spectrum of the latter gives access to the vector spectral
function of the thermal medium McLerran:1984ay  via a Kubo-Martin-Schwinger
relation. In the hadronic phase, dileptons are primarily produced through the
decay of vector mesons Braun-Munzinger:2015hba ; Rapp:1999ej . At sufficiently
high temperatures, one expects to be able to compute the vector spectral
function using QCD weak-coupling methods Laine:2013vma . Presently, it is
however not certain in what range of temperatures mesonic excitations
analogous to the $\rho,\omega$ and $\phi$ disappear from the medium. Here we
address this question in the isovector vector channel.

We note that, via a timelike $Z$ boson, dileptons can in principle also be
produced through the axial current Laine:2013vma , though at a very low rate
for $T\lesssim 1\,$GeV, which would probe the axial-vector spectral function
of the thermal medium. A restored chiral symmetry implies the vanishing of the
difference of the isovector vector and axial-vector spectral functions, $V-A$.
Determining in what temperature window this happens for physical quark masses
is one goal of this paper. In fact, we observe larger thermal effects in the
axial-vector spectral function than in the vector case. There are a number of
other probes of chiral symmetry restoration that can be studied. In Refs.
Aarts:2015mma ; Aarts:2017rrl  consequences of chiral symmetry restoration
have been investigated in the baryonic sector in the context of parity
doubling.

In the low-temperature phase, the dispersion relation of the pion
quasiparticle can be described at low momenta in terms of a renormalization-
group invariant (RGI) parameter $u(T)$. This temperature-dependent parameter
$u(T)$ is directly accessible on the lattice from screening quantities and
fulfills a dual role: not only does this parameter predict the modified
dispersion relation, but it also relates the static pion screening mass to the
(dynamic) pole mass. This has been done in Refs. Brandt:2014qqa ;
Brandt:2015sxa  and Refs. Ce:2022dax ; Krasniqi:2022djb  for the
$N_{\mathrm{f}}=2$ and $N_{\mathrm{f}}=2+1$ case, respectively. Here we extend
these studies to physical quark masses and higher statistical precision.

The paper is structured as follows: In Sec. II we start with the introduction
of some basic definitions which play a key role in the description of the pion
quasiparticle. We continue with the implementation of the lattice correlators,
followed by a brief description of the numerical setup. Our results on the
pseudoscalar sector, divided into subsections, are presented in Sec. III.
First we extract the mass and decay constant of the screening pion
(III.1-III.2). Next, we determine the pion velocity $u$ and examine to what
extent its extraction depends on a finite pion thermal width $\Gamma(T)$
(III.3-III.4). Subsequently, we reconstruct a smeared version of the spectral
function associated with the two-point function of the axial charge and
compare our results with the literature (III.5-III.6). We provide the static
pseudoscalar, vector and axial-vector screening masses along our temperature
scan in Sec. IV. Thereafter, we compare our lattice estimate for the isovector
quark number susceptibility with the prediction from the hadron resonance gas
model (HRG) in Sec. V. We look at several order parameters for chiral symmetry
restoration in Sec. VI, including the $V-A$ difference. Sec. VII contains our
results on the thermal modification of the vector and axial-vector spectral
functions. Finally, we collect our conclusions in Sec. VIII.

##  II Preliminaries

In this section we introduce the notation and some basic definitions as well
as the key quantities characterizing the pion quasiparticle that we will use
throughout the paper. Furthermore, the lattice implementation of the
correlators and the numerical setup are described briefly.

###  II.1 Definition of operators and correlation functions

The notation and conventions used in this work are adapted from Ref.
Ce:2022dax . Our framework is the up, down and strange-quark sector of
Euclideanized QCD on the space $S^{1}\times\mathbb{R}^{3}$, $S^{1}$ denoting
the Matsubara cycle of length $\beta\equiv 1/T$. We define the (isovector)
pseudoscalar density, the vector current and the axial-vector current as

$\displaystyle P^{a}(x)=\bar{\psi}(x)\gamma_{5}\frac{\tau^{a}}{2}\psi(x)\,,\ \
\ V^{a}_{\mu}(x)=\bar{\psi}(x)\gamma_{\mu}\frac{\tau^{a}}{2}\psi(x)\,,\ \
A^{a}_{\mu}(x)=\bar{\psi}(x)\gamma_{\mu}\gamma_{5}\frac{\tau^{a}}{2}\psi(x)\,,$
(1)

where $a\in\\{1,2,3\\}$ is an adjoint $SU(2)_{\text{isospin}}$ index,
$\tau^{a}$ is a Pauli matrix and
$\bar{\psi}(x)=(\bar{u}(x)\leavevmode\nobreak\ \bar{d}(x))$ is a Dirac field
flavor doublet. The partially-conserved axial current (PCAC) relation is an
operator identity that holds in Euclidean space when inserted in expectation
values. It relates the divergence of the axial vector current $A^{a}_{\mu}(x)$
to the pseudoscalar density $P^{a}(x)$,

$\displaystyle\partial_{\mu}A_{\mu}^{a}(x)=2m_{\text{PCAC}}\,P^{a}(x)\,.$ (2)

In the path integral formulation, this relation results from performing a
chiral rotation $\delta_{A}^{a}\psi(x)=\frac{\tau^{a}}{2}\gamma_{5}\psi(x)$ of
the fields (see Ref. Luscher:1998pe ). Applying the pseudoscalar density
operator on both sides and taking the expectation value one can obtain the
bare PCAC quark mass,

$\displaystyle 2m_{\text{PCAC}}{\langle
P^{b}(x)P^{b}(0)\rangle}={\partial_{\mu}\langle
A_{\mu}^{a}(x)P^{a}(0)\rangle}\,.$ (3)

It is convenient to integrate over three space-time directions, leaving only
one direction in which the derivative acts on the axial current. Since the
PCAC relation is an operator identity, we are free to choose this direction.
In our thermal system, the spatial direction is at least four times larger
than the temporal one. As a consequence, measuring along the spatial direction
results in a longer plateau and thus, smaller errors. Therefore, we will
extract the PCAC quark mass from the relation

$\displaystyle
m_{\text{PCAC}}(x_{3})=\frac{1}{2}\frac{\int\mathrm{d}x_{0}\mathrm{d}^{2}x_{\perp}\,\langle\widetilde{\partial}_{3}A_{3}^{a}(x)P^{a}(0)\rangle}{\int\mathrm{d}x_{0}\mathrm{d}^{2}x_{\perp}\,\langle
P^{b}(x)P^{b}(0)\rangle}\,,\ \ \ \ \ \ x_{\perp}=(x_{1},x_{2})\,.$ (4)

We introduce the static screening axial correlator, given by

$\displaystyle\delta^{ab}G_{A}^{s}(x_{3},T)=\int\,\mathrm{d}x_{0}\mathrm{d}^{2}x_{\perp}\langle
A_{3}^{a}(x)A_{3}^{b}(0)\rangle\overset{\absolutevalue{x_{3}}\rightarrow\infty}{=}\delta^{ab}\frac{f_{\pi}^{2}m_{\pi}}{2}e^{-m_{\pi}\absolutevalue{x_{3}}}\,,$
(5)

where we have specified the asymptotic form of the correlator, which defines
the pion screening mass $m_{\pi}$ and decay constant $f_{\pi}$. Analogously,
we define111Please note that in the case of identical operators at source and
sink, we do not use double labels for the correlators, e.g.
$G_{P}^{s}(x_{3},T)$, instead of $G_{PP}^{s}(x_{3},T)$ the following static
screening correlators:

$\displaystyle\delta^{ab}G_{P}^{s}(x_{3},T)$
$\displaystyle=\int\,\mathrm{d}x_{0}\mathrm{d}^{2}x_{\perp}\langle
P^{a}(x)P^{b}(0)\rangle$ (6) $\displaystyle\delta^{ab}G_{AP}^{s}(x_{3},T)$
$\displaystyle=\int\,\mathrm{d}x_{0}\mathrm{d}^{2}x_{\perp}\langle
A_{3}^{a}(x)P^{b}(0)\rangle\,.$ (7)

In order to probe the dynamical properties of the thermal system, we define
time-dependent correlators, projected to a definite spatial momentum,

$\displaystyle\delta^{ab}G_{A_{0}}(x_{0},T)$
$\displaystyle=\int\,\mathrm{d}^{3}x\,\langle A_{0}^{a}(x)A_{0}^{b}(0)\rangle$
(8) $\displaystyle\delta^{ab}G_{P}(x_{0},T)$
$\displaystyle=\int\,\mathrm{d}^{3}x\,\langle P^{a}(x)P^{b}(0)\rangle$ (9)
$\displaystyle\delta^{ab}G_{PA_{0}}(x_{0},T)$
$\displaystyle=\int\,\mathrm{d}^{3}x\,\langle
P^{a}(x)A_{0}^{b}(0)\rangle=-\int\,\mathrm{d}^{3}x\,\langle
A_{0}^{a}(x)P^{b}(0)\rangle$ (10) $\displaystyle\delta^{ab}G_{A}(x_{0},{\bf
p},T)$
$\displaystyle=-\frac{1}{3}\sum_{i=1}^{3}\int\,\mathrm{d}^{3}x\,e^{-i{\bf
p}\cdot{\bf x}}\,\langle A_{i}^{a}(x)A_{i}^{b}(0)\rangle$ (11)
$\displaystyle\delta^{ab}G_{V}(x_{0},{\bf p},T)$
$\displaystyle=-\frac{1}{3}\sum_{i=1}^{3}\int\,\mathrm{d}^{3}x\,e^{-i{\bf
p}\cdot{\bf x}}\,\langle V_{i}^{a}(x)V_{i}^{b}(0)\rangle\,.$ (12)

The time-dependent correlators are related to the corresponding spectral
function (see e.g. the review Meyer:2011gj ) via

$G_{X}(x_{0},\bm{p},T)=\int_{0}^{\infty}\,\mathrm{d}\omega\,\rho_{X}(\omega,{\bf
p},T)\,\frac{\text{cosh}(\omega(\beta/2-x_{0}))}{\text{sinh}(\omega\beta/2)}\,,\quad
X=A_{0},P,A,V.$ (13)

In Sec. III.5 we will analyze the axial spectral function using the Backus-
Gilbert method. The temperature argument of the correlators and spectral
functions will not always be shown explicitly.

###  II.2 Pion properties at finite temperature

It has been established within several frameworks Son:2001ff ; Son:2002ci
that at temperatures below the chiral phase transition a pion quasiparticle
persists, with the real part of the dispersion relation of sufficiently soft
pions given by

$\displaystyle\omega_{\mathbf{p}}=u(T)\sqrt{m_{\pi}^{2}+\mathbf{p}^{2}}\,,\ \
\ \text{for any}\ T\lesssim T_{c}\,.$ (14)

In the chiral limit, the parameter $u(T)$ can be interpreted as the group
velocity of a massless pion excitation. While the quasiparticle mass
$\omega_{\mathbf{0}}$ is the real-part of a pole of the retarded correlator
$G^{R}_{P}(\omega,{\bf p}=0,T)$ of the pseudoscalar density in the frequency
variable, the screening mass $m_{\pi}$ is a pole of
$G^{R}_{P}(\omega=0,\mathbf{p},T)$ in the spatial momentum
$\absolutevalue{\mathbf{p}}$ and represents an inverse spatial correlation
length. A simple interpretation of the dispersion relation (14) was given in
Ref. Brandt:2015sxa  in terms of the poles of the screening and the time-
dependent correlators. Son and Stephanov Son:2001ff ; Son:2002ci  showed that
the pion velocity $u$ in the chiral limit is given by the ratio of
$f_{\pi}^{2}$ to the axial charge susceptibility. The latter, however,
contains an ultraviolet divergence at any non-vanishing quark mass
Brandt:2014qqa  and is therefore not practical for lattice calculations. As an
alternative, in Refs. Brandt:2014qqa ; Brandt:2015sxa  the parameter $u$ was
determined using lattice correlation functions at vanishing spatial momentum
via two independent estimators,

$\displaystyle u_{m}$
$\displaystyle=\left[\left.-\frac{4m_{q}^{2}}{m_{\pi}^{2}}\,\frac{G_{P}(x_{0},T)}{G_{A_{0}}(x_{0},T)}\right|_{x_{0}=\beta/2}\right]^{1/2}\,,$
(15) $\displaystyle u_{f}$
$\displaystyle=\frac{f_{\pi}^{2}m_{\pi}}{2G_{A_{0}}(\beta/2,T)\,\text{sinh}(u_{f}m_{\pi}\beta/2)}\,,$
(16)

that agree in the chiral limit. We adopt them for the present study. In doing
so, for the estimator $u_{m}$, the parametric dominance of the pion in the
time-dependent Euclidean axial as well as the pseudoscalar density correlator
at small quark masses is exploited. The estimator $u_{f}$ only exploits the
pion dominance in the axial correlator; on the other hand, it relies on the
residue determined from the static screening correlator. The pion contribution
to the spectral function $\rho_{A_{0}}$ is expected to take the form of a
sharp peak,

$\rho_{A_{0}}(\omega,T)=\text{sgn}(\omega)\text{Res}(\omega_{\mathbf{0}})\delta(\omega^{2}-\omega_{\mathbf{0}}^{2})+\dots\,,$
(17)

where in Ref. Son:2002ci  (see also Ref. Brandt:2015sxa ) the residue was
predicted to have the form

$\text{Res}(\omega_{\mathbf{0}})\equiv(f_{\pi}^{t})^{2}\omega_{0}^{2}=f_{\pi}^{2}m_{\pi}^{2}\,,$
(18)

such that we can access the quasiparticle decay constant via
$f_{\pi}^{t}=f_{\pi}/u_{m}$.

We note that the equations of this subsection unambiguously define the
observables $\omega_{\bm{0}}$, $u_{m}$, $u_{f}$, $m_{\pi}$ and $f_{\pi}$, even
for temperatures above the chiral crossover, even though the notation reflects
their interpretation in the chirally broken phase. In the next subsection, we
elaborate on the expected pole structure at low energies in the axial-charge
correlator.

###  II.3 On the pole in the axial-charge correlator above the chiral
crossover

From the static pseudoscalar correlator, we define the overlap $G_{s}$ onto
the (screening) ground state,

$\int_{0}^{\beta}\mathrm{d}x_{0}\;\langle
P^{a}(x)\;P^{b}(0)\rangle\stackrel{{\scriptstyle|\bm{x}|\to\infty}}{{=}}-\delta^{ab}\,G_{s}^{2}\;\frac{e^{-m_{s}|\bm{x}|}}{4\pi|\bm{x}|}.$
(19)

Note that the engineering dimension of $G_{s}$ is GeV2. In this and the
following two paragraphs, we use Minkowski-space notation with mostly minus
metric $\eta^{\mu\nu}$ and $\\{\gamma^{\mu},\gamma^{\nu}\\}=2\eta^{\mu\nu}$.
In Ref. Son:2002ci , the spectral functions of the operators $A^{0,a}$ and
$P^{a}$ were determined in the chirally broken phase. Both are hydrodynamic
fields at small quark masses, $A^{0}$ because it is the density of the almost
conserved axial charge and $P$ because
$\varphi^{a}=2i\frac{P^{a}}{\langle\bar{\psi}\psi\rangle}$ corresponds to the
phase of the condensate [see e.g. Ref. Sharpe:2006pu , Sec. 2.2], which as a
Goldstone mode should be included in the list of hydrodynamic fields
Son:2002ci .

Here we are concerned with the chirally symmetric phase. Our main goal is to
determine the form of the spectral function of the axial charge density, which
is a hydrodynamic field, at low frequencies and momenta. The corresponding
current, up to higher-order terms, can be expressed via Fick’s law,

$\bm{A}^{a}=-D_{A}\nabla A^{0}+\dots$ (20)

In the chirally restored phase, the pseudoscalar density is no longer a
hydrodynamic field. Thus, it should be expressed via the axial charge density
and its spatial gradients. To leading order, we have the constitutive equation

$\phi^{a}\equiv iP^{a}=-\frac{\bar{\lambda}}{2}m_{q}G_{s}^{2}A^{0,a}+\dots,$
(21)

where the suppression by the quark mass reflects the fact that $P$ and $A^{0}$
belong to different chiral multiplets and as such decouple in the chiral limit
at the linear level. Note that $\bar{\lambda}$ is a hydrodynamic coefficient
and $(T^{5}\bar{\lambda})$ is renormalization-group invariant.

The PCAC relation can then be interpreted as an equation of motion for the
axial charge density,

$\partial_{0}A^{0,a}=-\nabla\cdot\bm{A}^{a}+2m_{q}\phi^{a}=D_{A}\triangle
A^{0,a}-{\bar{\lambda}}m_{q}^{2}G_{s}^{2}A^{0,a}+\dots\,.$ (22)

Linear response theory then leads to a prediction for the form of the spectral
function for the axial charge density (see e.g. Ref. Meyer:2011gj ),

$\frac{\rho_{A}(\omega,\bm{k})}{\omega}=\frac{\chi_{A}(\bm{k})}{\pi}\;\frac{\bar{\lambda}m_{q}^{2}G_{s}^{2}+D_{A}\bm{k}^{2}}{\omega^{2}+(\bar{\lambda}m_{q}^{2}G_{s}^{2}+D_{A}\bm{k}^{2})^{2}}.$
(23)

Close to the chiral limit, $\chi_{A}(\bm{k})\simeq\chi_{s}$ is close to the
isospin susceptibility corresponding to the vector current, and similarly,
$D_{A}=D_{I}$ corresponds to the diffusion coefficient of isospin. Thus, in
the chiral regime the only new and independent transport coefficient is
$\bar{\lambda}$. The qualitative difference in the spectral functions of the
broken phase and the restored phase at $\bm{k}=0$ is illustrated in Fig. 1. We
also note that for non-interacting quarks, $\bar{\lambda}=0$; see Appendix F.

Figure 1:  Qualitative features of the low-energy poles in the spectral
function of the axial charge, in the broken and in the restored phase.

From the result (23), one deduces the spectral funtions for the $\langle
A^{0}P\rangle$ and $\langle PP\rangle$ correlators at vanishing spatial
momentum,

$\displaystyle\rho_{AP}(\omega,\bm{0})$ $\displaystyle=$
$\displaystyle\frac{\omega}{2m_{q}}\rho_{A}(\omega,\bm{0})=\frac{\chi_{A}(0)}{2\pi}\;\frac{\bar{\lambda}m_{q}G_{s}^{2}\omega^{2}}{\omega^{2}+(\bar{\lambda}m_{q}^{2}G_{s}^{2})^{2}},$
(24) $\displaystyle\frac{\rho_{P}(\omega,\bm{0})}{\omega}$ $\displaystyle=$
$\displaystyle-\frac{1}{2m_{q}}\,\rho_{AP}(\omega,\bm{0})=-\frac{\chi_{A}(0)}{4\pi}\;\frac{\bar{\lambda}G_{s}^{2}\omega^{2}}{\omega^{2}+(\bar{\lambda}m_{q}^{2}G_{s}^{2})^{2}}.$
(25)

We note that in the chiral limit,

$\lim_{\omega\to
0}\frac{\rho_{P}(\omega,\bm{0})}{\omega}=-\frac{\chi_{s}}{4\pi}\;\bar{\lambda}\,G_{s}^{2}.$
(26)

By contrast, at large frequencies we have
$\frac{\rho_{P}(\omega,\bm{0})}{\omega}\stackrel{{\scriptstyle|\omega|\to\infty}}{{\sim}}-\frac{N_{c}}{16\pi^{2}}|\omega|$.
Eq. (26) represents a Kubo formula for the transport coefficient
$\bar{\lambda}$.

Finally, it is useful to record the parametric size of the Euclidean
correlators $G_{A_{0}}$ and $G_{P}$ at $x_{0}=\beta/2$; see Table 1. This
information determines whether the real part of the pole can be determined
without having to address an inverse problem. Clearly, this is possible in the
broken phase near the chiral limit, because both $G_{A_{0}}$ and $G_{P}$ are
dominated by the pole contribution. In the restored phase, however, the
pseudoscalar correlator is not dominated by this contribution.

pion pole rest diffusion pole rest $G_{A_{0}}$ O(1) O($m_{q}^{2}$) O(1)
O($m_{q}^{2}$) $G_{P}$ O($m_{q}^{-1}$) O(1) O(1) O(1)

Table 1:  Parametric size, in terms of the average up-down quark mass $m_{q}$,
of the low-energy contribution vs. the rest in the axial charge density and in
the pseudoscalar density Euclidean-time dependent correlators at zero spatial
momentum and $x_{0}=\beta/2$. The left part of the table concerns the chirally
broken phase, and the right-hand part the restored phase.

###  II.4 Lattice implementation of the correlators

In this work we use exclusively the local discretizations of the operators
introduced in the previous subsection. Therefore, the expression of the bare
operators in the lattice theory is the same as in Eq. (1). These bare
operators are first O($a$)-improved and then renormalized. While the bare
pseudoscalar density is by itself O($a$)-improved, the details of the
improvement of the vector and axial-vector currents are given in Appendix A.

The finite renormalization of the vector and the axial-vector currents is
performed with the non-perturbatively determined renormalization factors
$Z_{V}(g_{0}^{2})$ and $Z_{A}(g_{0}^{2})$, supplemented by a quark-mass
dependent factor in order to fully realize O($a$) improvement; details are
provided in Appendix B.

The pseudoscalar density $P^{a}(x)$ acquires a scale (and scheme) dependence
via the process of renormalization. The renormalization factor is notated
$Z_{P}(g_{0}^{2},a\mu)$. Here, we renormalize $P^{a}(x)$ in the (non-
perturbative) gradient-flow (GF) scheme at the renormalization scale $\mu$
where the corresponding coupling $\bar{g}^{2}_{\rm GF}=9.25$; this corresponds
to a low scale of $\mu\approx 230\;$MeV Campos:2018ahf . While none of our
physics applications relies on the choice of a specific scheme, we note that
in the latter publication, the scale dependence of the renormalization factor
has been computed up to perturbative scales $\mu$; thereby the connection to
the renormalization-group invariant (RGI) operator is known.

The chosen renormalization of the PCAC mass preserves the axial Ward identity
Eq. (2). Thus, all renormalization-scale dependent quantities in this paper
are quoted in the aforementioned gradient-flow scheme. In particular the PCAC
mass is renormalized by multiplying it with $Z_{A}/Z_{P}$, and the combination
$m_{\pi}^{2}f_{\pi}^{2}/m_{\rm PCAC}$ considered in Sec. VI.1 via the factor
$Z_{A}Z_{P}$. The numerical values of the renormalization factors are
collected in Appendix B.

###  II.5 Numerical setup

Our calculations are performed on three $N_{\mathrm{f}}=2+1$ ensembles with
tree-level $\mathcal{O}(a^{2})$-improved Lüscher-Weisz gauge action and non-
perturbatively $\mathcal{O}(a)$-improved Wilson fermions Bulava:2013cta . The
action corresponds to the choice of the Coordinated Lattice Simulations (CLS)
initiative Bruno:2014jqa  and the bare parameters match those of the CLS zero-
temperature ensemble E250 Mohler:2017wnb . The latter are listed in Table 2,
together with the lattice spacing as determined in Ref. Bruno:2016plf . The
spatial extent is $N_{s}=96$ on all three ensembles and the time direction
admits thermal boundary conditions with $N_{\tau}\in\\{24,20,16\\}$, which is
the only difference relative to the zero-temperature ensemble, resulting in
the temperatures

$\displaystyle T_{24}$ $\displaystyle=\frac{1}{24a}=127.9(1.5)\,{\rm MeV}\,,$
(27) $\displaystyle T_{20}$ $\displaystyle=\frac{1}{20a}=153.5(1.8)\,{\rm
MeV}\,,$ (28) $\displaystyle T_{16}$
$\displaystyle=\frac{1}{16a}=191.9(2.3)\,{\rm MeV}\,.$ (29)

With a pseudocritical temperature $T_{pc}=156.5(1.5)$ MeV obtained in
$(2+1)$-flavor QCD HotQCD:2018 , our temperature scan covers the range
$\\{T_{24},T_{20},T_{16}\\}/T_{pc}\approx\\{0.82,0.98,1.23\\}$. Notably, these
temperature choices place one ensemble in the hadronic phase, another in the
high-temperature or quark-gluon plasma phase, and the third at the chiral
crossover.

For reference, we also quote the zero-temperature pseudoscalar masses and pion
decay constant, determined in Ref. Ce:2022kxy ,

$\displaystyle T=0:\qquad m_{\pi}^{0}$ $\displaystyle=130.0(0.7)(1.5)\,{\rm
MeV},\qquad m_{K}^{0}=489.36(0.25)(5.8)\,{\rm MeV},$ (30) $\displaystyle
f_{\pi}^{0}$ $\displaystyle=87.1(0.3)(1.0)\,{\rm MeV}$ (31)

where the first error is from the corresponding quantity in lattice units, and
the second is from the lattice spacing determination of Ref. Bruno:2016plf .

Table 2:  Parameters and lattice spacing of the ensembles analyzed in this work. The lattice spacing determination is from Ref. Bruno:2016plf .  $N_{\tau}/a$ |  $L/a$ |  $6/g_{0}^{2}$ |  $\kappa_{l}$ |  $\kappa_{s}$ |  $a\,[{\rm fm}]$ |  $N_{\text{conf}}$ |  $\frac{\text{MDUs}}{\text{conf}}$ | Label
---|---|---|---|---|---|---|---|---
24 | 96 | 3.55 | 0.137232867 | 0.136536633 | 0.06426(76) | 1200 | 4 | E250Nt24
20 | 96 | 3.55 | 0.137232867 | 0.136536633 | 0.06426(76) | 1000 | 4 | E250Nt20
16 | 96 | 3.55 | 0.137232867 | 0.136536633 | 0.06426(76) | 1500 | 4 | E250Nt16

The ensembles have been generated using version 2.0 of the openQCD package
Luscher:2012av , applying a small twisted mass to the light quark doublet for
algorithmic stability. The correct QCD expectation values are obtained
including the reweighting factors for the twisted mass and for the rational
approximation of the propagator used to simulate the strange quark
Clark:2006fx ; Luscher:2012av ; Mohler:2020txx ; Kuberski:2023zky . Errors are
estimated using the $\Gamma$ method Madras:1988ei ; Wolff:2003sm ;
Ramos:2018vgu ; Joswig:2022qfe .

Measurements are performed using stochastic wall sources in the implementation
of Renwick J. Hudspith’s Witnesser code. One improvement of this calculation
over our previous one Ce:2022dax  is the use of stochastic wall sources for
the two-point functions. Here we use $Z_{2}\otimes Z_{2}$ noise Foster:1998vw
, inverting on the spin indices only. These sources, also known as “Z2SEMWall”
Boyle:2008rh  or “linked” sources ETM:2008zte  only require 4 inversions of
the Dirac matrix as compared to the usual 12 for a point source. As the wall-
source has support on a large volume, it is less sensitive to local
fluctuations of the gauge field, and typically gives a smaller variance to
two-point functions at fixed cost in comparison to point sources. On top of
this we make heavy use of the Truncated Solver Method Bali:2009hu  with a
large number of low-precision solves, making the calculations presented here
fairly cheap to perform.

The results for the E250Nt24 ensemble presented in Ref. Ce:2022dax  were
obtained using point sources. Comparing the result for the pion screening
decay constant, $f_{\pi}/T=0.559(11)$, therein to our result quoted here,
$f_{\pi}/T=0.569(4)$, we have gained a factor of $\approx 8$ in statistics
[see Sec. III.2].

##  III Results on the pseudoscalar sector

In this section, we present our lattice results on observables in the
pseudoscalar sector, i.e. those related to pion properties. As an important
benchmark, we begin with the determination of the light quark doublet PCAC
masses.

###  III.1 The PCAC mass: a control quantity

Figure 2: Renormalized PCAC mass in the E250Nt24 ensemble along the $x_{0}$
and $x_{3}$ direction. The final result — obtained from a fit along the
$x_{3}$ direction — is also shown with a $1-\sigma$ band. We have used the
improved axial current together with the symmetrized derivative (see Eqs.
(78-76)).

Figure 3:  Renormalized PCAC masses in the E250Nt20 ensemble (left) and
E250Nt16 ensemble (right) along the $x_{0}$ and $x_{3}$ direction. For the
E250Nt20 ensemble we have used $\mathcal{O}(a)$-improvement of the axial
current and the symmetrized derivative (see Eqs. (78-76)). In contrast, the
PCAC mass on E250Nt16 (right) was determined using the unimproved version of
the axial current.

The extraction of the PCAC masses as defined in Eq. (4) is carried out by
performing a fit to a constant in the range where a plateau is observed; see
Figs. 2-3. Here we benefit from the long plateaus available in the $x_{3}$
direction. We obtain

$\displaystyle\begin{split}a\cdot m_{\text{PCAC}}^{vac}&=0.001428(16)\,,\\\
a\cdot m_{\text{PCAC}}^{N_{t}24}&=0.001466(41)\,,\\\ a\cdot
m_{\text{PCAC}}^{N_{t}16}&=0.001335(86)\,.\end{split}$ (32)

For the boxes in the hadronic phase and at the chiral phase transition the
extracted PCAC masses obtained from the $x_{0}$ direction are compatible with
the ones obtained from the $x_{3}$ direction, pointing to small cutoff effects
at this value of the lattice spacing. However, the time extent222Due to
periodic boundary conditions, we are dealing with effectively eight points in
time direction. of the high temperature E250Nt16 box is too small to match the
value of the PCAC mass obtained from the spatial direction. Furthermore, for
the E250Nt16 ensemble we have not used the improved axial current. The
$AP$-correlator is an order parameter for chiral symmetry restoration and
consequently approaches zero in the chirally symmetric phase. Thus, the
correction term proportional to $c_{A}$ in Eq. (78) becomes of the same order
of magnitude as the correlator itself in the chirally restored phase.

Since the PCAC mass is based on an operator identity, it should be independent
of the temperature. However, in lattice units the value of the PCAC mass,
$a\cdot m_{\text{PCAC}}^{N_{t}20}=0.001003(90)$, is off by a factor of
$\approx 1.5$ relative to the corresponding vacuum box. As the PCAC mass is
proportional to the derivative of the $AP$-correlator, we have compared the
integrated autocorrelation time $\tau_{int}$ at fixed source-sink seperation
$x_{3}/a=10$ on the E250Nt20 box with the E250Nt24 box. As can be seen in Fig.
4, we are dealing with severe autocorrelation in this observable at the chiral
crossover (More details in Appendix C). As a consequence, we are going to use
the value of the PCAC mass extracted from the E250Nt24 ensemble in our further
analysis, since it is compatible within errors with the PCAC mass on the
vacuum ensemble.

Figure 4:  Comparison of the integrated autocorrelation times $\tau_{int}$ for
the $AP$-correlator at fixed source sink seperation $x_{3}/a=10$. The l.h.s.
shows the correlator on the hadronic E250Nt24 ensemble and the r.h.s. shows
the correlator on the chiral crossover ensemble E250Nt20.

###  III.2 Static correlators: the pion screening mass and decay constant

In order to extract the pion screening mass and decay constant, we fit
simultaneously the static screening pseudoscalar and the axial-pseudoscalar
correlator In Appendix E, we provide the details of the calculation, including
effective mass plots.

The results are listed in Table 3 in units of the respective temperature. In
physical units, they correspond to

$\displaystyle m_{\pi}^{Nt24}$ $\displaystyle=$ $\displaystyle 145(2)\,{\rm
MeV}\,,$ $\displaystyle m_{\pi}^{Nt20}$ $\displaystyle=$ $\displaystyle
196(4)\,{\rm MeV}\,,$ (33) $\displaystyle m_{\pi}^{Nt16}$ $\displaystyle=$
$\displaystyle 555(8)\,{\rm MeV}\,.$

The pion screening decay constant is significantly reduced on all three
ensembles compared to the decay constant $f_{\pi}^{0}=87.4(1.0)$ MeV on the
corresponding zero-temperature ensemble. We obtain

$\displaystyle f_{\pi}^{Nt24}$ $\displaystyle=$ $\displaystyle 73(1)\,{\rm
MeV}\,,$ (34) $\displaystyle f_{\pi}^{Nt20}$ $\displaystyle=$ $\displaystyle
26(2)\,{\rm MeV}\,,$ (35) $\displaystyle f_{\pi}^{Nt16}$ $\displaystyle=$
$\displaystyle 5.9(4)\,{\rm MeV}\,.$ (36)

The observables characterizing the properties of the pion quasiparticle across
the chiral phase transition are summarized in Table 3.

###  III.3 Properties of the pion quasiparticle

The results for the estimators $u_{m}$ and $u_{f}$ defined in Eqs. (15,16)
together with those of the screening quantities are presented in Table 3. Good
agreement is found for the two independent estimators $u_{f}$ and $u_{m}$ of
the pion velocity, in the hadronic as well as high-temperature phase. Both of
them differ significantly from unity, which clearly represents a breaking of
Lorentz invariance due to thermal effects. However, well in the chirally
restored phase, the interpretation of $u$ as the pion velocity is no longer
valid: still, the ratio $u_{f}/u_{m}$ is expected to remain finite in the
chiral limit, since both $u_{f}$ and $u_{m}$ are of order the quark
mass333Indeed, $m_{\pi}$ is of order unity but $f_{\pi}$ is of order $m_{q}$.
In this point, we rectify an uncorrect remark written in Ref. Brandt:2014qqa
.. At the chiral crossover, $u_{f}$ and $u_{m}$ differ significantly,
$u_{f}/u_{m}\approx 0.64$; in the latter estimate, we use the PCAC mass
obtained from E250Nt24, since it is more reliable. Thus, the interpretation as
a pion velocity is already questionable at the chiral crossover.

Additionally, we found that the zero-temperature pion mass given in Eq. (30)
‘splits’ into a lower pion quasiparticle mass and a higher pion screening
mass. For the former, we obtain from the estimator
$\omega_{\bm{0}}=u_{m}\,m_{\pi}$

$\displaystyle\omega_{\mathbf{0}}^{Nt24}$ $\displaystyle=$ $\displaystyle
117(2)\,{\rm MeV}\,,$ (37) $\displaystyle\omega_{\mathbf{0}}^{Nt20}$
$\displaystyle=$ $\displaystyle 81(2)\,{\rm MeV}\,,$ (38)
$\displaystyle\omega_{\mathbf{0}}^{Nt16}$ $\displaystyle=$ $\displaystyle
26(1)\,{\rm MeV}\,;$ (39)

these values can be compared to the pion screening masses given in Eq. (33).
Clearly, the drop in $\omega_{\mathbf{0}}$ is particularly rapid between the
crossover point and our ensemble in the high-temperature phase.

The quasiparticle decay constant in the hadronic phase,
$f^{t}_{\pi}=f_{\pi}/u_{m}=90(1)$ MeV, is much closer to the vacuum decay
constant $f_{\pi}^{0}=87(1)$ MeV than $f_{\pi}^{Nt24}$. At the chiral
crossover, we obtain $f^{t}_{\pi}=62(3)$ MeV ($f^{t}_{\pi}=96(3)$ MeV using
$u_{m}$ determined from the PCAC mass on the E250Nt20 ensemble) and in the
high-temperature phase $f^{t}_{\pi}=127(10)$ MeV.

Table 3: Summary of the results of the E250 thermal ensembles with $N_{\tau}\in\\{24,20,16\\}$. The pion quasiparticle mass $\omega_{\mathbf{0}}$ is calculated using $\omega_{\mathbf{0}}=u_{m}m_{\pi}$. Ensemble | E250Nt24 | E250Nt20 | E250Nt16
---|---|---|---
$m_{\pi}/T$ | 1.136(16) | 1.250(24) | 2.891(43)
$f_{\pi}/T$ | 0.569(4) | 0.172(12) | 0.0305(21)
$u_{f}$ | 0.831(3) | 0.274(11) | 0.0462(31)
$u_{m}$ | 0.807(9) | 0.425(13) | 0.0462(16)
$u_{f}/u_{m}$ | 1.028(10) | 0.644(28) | 1.001(76)
$\omega_{\mathbf{0}}/T$ | 0.917(15) | 0.531(16) | 0.133(4)
$f_{\pi}^{t}/T$ | 0.705(7) | 0.405(18) | 0.661(50)
$\text{Res}(\omega_{\mathbf{0}})/T^{4}$ | 0.418(14) | 0.0463(60) | 0.0078(11)

###  III.4 Dependence of the pion velocity $u_{f}$ on a finite pion thermal
width $\Gamma(T)$

The analysis of Son and Stephanov Son:2001ff ; Son:2002ci  concluded that at
temperatures below the chiral phase transition, the imaginary part of the pion
pole is parametrically small compared to its real part. In this subsection, we
investigate the sensitivity of our results for the pion quasiparticle mass and
velocity parameter $u$ to the assumption of a negligible thermal width of this
quasiparticle. In order to examine the consequences of a finite thermal pion
width on the pion velocity $u$ we replace the $\delta$-distribution in Eq.
(17) by a Breit-Wigner peak of width $\Gamma(T)$ resulting in

$\rho_{A_{0}}(\omega,T)=f_{\pi}^{2}m_{\pi}^{2}\frac{\Gamma(T)}{\pi}\frac{1}{2\omega_{\mathbf{0}}}\left(\frac{1}{(\omega-\omega_{\mathbf{0}})^{2}+\Gamma(T)^{2}}-\frac{1}{(\omega+\omega_{\mathbf{0}})^{2}+\Gamma(T)^{2}}\right)+\dots\,,$
(40)

where the second term is needed to ensure the antisymmetry of the spectral
function in $\omega$ Meyer:2011gj . Expressing the correlator midpoint of the
time-dependent Euclidean correlator $G_{A_{0}}(\beta/2,T)$ in terms of the
spectral function $\rho_{A_{0}}$ with the help of Eq. (13) and using
$\omega_{\mathbf{0}}=u_{f}m_{\pi}$ one can extract the pion velocity $u_{f}$
for different thermal pion widths. The results are shown in Table 4.

Table 4: Dependence of the pion velocity $u_{f}$ on a finite pion thermal width $\Gamma(T)$ on the E250 thermal ensembles with $N_{\tau}\in\\{24,20,16\\}$. $\Gamma(T)\,[\text{MeV}]$ | E250Nt24 | E250Nt20 | E250Nt16
---|---|---|---
5 |  |  | 0.0454(31)
10 |  |  | 0.0426(33)
15 | 0.825(3) | 0.262(20) | 0.0375(38)
20 | 0.819(3) | 0.251(21) | 0.0289(50)
30 | 0.806(3) | 0.224(23) |
40 | 0.785(3) | 0.177(28) |
80 | 0.626(6) |  |
100 | 0.472(10) |  |

We find that the extracted estimator of the pion velocity in the hadronic
phase, $u_{f}=0.831(3)$ (assuming the presence of a discrete delta term in the
spectral function), is consistent with a Breit-Wigner approach up to pion
thermal widths $\Gamma(T)\approx 15\,\text{MeV}$. At the chiral crossover the
estimate assuming a vanishing width, $u_{f}=0.274(11)$, is consistent with the
estimation of $u_{f}$ assuming a pion thermal width of up to $\approx
20\,\text{MeV}$.

###  III.5 Spectral function reconstruction with the Backus-Gilbert method

In order to extract the spectral function $\rho_{A_{0}}(\omega)$ at zero
momentum from the corresponding temporal Euclidean correlator,
$G_{A_{0}}(\tau_{i},T),\ \tau\equiv x_{0}$, one has to invert the analogue for
$G_{A_{0}}$ of Eq. (13) with a kernel
$K(\tau_{i},\omega)=\text{cosh}(\omega(\beta/2-\tau_{i}))/\text{sinh}(\omega\beta/2)$,
encountering a numerically ill-posed problem. A possible approach to surmount
this is the Backus-Gilbert method Backus:1968 . We adopt the notation of Ref.
Brandt:2015sxa , where the method was applied to lattice QCD for the first
time. It should be emphasized that, with this approach, no particular ansatz
is needed for the spectral function. The Backus-Gilbert method provides an
estimator for the smeared axial spectral function, $\hat{\rho}$. It is useful
to introduce a rescaling function $f(\omega)$, with $\lim_{\omega\to
0}f(\omega)/\omega$ finite and $f(\omega)>0$ for $\omega>0$, and to write
$K^{f}(\tau,\omega)=f(\omega)K(\tau,\omega)$. The values of $\hat{\rho}$ are
then constructed linearly from the lattice correlator data $G(\tau_{i})$,

$\frac{\hat{\rho}(\bar{\omega})}{f(\bar{\omega})}=\sum_{i=1}^{N_{p}}q_{i}(\bar{\omega})G(\tau_{i})=\int_{0}^{\infty}\mathrm{d}\omega\,\hat{\delta}(\bar{\omega},\omega)\frac{\rho(\omega)}{f(\omega)}\,.$
(41)

Note that the coefficients $q_{i}$ depend on some reference value
$\bar{\omega}$ around which the resolution function,

$\hat{\delta}(\bar{\omega},\omega)=\sum_{i=1}^{N_{\tau}}q_{i}(\bar{\omega})K^{f}(\tau_{i},\omega)\,,$
(42)

is concentrated. The coefficients $q$ will be chosen such that the resolution
function is normalized according to

$\int_{0}^{\infty}\mathrm{d}\omega\,\hat{\delta}(\bar{\omega},\omega)=1\,.$
(43)

Although the final coefficients $q$ we use depend on the choice of $f$, and
the functions $\hat{\rho}$ and $\hat{\delta}$ depend on this choice as well,
we do not indicate this explicitly, in order to avoid overcluttering the
notation. The desirable resolution function would be a Dirac delta
distribution centered at $\bar{\omega}$. In order to make the resolution
function as sharply centered around $\bar{\omega}$ as possible, we follow the
strategy of minimizing the second moment of its square, subject to the
constraint in Eq. (43). However, in order to avoid widely oscillating
coefficients $q_{i}(\bar{\omega})$, the quadratic form to be minimized is
regulated by the covariance matrix $\text{Cov}[G]$ of the Euclidean
correlator: defining

$\displaystyle R_{i}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\mathrm{d}\omega\,K^{f}(\tau_{i},\omega),$ (44)
$\displaystyle W_{ij}(\bar{\omega})$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\mathrm{d}\omega\,K^{f}(\tau_{i},\omega)(\omega-\bar{\omega})^{2}K^{f}(\tau_{j},\omega),$
(45) $\displaystyle W^{\text{reg}}(\bar{\omega})$ $\displaystyle=$
$\displaystyle\lambda
W(\bar{\omega})+(1-\lambda)\text{Cov}[G]\phantom{\Big{|}}$ (46)

with $0<\lambda<1$, the function to be minimized is ${\cal
F}_{\alpha}(q)=(q,W^{\rm reg}q)-\alpha((q,R)-1)$, with the result

$\displaystyle q(\bar{\omega})$ $\displaystyle=$
$\displaystyle\frac{1}{(R,W^{\text{reg}}(\bar{\omega})^{-1}R)}\;W^{\text{reg}}(\bar{\omega})^{-1}R.$
(47)

For the analysis of $\rho_{A_{0}}$, we make the choice

$f(\omega):=\tanh(\beta\omega/2).$ (48)

The values of $\lambda$ quoted below refer to units in which all dimensionful
quantities are turned into dimensionless ones by appropriate powers of the
temperature. Some examples for the resolution function
$T\hat{\delta}(\bar{\omega},\omega)$ for different values of $\lambda$ are
shown in the left panel of Fig. 5. The right panel of Fig. 5 shows the smeared
and rescaled axial spectral function
$\hat{\rho}_{A_{0}}(\bar{\omega},T)/(f\,T^{2})$ with $\lambda=10^{-3}$ for the
hadronic E250Nt24 ensemble, while Fig. 6 shows the results on the E250Nt20
(left panel) and E250Nt16 ensemble (right panel). It demonstrates model-
independently that the axial-charge correlator is dominated by low
frequencies. On the other hand, for $\omega/T>8$ we observe the expected
perturbative slight increase proportional to the quark mass squared.
Furthermore, the predicted position of the quasiparticle mass
$\omega_{\mathbf{0}}$ is close to the peak of the smeared spectral function in
the hadronic regime.

Figure 5:  Left panel: Some examples of resolution functions for different
values of $\lambda$, centered around $\bar{\omega}/T$. Right panel: Estimator
of the spectral function
$\hat{\rho}_{A_{0}}(\omega,T=128\,\text{MeV})/(f\,T^{2})$. The blue dashed
line corresponds to the location of the expected position of the pole
$\omega_{\mathbf{0}}$ according to Eq. (14). The temporal
$G_{A_{0}}(x_{0},T)$-correlator is $\mathcal{O}(a)$-improved.

Figure 6:  Left panel: Estimator of the spectral function
$\hat{\rho}_{A_{0}}(\omega,T=154\,\text{MeV})/(f\,T^{2})$. The blue dashed
line corresponds to the location of the expected position of the pole
$\omega_{\mathbf{0}}$ according to Eq. (14). The temporal
$G_{A_{0}}(x_{0},T)$-correlator is $\mathcal{O}(a)$-improved. Right panel:
Estimator of the spectral function
$\hat{\rho}_{A_{0}}(\omega,T=192\,\text{MeV})/(f\,T^{2})$. The blue dashed
line corresponds to the location of the expected position of the pole
$\omega_{\mathbf{0}}$ according to Eq. (14).

Up until now, we have determined the pion quasiparticle mass
$\omega_{\mathbf{0}}$ indirectly with the aid of screening quantities and
compared it to the location of the pion pole position in the smeared axial
spectral funtion $\hat{\rho}_{A_{0}}$. As a next step, we want to predict the
temporal axial correlator $G_{A_{0}}(x_{0},T=128\,\text{MeV})$ in terms of the
screening quantities $f_{\pi}$ and $m_{\pi}$. Using Eqs. (13) and (17-18) we
obtain:

$G_{A_{0}}(x_{0},T)=\frac{\text{Res}(\omega_{\mathbf{0}})}{2\omega_{\mathbf{0}}}\frac{\text{cosh}(\omega_{\mathbf{0}}(\beta/2-x_{0}))}{\text{sinh}(\omega_{\mathbf{0}}\beta/2)}+\cdots.$
(49)

The result of this prediction together with the lattice data of the temporal
$G_{A_{0}}(x_{0},T)$-correlator and its fit result are displayed in the right
panel of Fig. 19 in Appendix G. As can be seen the prediction shows excellent
agreement with the lattice data for $x_{0}/a\geq 7$. The direct fit result
reads $\omega_{\mathbf{0}}/T=0.91(9)$. As a further crosscheck we also have
determined the effective cosh mass directly from the temporal correlator,
obtaining $m_{\text{cosh}}(x_{0})/T=0.83(6)$ [see left panel of Fig. 19].

We want to emphasize that all three estimators for the pion pole mass in the
hadronic phase, namely

  * •

$\omega_{\mathbf{0}}=u_{m}m_{\pi}=117(2)$ MeV, obtained from screening
quantities,

  * •

$\omega_{\mathbf{0}}=116(12)$ MeV, obtained from a direct fit of the temporal
correlator,

  * •

$m_{\text{cosh}}(x_{0})/T=106(8)$ MeV, obtained as the effective mass of the
temporal correlator,

are all consistent within errors of each other. They all give significantly
reduced values in comparison to the in vacuo pion mass value of
$m_{\pi}^{0}=130(2)$ MeV.

###  III.6 Comparison with results from the literature

Comparing our pion quasiparticle mass $\omega_{\mathbf{0}}(T)$ and
quasiparticle decay constant $f_{\pi}^{t}(T)$ at $T=128$ MeV to the matching
quantities at the corresponding zero temperature ensemble we get
$\omega_{\mathbf{0}}(T)/m_{\pi}^{0}=0.892(14)$ and
$f_{\pi}^{t}(T)/f_{\pi}^{0}=1.023(10)$. Thus, the quasiparticle mass decreases
at finite temperature while the quasiparticle decay constant experiences a
slight increase. This behavior is similar to what is found in a ChPT
calculation (see Ref. Toublan:1997rr , Fig. 3 and Fig. 4) at two loops, where
the reduction of the quasiparticle mass therein is $\approx 0.9$. On the other
hand the quasiparticle decay constant increases by a factor of approximately
1.03. That the pion pole mass decreases with temperature is also supported by
a recent lattice-improved soft-wall AdS/QCD model study Wen:2024hgu . On the
chiral crossover ensemble, we obtain
$\omega_{\mathbf{0}}(T)/m_{\pi}^{0}=0.63(2)$ and
$f_{\pi}^{t}(T)/f_{\pi}^{0}=0.69(3)$. In the high temperature phase, we obtain
$\omega_{\mathbf{0}}(T)/m_{\pi}^{0}=0.17(1)$ and
$f_{\pi}^{t}(T)/f_{\pi}^{0}=1.57(5)$.

Regarding the screening pion mass $m_{\pi}$, we found that it increases with
temperature compared to $m_{\pi}^{0}$. The ratio is
$m_{\pi}/m_{\pi}^{0}=1.105(14)$ at $T=128\,$MeV. This increase is also
supported qualitatively by the study of Son and Stephanov near the chiral
phase transition Son:2001ff . On the chiral crossover ensemble we obtain
$m_{\pi}/m_{\pi}^{0}=1.460(27)$, while the high temperature phase ensemble
gives $m_{\pi}/m_{\pi}^{0}=4.218(63)$.

In Ref. Goderidze:2022vlm , the authors work out the pion damping width and
the pion spectral function in the framework of a $SU(2)$ Nambu-Jona-Lasinio
(NJL) model for a few temperatures below the critical temperature
$T_{c}^{\text{NJL}}$ = 190 MeV. They observe that the position of the peak of
the pion spectral function at vanishing momentum $\mathbf{p}$ moves to the
right for increasing temperatures
$T/T_{c}^{\text{NJL}}\in\\{0,0.79,0.89,0.97\\}$ (see Fig. (3) in Ref.
Goderidze:2022vlm ). This contradicts our observation of the pion pole mass
being reduced at finite temperature.

##  IV Static mesonic screening masses

Figure 7:  Top panel: Screening vector and axial-vector effective masses at
$T=128\,$MeV. Middle panel: Screening vector and axial-vector effective masses
at $T=154\,$MeV. Bottom panel: Screening vector and axial-vector effective
masses at $T=192\,$MeV. The bands represent correlated fits to the plateau,
while the green dashed line represents the $2\pi T$ limit. Furthermore, the
value for the screening masses obtained from a direct fit of the axial-vector
respectively vector correlator are also included and represented by a single
light bar. The violet dotted line represents the free two pion states (for
$T=128\,$MeV and $T=154\,$MeV).

In this section we present our results for the screening masses in the
pseudoscalar, vector and axial-vector channels and compare them to those of
the HotQCD collaboration Bazavov:2019www . The latter have been obtained using
the highly improved staggered quark (HISQ) action and a continuum
extrapolation has been performed from a range of lattice spacings
corresponding to $N_{\tau}=6-16$. In contrast to Ref. Bazavov:2019www  we only
have results for mesons of the type $u\bar{d}$. A direct comparison of the
numerical values in the three channels is summarized in Table 5.

At $T=128\,$MeV, we obtain a pseudoscalar screening mass of $145(2)\,$MeV,
which is increased by $\approx 12\,\%$ compared to the pion mass on the
corresponding vaccuum ensemble [see Eq. (30)] and also compared to the HotQCD
result of $129(5)\,$MeV. The latter has been obtained at a slightly higher
temperature $T=132\,$MeV. The results of the HotQCD collaboration for the
vector and axial vector screening masses at the aforementioned temperature
read $0.7(2)\,$GeV and $1.0(2)\,$GeV, respectively and are compatible within
errors with our results of $614(19)\,$MeV and $832(35)\,$MeV, respectively.
The lightest vector and axial vector mesons with $u\bar{d}$ valence quarks are
the $\rho$ meson with a physical mass of $775\,$MeV, respectively the $a_{1}$
meson with a physical mass of $1.23\,$GeV. It is notable that both masses are
significantly larger $[m_{\rho}/m_{scr}^{V}\approx
1.26,m_{a_{1}}/m_{scr}^{A}\approx 1.48]$ than the screening masses extracted
on our $N_{\tau}=24$ ensemble at $T=128\,$MeV in the hadronic phase. However,
since the vector and axial-vector screening spectra start with a continuum,
not with an isolated pole, we also show the free two pion states with a violet
dotted line in Fig. 7.

Table 5: Comparison of the screening masses in the $u\bar{d}$ channel extracted at temperatures $T\in\\{128,154,192\\}\,$MeV to the continuum-extrapolated results obtained by the HotQCD collaboration [see Ref. Bazavov:2019www , Table X]. |  $m_{scr}^{PS}\,$[MeV] |  $m_{scr}^{V}\,$[MeV] |  $m_{scr}^{A}\,$[MeV]
---|---|---|---
Temperature [MeV] | HotQCD | this work | HotQCD | this work | HotQCD | this work
128 |  | 145(2) |  | 614(19) |  | 832(35)
132 | 129(5) |  | 700(200) |  | 1000(200) |
152 | 187(2) |  | 730(50) |  | 850(40) |
154 |  | 196(4) |  | 625(28) |  | 741(41)
156 | 202(3) |  | 750(60) |  | 830(60) |
192 | 540(10) | 555(8) | 1020(40) | 1069(5) | 1040(30) | 1068(12)

Next, we compare our results of the screening masses obtained on our
$N_{\tau}=20$ ensemble near the pseudocritical temperature to the continuum
extrapolated HotQCD results. As far as the pseudoscalar screening mass is
concerned our result $m_{scr}^{P}(T=154\,\text{MeV})=196(4)\,$MeV is in very
good agreement with the HotQCD results $187(2)\,$MeV and $202(3)\,$MeV at the
temperatures $T=152\,$MeV and $T=156\,$MeV, respectively. However, our vector
screening mass $m_{scr}^{V}(T=154\,\text{MeV})=625(28)\,$MeV differs by
$\approx 2\sigma$ compared to the HotQCD results 730(50) MeV and 750(60) MeV
at the temperatures $T=152\,$MeV and $T=156\,$MeV, respectively. Furthermore,
it is compatible within errors with our result obtained in the hadronic phase.
Our result for the axial-vector screening mass
$m_{scr}^{A}(T=154\,\text{MeV})=741(41)\,$MeV differs again by $\approx
2\sigma$ compared to the HotQCD result and is reduced by a factor of $\approx
1.12$ compared to our result at $T=128\,$MeV. Thus, the ratio of axial vector
screening mass and vector screening mass decreases from $\approx 1.36$ at
$T=128\,$MeV to $\approx 1.19$ at $T=154\,$MeV.

Finally, we confront our screening masses extracted from a $N_{\tau=16}$
ensemble in the chirally-symmetric phase with the results of the HotQCD
collaboration. The pseudoscalar screening mass
$m_{scr}^{P}(T=192\,\text{MeV})=555(8)\,$MeV is compatible within errors with
the HotQCD result $540(10)\,$MeV obtained at the same temperature. For the
vector and axial-vector screening masses we obtain $1069(5)\,$MeV and
$1068(12)\,$MeV, respectively. Both screening masses are compatible within
errors with the HotQCD collaboration results $1020(40)\,$MeV and
$1040(30)\,$MeV and their ratio is precisely one. Inspecting effective mass
plots displayed in the lower panel of Fig. 7, we notice that the degeneracy of
the vector and axial-vector screening mass at $T=192\,$MeV is not only evident
in the plateau region, but can also be seen (within errors) individually over
a long range of source-sink seperations. However, at this temperature, the
vector and axial vector screening masses are still $\approx 11\%$ below the
$2\pi T$ limit of screening masses. This high-temperature limit would
correspond to the propagation of two quark quasiparticles in the thermal
medium, each carrying a ‘mass’ of $\pi T$ corresponding to the first Matsubara
mode.

##  V Isovector quark number susceptibility

The quark number susceptibility (QNS) measures the response of the net quark
number density to an infinitesimal change in the quark chemical potential. On
the lattice, we determine the QNS via

$\displaystyle\delta^{ab}\chi_{q}(x_{0},T)=2Z_{V}^{2}(g_{0}^{2})\,\beta\,\int\mathrm{d}^{3}x\,\langle
V_{0}^{a}(x_{0},\mathbf{x})\,V_{0}^{b}(0,\mathbf{0})\rangle\,,\quad\quad
x_{0}\neq 0\,.$ (50)

For the QNS, no additive improvement of the vector current is needed.

Figure 8: Quark number susceptibility (QNS) extracted from the local vector
current correlator, Eq. (50). For comparison also $G_{A_{0}}(x_{0},T)$ and its
prediction in terms of screening quantities is shown. The top panel shows the
QNS for a temperature $T=128\,$MeV, the middle panel for a temperature
$T=154\,$MeV and the bottom panel for a temperature $T=192\,$MeV. The
correlators in these plots follow the normalization given in Eq. (50), with
the exception of $G_{A_{0}}(x_{0},T=128\,\text{MeV})$, as indicated by the
legend, for better visibility.

The results are shown in Fig. 8 in temperature units. Note that computing the
isospin susceptibility does not involve any contributions from disconnected
diagrams. In Ref. Borsanyi:2011sw  (see Table I) the QNS was determined as a
function of the temperature using 2+1 dynamical staggered quark flavors and,
additionally, a continuum extrapolation was done. Their result is
$\chi_{q}(T)/T^{2}=0.432(92)$ and $\chi_{q}(T)/T^{2}=0.481(87)$ for the
temperatures $T=125$ MeV and $T=130$ MeV, respectively (see Table I of Ref.
Borsanyi:2011sw ). Our result,

$\chi_{q}(T)/T^{2}=0.485(7)\,,\qquad(T=128\,{\rm MeV})$ (51)

is compatible with both of these results. At the crossover, our estimate of
the QNS

$\chi_{q}(T)/T^{2}=0.744(7)\,,\qquad(T=154\,{\rm MeV})$ (52)

overshoots theirs, $\chi_{q}(T)/T^{2}=0.669(47)$ at a temperature $T=155\,$MeV
by $\approx 1.6\sigma$. Also in the high temperature phase our estimate

$\chi_{q}(T)/T^{2}=0.856(3)\,,\qquad(T=192\,{\rm MeV})$ (53)

overshoots their result of $0.810(17)$ by $\approx 2.7\sigma$. Note however
that our results are not continuum-extrapolated. For comparison also the axial
charge correlator $G_{A_{0}}(x_{0},T)$ is shown in the plots. Already at
$T=128\,$MeV the correlator is quite flat, but a factor of $\approx 2$ bigger
in magnitude compared to the isospin charge correlator. Around the crossover
at $T=154\,$MeV the axial charge correlator is only $6\%$ bigger than the
isospin charge correlator. Finally, at $T=192\,$MeV the two correlators are
near-degenerate. Next, we are going to compare our lattice estimate for the
QNS with the hadron resonance gas (HRG) model and also test an alternative
‘quasiparticle gas’ model employing our modified dispersion relation for the
pion quasiparticle.

###  V.1 Comparison with the hadron resonance gas model

The HRG model Hagedorn:1984hz ; Karsch:2003vd  describes the thermodynamic
properties and the quark number susceptibilities of the low-temperature phase
rather well. It assumes that the thermodynamic properties of the system are
given by the sum of the partial contributions of non-interacting hadron
species, i.e.

$\displaystyle\text{ln}[Z(T,V)]=-\frac{V}{2\pi^{2}}\sum_{i}\int_{0}^{\infty}\,\mathrm{d}p\,p^{2}\,\text{ln}[1-\eta_{i}\,e^{-\sqrt{m_{i}^{2}+\mathbf{p}^{2}}/T}]\,,$
(54)

where $\eta_{i}=\pm 1$ takes into account bosons and fermions respectively.
The sum extends over all resonances up to a mass of 2.0 GeV, since for most of
them the width is not large compared to the temperature.

The QNS can be obtained as the sum Brandt:2015aqk

$\chi_{q}(T)=(\chi_{q})_{\text{mesons}}+(\chi_{q})_{\text{baryons}}\,,$ (55)

where

$\displaystyle\frac{(\chi_{q})_{\text{mesons}}}{T^{2}}$
$\displaystyle=\frac{2\beta^{3}}{3}\sum_{\text{multiplets}}(2J+1)I(I+1)(2I+1)\int\frac{\mathrm{d^{3}}\mathbf{p}}{(2\pi)^{3}}\,f_{\mathbf{p}}^{B}(1+f_{\mathbf{p}}^{B})\,,$
(56) $\displaystyle\frac{(\chi_{q})_{\text{baryons}}}{T^{2}}$
$\displaystyle=\frac{2\beta^{3}}{3}\sum_{\text{multiplets}}(2J+1)I(I+1)(2I+1)\int\frac{\mathrm{d^{3}}\mathbf{p}}{(2\pi)^{3}}\,f_{\mathbf{p}}^{F}(1-f_{\mathbf{p}}^{F})\,,$
(57)

and $f_{\mathbf{p}}^{B/F}=1/[e^{\beta\omega_{\mathbf{p}}}\mp 1]$ are the Bose-
Einstein and Fermi-Dirac distributions, respectively. The sums are carried out
over all multiplets of spin $J$ and isospin $I$ that are not identical.
Especially particles and antiparticles have to be considered separately. This
results in an additional factor of two in the baryon case and for mesons with
strange quark constituents.

Employing Eqs. (56-57) within the HRG model and summing all resonances up to a
mass of 2 GeV, we obtain $\chi_{q}(T)/T^{2}=0.486$ which agrees perfectly with
our lattice estimate, Eq. (51). The relative composition of the total QNS is
summarized in Fig. 9 It is questionable whether resonances like the
$K_{0}^{*}(700)$, whose full Breit-Wigner width $(478\pm 50\,\text{MeV})$ is
higher than the temperature considered should be taken into account.

An alternative to the HRG model is to only include the pion contribution,
however taking into account the modified dispersion relation (14) at low
momenta,

$\displaystyle\frac{\chi_{q}}{T^{2}}=4\beta^{3}\int_{\absolutevalue{\mathbf{p}}<\Lambda_{p}}\frac{\mathrm{d^{3}}\mathbf{p}}{(2\pi)^{3}}\,f_{\mathbf{p}}^{B}(\omega_{\mathbf{p}})(1+f_{\mathbf{p}}^{B}(\omega_{\mathbf{p}}))\,,$
(58)

At this point we only integrate up to the momentum cut off
$\Lambda_{p}=400\,\text{MeV}$ since it is not clear if the thermal width of
the pion is still negligible for $\absolutevalue{\mathbf{p}}>\Lambda_{p}$ and,
as a consequence, including contributions from higher momenta may not be
justified. In Ref. Brandt:2015sxa  the authors have performed fits to the
axial-charge density correlator at nonvanishing momenta. They then have
compared the residue obtained from the fit parameters to the chiral prediction

$\text{Res}(\omega_{\mathbf{p}})=f_{\pi}^{2}(m_{\pi}^{2}+\mathbf{p}^{2})(1+b(\mathbf{p}))\,.$
(59)

For $\absolutevalue{\mathbf{p}}\approx 400\,$MeV they obtained a small
$b=-0.08(3)$ indicating that the chiral effective theory description still
works. To estimate our error, we thus quote the results for
$\Lambda_{p}\in\\{400,500\\}\,\text{MeV}$. Making use of Eq. (58) with
$u_{m}=0.807$, a screening pion mass $m_{\pi}=145\,\text{MeV}$ and
$\Lambda_{p}=400\,$MeV we obtain $\chi_{q}(T)/T^{2}=0.402$ which is $\approx
17\%$ below the lattice estimate. On the other hand, for
$\Lambda_{p}=500\,$MeV we obtain $\chi_{q}(T)/T^{2}=0.512$ which is $\approx
5\%$ above the lattice estimate. Note that in this model the sum over the
resonances is absent. The contributions of the other hadrons are taken into
account indirectly via the modified dispersion relation, since the collisions
of the pions among themselves and with other hadrons give rise to the modified
pion dispersion relation.

Employing the HRG for the ensemble at the chiral crossover, we get
$\chi_{q}(T)/T^{2}=0.701$, a result that undershoots the lattice estimate by
$5.8\%$. The relative composition of the total QNS predicted by the HRG model
in this case is given in Fig. 10. At both temperatures the pion clearly
dominates the contribution to the total QNS.

77.9%8.1%7.3%3.8%2.9%pionvector- and pseudoscalar meson octet without
$\rho$$\rho$ vector mesonbaryon octet and decupletheavier meson and baryon
resonances up to $2.0\,\text{GeV}$ Figure 9: Relative composition of the total
quark number susceptibility predicted by the hadron resonance gas model in the
hadronic phase.  59.7%10.1%11.3%10.0%8.9%pionvector- and pseudoscalar meson
octet without $\rho$$\rho$ vector mesonbaryon octet and decupletheavier meson
and baryon resonances up to $2.0\,\text{GeV}$ Figure 10: Relative composition
of the total quark number susceptibility predicted by the hadron resonance gas
model at the chiral crossover.

##  VI Order parameters for chiral symmetry restoration

In this section, several order parameters for chiral symmetry restoration are
investigated. Based on the pion screening quantities $m_{\pi}^{2}$ and
$f_{\pi}^{2}$ presented in Sec. III, we first evaluate an ‘effective chiral
condensate’ based on the Gell-Mann–Oakes–Renner relation. Additionally, we
explore two Euclidean-time dependent thermal correlation functions that are
order parameters for chiral symmetry and compare them to their zero-
temperature counterparts. We begin with the $(PA_{0})$-correlator, which
contains the pion pole that we have studied in Sec. III. We then consider the
difference of the (isovector) vector and axial-vector correlators. In the QCD
vacuum, the corresponding spectral functions are measured experimentally in
$\tau$ decays Davier:2005xq . They become degenerate in the chirally restored
phase of QCD. Their temperature dependence in the chirally broken phase has
been studied extensively in the framework of hadronic models supplemented by
sum rules Rapp:1999ej ; Kapusta:1993hq ; Hohler:2013eba .

###  VI.1 The Gell-Mann–Oakes–Renner relation

Following Ref. Brandt:2014qqa , we introduce an ‘effective chiral condensate’
based on the Gell-Mann–Oakes–Renner (GOR) relation,

$C_{\rm GOR}^{\rm GF}\equiv-\frac{f_{\pi}^{2}m_{\pi}^{2}}{m_{\text{q}}}\,.$
(60)

In the chiral limit, $C_{\rm GOR}^{\rm GF}\to\langle\bar{\psi}\psi\rangle$.
Additionally, since above $T_{c}$, $m_{\pi}\sim T$ and $f_{\pi}\sim
m_{\text{q}}$, $C_{\rm GOR}^{\rm GF}$ is of
$\mathcal{O}(m_{\text{q}}\,T^{2})$. Thus, it serves as an order parameter for
chiral symmetry. Using $m_{\text{q}}=m_{\text{PCAC}}$ and the screening
quantities of Table 3 we obtain

$\displaystyle\left|C_{\rm GOR}^{\rm GF}\right|^{1/3}_{T=128\,\text{MeV}}$
$\displaystyle=291(3)\,\text{MeV}\,,$ $\displaystyle\left|C_{\rm GOR}^{\rm
GF}\right|^{1/3}_{T=154\,\text{MeV}}$ $\displaystyle=180(8)\,\text{MeV}\,,$
(61) $\displaystyle\left|C_{\rm GOR}^{\rm
GF}\right|^{1/3}_{T=192\,\text{MeV}}$ $\displaystyle=133(6)\,\text{MeV}\,.$

The value of the chiral condensate has been extracted in the gradient flow
scheme just like the PCAC mass (see Sec. II.4). Comparing with the chiral
condensate on the corresponding zero-temperature ensemble Ce:2022kxy , we get

$\displaystyle\left[\frac{C_{T\approx 128\,\text{MeV}}}{C_{T\approx
0\,\text{MeV}}}\right]_{\text{GOR}}$
$\displaystyle\equiv\frac{(f_{\pi}^{2}m_{\pi}^{2})_{T\approx
128\,\text{MeV}}}{(f_{\pi}^{2}m_{\pi}^{2})_{T\approx
0\,\text{MeV}}}=0.84(3)\,,$ $\displaystyle\left[\frac{C_{T\approx
154\,\text{MeV}}}{C_{T\approx 0\,\text{MeV}}}\right]_{\text{GOR}}$
$\displaystyle\equiv\frac{(f_{\pi}^{2}m_{\pi}^{2})_{T\approx
154\,\text{MeV}}}{(f_{\pi}^{2}m_{\pi}^{2})_{T\approx
0\,\text{MeV}}}=0.21(3)\,,$ (62) $\displaystyle\left[\frac{C_{T\approx
192\,\text{MeV}}}{C_{T\approx 0\,\text{MeV}}}\right]_{\text{GOR}}$
$\displaystyle\equiv\frac{(f_{\pi}^{2}m_{\pi}^{2})_{T\approx
192\,\text{MeV}}}{(f_{\pi}^{2}m_{\pi}^{2})_{T\approx
0\,\text{MeV}}}=0.08(1)\,.$

which corresponds to a reduction by 16% in the hadronic phase. Notably, during
the temperature increase from $T=128\,$MeV to $T=154\,$MeV the effective
chiral condensate (normalized to the zero-temperature equivalent) decreases by
a factor of four down to $0.21(3)$. Finally, at a temperature of $T=192\,$MeV
the effective chiral condensate reaches a value of $\approx 8\%$ of its vacuum
value. The reduction in the low-temperature phase is compatible within the
scope of the error with a three-loop result of Gerber and Leutwyler (see Ref.
Gerber:1988tt , Fig. 5).

###  VI.2 The $(PA_{0})$-correlator

In Ref. Brandt:2014qqa  it was shown that the temporal $(PA_{0}$)-correlator
can be predicted exactly in the chiral limit,

$G_{PA_{0}}(x_{0},T)=\frac{\langle\bar{\psi}\psi\rangle}{2\beta}\left(x_{0}-\frac{\beta}{2}\right)\,.$
(63)

As can be seen from Eq. (63) the $(PA_{0})$-correlator is antisymmetric around
$\beta/2$. Consequently, we set the point $x_{0}=\beta/2$ to zero. In our
analysis we have averaged over the $(PA_{0})$\- and $(A_{0}P)$-correlator. The
latter can be obtained interchanging source and sink.444This is only true for
our vacuum correlators, which are obtained from point-sources. For our finite
temperature correlators obtained from stochastic wall-sources these are
equivalent, as we use a symmetric source-sink setup. We find that the
estimator based on the last expression of Eq. (10) has the better signal-to-
noise ratio. Since this correlator is proportional to the chiral condensate
$\langle\bar{\psi}\psi\rangle$, it can serve as an order parameter for chiral
symmetry restoration as well. Looking at the ratio of the thermal over the
reconstructed correlator obtained via Eq. (67), we find the reductions of the
$(PA_{0})$ correlator to be in satisfactory agreement with the reductions we
have obtained for the chiral condensate using the Gell-Mann–Oakes–Renner
relation (see Eq. (62)) for all three temperatures, namely

$\lim_{x_{0}\to\beta/2}\frac{G_{PA_{0}}(x_{0},T)}{G_{PA_{0}}^{\rm
rec}(x_{0},T)}=\left\\{\begin{array}[]{ll}0.817(21)\,,&N_{\tau}=24\\\
0.276(23)\,,&N_{\tau}=20\\\ 0.088(5)\,,&N_{\tau}=16.\end{array}\right.$ (64)

The corresponding plots can be found in Fig. 18 in Appendix G.

Figure 11:  Top panel: Ratio of the difference ‘$V-A$’ and the difference of
the reconstructed correlator ‘$(V-A)_{\text{rec}}$’ at $T=128\,$MeV. Middle
panel: Same quantity at $T=154\,$MeV. Bottom panel: Same quantity at
$T=192\,$MeV. The green band shows the result of a correlated fit to the
plateau.

Figure 12:  The Backus-Gilbert spectral function
$\hat{\rho}_{V-A}^{\text{rec}}(\omega)/(f\,T^{2})$ obtained from the
reconstructed $V-A$ correlator is displayed as a blue band, at a temperature
of $T=128\,$MeV and $154\,$MeV respectively for the left and right panel. The
corresponding estimator of the thermal spectral function
$\hat{\rho}_{V-A}^{\text{ther}}(\omega)/(f\,T^{2})$ is shown as a red band.
For orientation, the blue vertical line corresponds to the physical mass of
the $\rho$ meson, while the vertical red line corresponds to the vector
screening mass at the respective temperature [see Table 5]. The underlying
temporal vector correlators are $\mathcal{O}(a)$ improved.

###  VI.3 Dey-Eletsky-Ioffe mixing theorem at finite quark mass

In Refs. Dey:1990ba ; Eletsky:1992ay ; Eletsky:1994rp  it was demonstrated,
using PCAC current algebra, that near the chiral limit the finite-temperature
vector and axial-vector correlators can be described with the help of their
vacuum counterparts. In terms of the corresponding spectral functions, this
statement reads

$\displaystyle\rho_{V}(\omega,\mathbf{p},T)-\rho_{A}(\omega,\mathbf{p},T)=(1-2\epsilon(T))\left[\rho_{V}(\omega,\mathbf{p},T=0)-\rho_{A}(\omega,\mathbf{p},T=0)\right]\,,$
(65)

while the sum $(\rho_{V}+\rho_{A})$ receives no thermal correction at this
order. The coefficient

$\epsilon(T)\equiv\frac{T^{2}}{6(f_{\pi}^{0})^{2}}$ (66)

is a temperature dependent expansion parameter in powers of the pion density.
Thus, the difference in Eq. (65) serves as an order parameter for chiral
symmetry restoration.

Inspired by the Dey-Eletsky-Ioffe mixing theorem, in the following we consider
the difference $(G_{V}(x_{0},\mathbf{p}=0,T)-G_{A}(x_{0},\mathbf{p}=0,T))$ at
vanishing spatial momentum. Using the corresponding quasi zero-temperature
E250 ensemble of size $192\times 96^{3}$, we calculate the ‘reconstructed’
correlator $G_{V}^{\text{rec}}-G_{A}^{\text{rec}}$ for the difference. This
correlator can be interpreted as the thermal Euclidean correlator that would
be realized if the spectral function was unaffected by thermal effects, during
the transition from $T\approx 0$ to $T\in\\{128,154,192\\}\,$MeV. Following a
method first proposed in Ref. Meyer:2010 , we obtain the reconstructed
correlator as

$\displaystyle
G_{J}^{\text{rec}}(x_{0},T,\mathbf{p})=\sum_{m\in\mathbb{Z}}\,G_{J}(\absolutevalue{x_{0}+m\beta},0,\mathbf{p})\
\ \ \ \ \ \ (J=\\{V,A,\dots\\}).$ (67)

It is based on the identity of the kernel function

$\displaystyle\frac{\text{cosh}(\omega(\beta/2-x_{0}))}{\text{sinh}(\omega\beta/2)}=\sum_{m\in\mathbb{Z}}\,e^{-\omega\absolutevalue{x_{0}+m\beta}}\,.$
(68)

The difference ‘$V-A$’ for the thermal ensembles are shown in the l.h.s. of
Fig. 20 in Appendix G, while the r.h.s. displays the corresponding
reconstructed correlators obtained from the $T\approx 0$ ensemble using Eq.
(67). In both cases renormalization factors are included following Eqs.
(81)-(82). Furthermore, the vector-current has been
$\mathcal{O}(a)$-improved555The axial-vector current is automatically
$\mathcal{O}(a)-$improved since the spatially integrated improvement term
vanishes in this case. On the other hand, for the vector current the term
$\propto T_{i0}$ contributes in the improvement process [see Eqs. (78)-(79)].
using the improved derivative defined in Eq. (A). The top two panels show the
relevant correlator difference ‘$V-A$’ at $T=128\,$MeV, the middle panels
‘$V-A$’ at $T=154\,$MeV and the bottom panels ‘$V-A$’ at $T=192\,$MeV. It is
notable that we are completely dominated by the error on the corresponding
vacuum ensemble. The ratios of the thermal and reconstructed correlators are
shown in Fig. 11. It is immediately apparent that the ratio differs
significantly from one at all three temperatures, even though one expects the
ratio to tend to unity for $x_{0}\to 0$ for $m_{q}\neq 0$. Remembering the
physical interpretation of the reconstructed correlator given above Eq. (67),
we can conclude that also the corresponding spectral function must have
changed during the transition from $T\approx 0$ to
$T\in\\{128,154,192\\}\,$MeV.

As Eq. (65) is obtained in the chiral limit, one would not necessarily expect
to see a drop to a constant value in the ratio even for finite quark masses.
However, as can be seen in Fig. 11, we observe a flat behavior around the
midpoint at $T=128\,$MeV and $T=154\,$MeV. Especially in the $N_{\tau}=20$
ensemble the ratio is almost constant in the range $6\leq x_{0}/a\leq 10$. In
the hadronic phase, we read off a reduction to a factor $0.69(2)$ pointing to
chiral symmetry restoration being already at an advanced stage in the hadronic
phase [see top panel of Fig. (11)]. During the transition from $T=128\,$MeV to
$T=154\,$MeV the ratio in the ‘$V-A$’-channel drops drastically by a factor of
$\approx 9$ down to a reduction factor of $\approx 0.08(1)$. This is confirmed
by inspecting the smoothed $(V-A)$ spectral function obtained from the Backus-
Gilbert method in Fig. 12. Increasing the temperature further to $T=192\,$MeV
well in the chirally symmetric phase, the ratio even becomes slightly negative
($\approx-0.02$). Due to the short extent in the time direction, cut-off
effects could be responsible for this behavior.

Summarizing our results, we find that at temperatures
$T\in\\{128,154,192\\}\,$MeV the ’effective chiral condensate’, which is
extracted from static screening quantities, is reduced to
{$0.84(3),0.21(3),0.08(1)\\}$ times the value at the corresponding vacuum
ensemble. These values are compatible within errors with the reduction to
$\\{0.817(21),0.276(23),0.088(5)\\}$ times the value at the corresponding
vacuum ensemble observed in the temporal $(PA_{0})$-correlator, which, in the
chiral limit, is proportional to the chiral condensate. Especially in the low-
and high-temperature phase these two predictions are in good agreement. This
reduction factors can be compared to the reduction in the temporal $(V-A)$
sector. While the reduction to $0.69(2)$ times the vacuum value at
$T=128\,$MeV is in the same ballpark, we find that during the transition to
the crossover region chiral symmetry restoration appears to happen more
rapidly. At $T=154\,$MeV the reduction factor in the $V-A$ sector is smaller
by a factor of $\approx 4$. Increasing the temperature further to $T=192\,$MeV
the ratio in the $PA_{0}$-channel drops down to $\approx 0.09$, while the
ratio in the temporal $(V-A)$-channel is compatible with zero.

The different rate at which different order parameters of chiral symmetry
breaking approach zero already shows up in predictions from chiral
perturbation theory at vanishing quark masses. Specifically the suppression of
the chiral condensate to leading order in $T^{2}$ goes like
$1-\frac{T^{2}}{8(f_{\pi}^{0})^{2}}$ (see Ref. Gasser:1986vb ; Gerber:1988tt )
and is therefore less pronounced relative to Eq. (65) with
$1-2\epsilon=1-\frac{T^{2}}{3(f_{\pi}^{0})^{2}}$ Dey:1990ba .

##  VII Vacuum subtracted vector and axial-vector spectral functions

By inspection of Eq. (13) it becomes evident that the Euclidean correlators
depend on the temperature in a dual way: implicitly through the occurrence of
the inverse temperature $\beta$ in the kernel function
$\text{cosh}(\omega(\beta/2-x_{0}))/\text{sinh}(\omega\beta/2)$ as well as
through the explicit temperature dependence of the spectral function. The
first effect leads to a temperature dependence of the Euclidean correlator
even if the spectral function does not change. It is important to isolate this
effect from a genuine temperature dependence resulting from a change in the
thermal medium (e.g. the quark-gluon plasma). Consequently, differences in the
thermal and reconstructed correlator can be traced back to a change in the
spectral function, and thus the physical situation. The inverse is not
necessarily true Aarts:2020dda .

Figure 13:  Top panels: The l.h.s. shows the smeared spectral function of the
difference between the finite temperature and reconstructed zero-temperature
temporal vector correlators projected to zero momentum at a temperature
$T=128\,$MeV. The r.h.s shows the same quantity for the temporal axial-vector
correlators projected to zero momentum. Middle panels: The same quantities at
a temperature $T=154$ MeV. Bottom panels: Again the same quantities at a
temperature $T=192$ MeV.

Using the Backus-Gilbert method described in Sec. III.5, we provide an
estimator for the smeared spectral function of the difference between the
finite temperature and reconstructed zero-temperature temporal vector
correlators, projected to zero momentum. For the analysis of vacuum subtracted
spectral functions we use

$f_{2}(\omega):=\omega$ (69)

as a rescaling function. The results of the smeared and rescaled spectral
functions are shown in Fig. 13.

The difference

$\Delta\rho_{V}(\omega,\bm{p},T)\equiv\rho_{V}(\omega,\bm{p},T)-\rho_{V}(\omega,\bm{p},T=0)\,,$
(70)

between the spatial components of the finite and zero-temperature spectral
function satisfies the sum rule Bernecker:2011gh

$\int_{0}^{\infty}\frac{\mathrm{d}\omega}{\omega}\,\Delta\rho_{V}(\omega,\bm{0},T)=0$
(71)

in the thermodynamic limit. From the operator-product expansion (OPE) it is
known that the difference in Eq. (70) falls of like $1/\omega^{2}$ for large
frequencies $\omega$ Caron-Huot:2009ypo . Applying the sum rule to our smeared
spectral function obtained from the Backus-Gilbert method, Eq. 71 evaluates to
$-0.30(29)$ at $T=128\,$MeV, $-0.46(37)$ at $T=154\,$MeV and $0.44(1)$ at
$T=192\,$MeV. Thus, the sum rule is satisfied within the large errors at the
lower two temperatures. In the high temperature phase the sum rule is not
satisfied. However, it should be stressed that due to the short time extent of
the E250Nt16 ensemble, we could only use $4\leq x_{0}/a\leq 8$ for the
spectral function reconstruction. Therefore, this result should be interpreted
with caution.

For the thermal modification of the axial-vector channel, we have the spectral
representation

$2\int_{0}^{\infty}\frac{\mathrm{d}\omega}{\omega}\,\Delta\rho_{A}(\omega,\bm{p},T)=\int_{0}^{\beta}\mathrm{d}x_{0}\,G_{A}(x_{0},\bm{p},T)-\int_{-\infty}^{\infty}\mathrm{d}x_{0}\,G_{A}(x_{0},\bm{p},0).$
(72)

Close to the chiral limit, we can saturate the correlation functions on the
r.h.s. at the zero four-momentum by the screening pion, leading to the sum
rule

$2\int_{0}^{\infty}\frac{\mathrm{d}\omega}{\omega}\,\Delta\rho_{A}(\omega,\bm{0},T)=f_{\pi}^{2}(T=0)-f_{\pi}^{2}(T)+{\mathcal{O}}(m_{q}^{2}).$
(73)

Obviously, this is similar to the first Weinberg sum rule Weinberg:1967kj  and
its finite-temperature generalisation Kapusta:1993hq , however, the integral
on the l.h.s. of Eq. (73) is convergent even at $m_{q}\neq 0$. We have seen
that $f_{\pi}^{2}(T)$ is already significantly reduced at $T=128$ MeV, and it
is of ${\mathcal{O}}(m_{q}^{2})$ in the chirally restored phase. Thus the
r.h.s of Eq. (73) soon becomes only weakly temperature-dependent. This
qualitative behaviour is confirmed by the numerical data shown in the right
panels of Fig. 13. Since $\Delta\rho(\omega,\bm{p},T)\sim 1/\omega^{2}$ at
large $\omega$, the sum rule suggests that thermal effects lead to the build-
up of spectral weight below the mass of the vacuum $a_{1}(1260)$ meson.

##  VIII Conclusion

Using lattice simulations at physical $(u,d,s)$ quark masses, we have shown
that the real part of the low-energy pole in the two-point function of the
axial charge is reduced as the temperature increases throughout the hadronic
phase of QCD. We interpret this as a reduction of the pion quasiparticle mass
in the thermal medium. At our lowest temperature of 128 MeV, the effect
amounts to a ten percent reduction. This finding is possible because the pole
describing the pion quasiparticle parametrically dominates both the Euclidean
two-point functions of the axial charge and of the pseudoscalar density at
$x_{0}=\beta/2$.

Above the crossover temperature, the pole in the axial-charge correlator is
expected to be of diffusive nature, i.e. its imaginary part dominates over the
real part. Since the curvature in the Euclidean correlator of the axial charge
is not parametrically dominated by the pole in the chirally restored phase,
one faces an inverse problem if one wants to determine the precise position of
the pole. What is certain is that the pole tends to the origin, $\omega=0$, in
the chiral limit. Indeed our calculation shows that the axial-charge
correlator is near-degenerate with the (temporally constant) isospin charge
correlator.

We have found that the static pseudoscalar correlation length, $m_{\pi}^{-1}$,
shrinks steadily with increasing temperature, consistently with previous
lattice calculations. The associated ‘decay constant’ $f_{\pi}$, which we
defined (in Eq. (5)) as the coupling of the spatial, longitudinal axial
current to this correlation length, decreases steadily with temperature.
Parametrically, it is of order the quark mass in the chirally restored phase.
Similarly, we have computed the spatially transverse vector and axial-vector
static screening lengths, finding mostly consistency with results by the
HotQCD collaboration.

The combination $f_{\pi}^{2}m_{\pi}^{2}/m_{q}$ is an estimator for the
condensate $\langle-\bar{\psi}\psi\rangle$ in the chiral limit, and we have
computed the ratio of $f_{\pi}^{2}m_{\pi}^{2}$ to its value in the vacuum as
an estimator for the ‘melting’ of the chiral condensate. That ratio was
compared to a different estimator based on the temporal correlation function
of the axial charge with the pseudoscalar density, yielding consistent
results.

We have constrained the thermal modification of the vector and the axial-
vector spectral functions by constructing the difference of the thermal and
the ‘reconstructed’ correlator, i.e. the correlator that would be realised at
finite temperature if no change took place in the spectral function. This
strategy has been used before in a two-flavour calculation with $m_{\pi}\simeq
270\,$MeV Brandt:2015aqk , but the present calculation is the first one to
operate at physical quark masses around the crossover temperature. Using the
Backus-Gilbert method to directly analyze the thermal change in the spectral
function, we conclude that in the vector channel, excess spectral weight
appears at low frequencies, while a depletion occurs around energies of one
GeV. Such a change is consistent with the sum rule (71). In the axial-vector
case, we observe a significant excess thermal spectral weight up to one GeV –
see Fig. 13.

Finally, we have investigated the difference of the thermal correlators for
spatial components of the vector and axial-vector currents (at zero spatial
momentum). Phenomenologically, it is often studied in the context of the
dilepton production by the thermal medium Hohler:2013eba . This difference,
too, is an order parameter for chiral symmetry, and it provides additional
information on its restoration. This correlator is already suppressed relative
to its vacuum analogue by about at factor of two thirds at our lowest
temperature, $T\simeq 128$ MeV, and by more than an order of magnitude at the
crossover, $T\simeq 154$ MeV. The flat behaviour of this ratio is consistent
with an overall suppression factor of the $V-A$ spectral function up to fairly
high energies. This qualitative behaviour was the one originally predicted at
low temperatures in the chiral limit Dey:1990ba , and we have found it to be a
good approximation at the correlator level up to the crossover.

Overall, we have found simulations at physical quark masses with
$\mathcal{O}(a)$-improved Wilson fermions to be entirely feasible with the
current algorithmic techniques. We have observed long auto-correlation times
for the two-point correlation function of the axial current with the
pseudoscalar density at the chiral crossover, which is not surprising for an
order parameter in the vicinity of a phase transition. Future simulations
could be performed at finer lattice spacing in order to perform continuum
extrapolations of the most important observables such as the pion
quasiparticle mass.

###### Acknowledgements.

We thank Tim Harris for many helpful discussions related to this project. This
work was supported by the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation program through Grant Agreement
No. 771971-SIMDAMA, as well as by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) through the Cluster of Excellence “Precision
Physics, Fundamental Interactions and Structure of Matter” (PRISMA+ EXC
2118/1) funded by the DFG within the German Excellence strategy (Project ID
39083149). The research of M.C. is funded through the MUR program for young
researchers “Rita Levi Montalcini”. The authors gratefully acknowledge the
Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this
project by providing computing time on the GCS Supercomputer HAWK at
Höchstleistungsrechenzentrum Stuttgart (www.hlrs.de), as well as through the
John von Neumann Institute for Computing (NIC) on the GCS Supercomputer JUWELS
at Jülich Supercomputing Centre (JSC) under project SIMCHIC. The measurements
of the two-point functions have been partly performed on our local machine
MogonII and partly on Noctua2 at PC2 in Paderborn as part of an NHR proposal.
A.K. wishes to thank Simon Kuberski for providing the pion decay constant on
the E250 vacuum ensemble, as well as for sharing high-precision vector
correlator data on the aforementioned box Kuberski:2023qgx . He also
acknowledges Konstantin Ottnad for sharing his data in the axial-pseudoscalar
channel on the E250 vacuum ensemble. We made use of the following libraries:
GNU Scientific Library Galassi:2019 , GNU MPFR 10.1145/1236463.1236468  and
GNU MP Granlund:2020 .

##  Appendix A Improvement process

Lattice derivatives lead to cutoff effects. Therefore, in this section we
specify the lattice discretized (improved) derivatives used in this work, as
well as the $\mathcal{O}(a)-$improvement process of the axial-vector and
vector current. Defining a forward and backward derivative as

$\displaystyle\partial_{\mu}f(x)$
$\displaystyle=\frac{1}{a}\left(f(x+a\hat{\mu})-f(x)\right)+\mathcal{O}(a)\,,$
(74) $\displaystyle\partial^{\,*}_{\mu}f(x)$
$\displaystyle=\frac{1}{a}\left(f(x)-f(x-a\hat{\mu})\right)+\mathcal{O}(a)\,,$
(75)

a symmetrized version with reduced discretization errors can be obtained by
combining forward and backward derivative

$\displaystyle\widetilde{\partial}_{\mu}f(x)=\frac{1}{2}\left(\partial_{\mu}f(x)+\partial^{\,*}_{\mu}f(x)\right)=\frac{1}{2a}\left(f(x+a\hat{\mu})-f(x-a\hat{\mu})\right)+\mathcal{O}(a^{2})\,.$
(76)

Note the absence of $\mathcal{O}(a)$ dicretization errors. The above
symmetrized derivative can be further improved by making the replacement
Guagnelli:2000jw

$\displaystyle\widetilde{\partial}_{\mu}f(x)$
$\displaystyle\rightarrow\widebar{\partial}_{\mu}f(x):=\widetilde{\partial}_{\mu}f(x)\left(1-\frac{1}{6}a^{2}\partial^{\,*}_{\mu}\partial_{\mu}f(x)\right)\,,$
$\displaystyle=\frac{1}{12a^{2}}\left[f(x-2a\hat{\mu})-8f(x-a\hat{\mu})+8f(x+a\hat{\mu})-f(x+2a\hat{\mu})\right]\,.$
(77)

When acting on smooth functions this improved lattice derivative has
$\mathcal{O}(a^{4})$ discretization errors.

Next, we define the $\mathcal{O}(a)-$improved axial-vector and vector current
as

$\displaystyle A_{\mu}^{\text{imp},b}(x)$
$\displaystyle=A_{\mu}^{b}(x)+ac_{A}(g_{0}^{2})\widetilde{\partial}_{\mu}P^{b}(x)\,,$
(78) $\displaystyle V^{\text{imp},b}_{\mu}(x)$
$\displaystyle=V_{\mu}^{b}(x)+ac_{V}(g_{0}^{2})\widebar{\partial}_{\nu}\,T_{\mu\nu}(x)\,,$
(79)

where

$T_{\mu\nu}^{a}(x)\equiv-\frac{1}{2}\bar{\psi}[\gamma_{\mu},\gamma_{\nu}]\frac{\tau^{a}}{2}\psi$
(80)

is the tensor current.

The non-perturbatively calculated coefficient $c_{A}$ was taken from Ref.
Bulava:2015bxa , and the coefficient $c_{V}$ from Ref. Heitger:2020mkp . For
the analysis of the $V-A$ correlators in Sec. VI.3, we use an improved version
of the symmetric derivative defined in Eq. A.

##  Appendix B Renormalization process

Following Korcyl:2016ugy , we renormalize the correlators in the following
way:

$\displaystyle G_{V}^{\text{ren.}}$
$\displaystyle=Z_{V}^{2}(g_{0}^{2})\left[1+2am_{l}b_{V}(g_{0}^{2})+6am_{\text{av}}\widetilde{b}_{V}(g_{0}^{2})+\mathcal{O}(a^{2})\right]G_{V}\,,$
(81) $\displaystyle G_{A}^{\text{ren.}}$
$\displaystyle=Z_{A}^{2}(g_{0}^{2})\left[1+2am_{l}b_{A}(g_{0}^{2})+6am_{\text{av}}\widetilde{b}_{A}(g_{0}^{2})+\mathcal{O}(a^{2})\right]G_{A}\,,$
(82) $\displaystyle G_{A_{0}}^{\text{ren.}}$
$\displaystyle=Z_{A}^{2}(g_{0}^{2})\left[1+2am_{l}b_{A}(g_{0}^{2})+6am_{\text{av}}\widetilde{b}_{A}(g_{0}^{2})+\mathcal{O}(a^{2})\right]G_{A_{0}}\,,$
(83) $\displaystyle G_{PA_{0}}^{\text{ren.}}$
$\displaystyle=Z_{P}(g_{0}^{2})Z_{A}(g_{0}^{2})\left[1+am_{l}b_{A}(g_{0}^{2})+3am_{\text{av}}\widetilde{b}_{A}(g_{0}^{2})+\mathcal{O}(a^{2})\right]G_{PA_{0}}\,,$
(84) $\displaystyle G_{P}^{\text{ren.}}$
$\displaystyle=Z_{P}^{2}(g_{0}^{2})G_{P}\,,$ (85)

with $g_{0}^{2}=6/\beta$ being the bare gauge coupling and

$\displaystyle am_{q}$
$\displaystyle=\frac{1}{2}\left(\frac{1}{\kappa_{\text{q}}}-\frac{1}{\kappa_{\text{cr.}}}\right)\,,\quad
q\in\\{l,s\\}\,,$ (86) $\displaystyle m_{\text{av}}$
$\displaystyle=\frac{1}{3}\left(2m_{l}+m_{s}\right)\,,$ (87)

denoting the bare subtracted quark mass and the ‘average quark mass’,
respectively. The values for the renormalization constants $Z_{J}$ and the
finite quark mass parameters $b_{J}$ are given in Table 6.

Table 6: Summary of the renormalization parameters. $Z_{V}(g_{0}^{2})$ |  Heitger:2020mkp  | $0.735\,71(10)$
---|---|---
$Z_{A}(g_{0}^{2})$ |  DallaBrida:2018tpn  | $0.769\,00(42)$
$Z_{P}(g_{0}^{2})$ |  Campos:2018ahf  | $0.347\,68$
$b_{V}(g_{0}^{2})$ |  Bali:2023sdi  | $1.388(17)$
$b_{A}(g_{0}^{2})$ |  Bali:2023sdi  | $1.232(15)$
$\widetilde{b}_{V}(g_{0}^{2})$ |  Bali:2023sdi  | $0.069(42)$
$\widetilde{b}_{A}(g_{0}^{2})$ |  Bali:2023sdi  | $-0.03(13)$
$\kappa_{\text{cr.}}$ |  Gerardin:2018kpy  | $0.137\,172\,6(13)$
$\kappa_{\text{l}}$ |  Ce:2022eix  | $0.137\,232\,867$
$\kappa_{\text{s}}$ |  Ce:2022eix  | $0.136\,536\,633$

##  Appendix C Error analysis

Figure 14:  Integrated autocorrelation time $(S=6)$ of the PCAC mass at fixed
source-sink seperation $x_{3}/a=20$ on the E250Nt20 at the chiral crossover.

For the estimation of the statistical errors of our (derived) observables we
use the $\Gamma$ method Madras:1988ei ; Wolff:2003sm ; Ramos:2018vgu  in the
implementation of the pyerrors package introduced in Ref. Joswig:2022qfe . If
$\tau_{\text{exp}}$ (exponential autocorrelation time) denotes the largest
autocorrelation time in the system, the autocorrelation function typically
exhibits an asymptotic exponential behavior,

$\frac{\Gamma_{F}(t)}{\Gamma_{F}(0)}\sim\text{exp}\left(-\frac{t}{\tau_{\text{exp}}}\right)\,.$
(88)

For any (derived) observable $F$, one defines

$\tau_{\text{int},F}(W)=\frac{1}{2}+\sum_{t=1}^{W}\frac{\Gamma_{F}(t)}{\Gamma_{F}(0)}\,,$
(89)

as the integrated autocorrelation time. In practice one needs to choose an
appropriate window $W$. In doing so, one is confronted with an trade-off
between capturing all the relevant information in the chosen window and
simultaneously not being dominated by noise contributions to
$\tau_{\text{int}}$ for large values of $t$. In Ref. Wolff:2003sm  an
automatic windowing procedure was introduced to determine $W$ in actual data
analysis. Based on the hypothesis
$\tau_{\text{exp}}=S\,\tau_{\text{int},F}(W)$ one chooses the factor $S$ such
that the sum of systematic and statistical error is minimized. The default
value is $S=2$. However, for the correct error estimation of the PCAC mass in
the E250Nt20 ensemble, we choose $S=6$ [see Fig. 14]. The final error based on
the integrated autocorrelation time reads

$\sigma_{F}=\sqrt{2\tau_{\text{int},F}\frac{\Gamma_{F}(0)}{N}}\,.$ (90)

We utilize the robust median and median absolute deviation (MAD) described to
identify two (out of 1000) outliers on the E250Nt20 ensemble Agadjanov:2023efe
; Krasniqi:2024inr .

##  Appendix D Chiral effective theory Lagrangian of Son and Stephanov

In the chiral effective theory approach of Son and Stephanov Son:2001ff
Son:2002ci  the dynamics of the pions at finite temperature is described by
the Lagrangian,

$\displaystyle\mathcal{L}_{\text{eff}}=\frac{f_{t}^{2}}{4}\langle\nabla_{0}\Sigma\nabla_{0}\Sigma^{\dagger}\rangle-\frac{f_{\pi}^{2}}{4}\langle\partial_{i}\Sigma\partial_{i}\Sigma^{\dagger}\rangle+\frac{m_{\pi}^{2}f_{\pi}^{2}}{2}\text{Re}\langle\Sigma\rangle\,,$
(91)

where $\Sigma$ denotes an $SU(2)$ matrix whose phase describes the pions,
$\nabla_{0}\Sigma=\partial_{0}\Sigma-\frac{i}{2}\mu_{I5}(\tau_{3}\Sigma+\Sigma\tau_{3})$
is the covariant derivative, $\mu_{I5}$ denotes the axial isospin chemical
potential and the trace is taken in flavor space. Note that in the presence of
a thermal medium Lorentz invariance is broken resulting in two independent
decay constants which are related through the pion velocity Son:2001ff ,

$\displaystyle u=\frac{f_{\pi}}{f_{t}}\,.$ (92)

##  Appendix E Determination of the static screening properties in the
pseudoscalar channel

Figure 15:  Effective mass plots for the cosh mass $m_{\text{cosh}}(x_{3})/T$
as a function of the $x_{3}$-coordinate in temperature units, obtained from
the pseudoscalar screening correlation function at zero spatial momentum
$G_{P}^{s}(x_{3},T)$. For comparison the value for the screening pion mass
$m_{\pi}/T$, obtained from a simultaneous fit of the pseudoscalar and axial-
pseudoscalar correlator is also included (green bars). Top panel: Hadronic
E250Nt24 ensemble. It is assumed that the effective mass plateau starts at
$x_{3}/a=25$. Middle panel: E250Nt20 ensemble at the chiral crossover. It is
assumed that the effective mass plateau starts at $x_{3}/a=18$. Bottom panel:
E250Nt16 ensemble in the high-temperature phase. It is assumed that the
effective mass plateau starts at $x_{3}/a=18$. For all three ensembles he
result of the fit to the effective mass values is represented by a $1-\sigma$
band.

In this appendix, we describe how the screening pion mass $m_{\pi}$ and the
screening decay constant $f_{\pi}$ can be calculated. In order to accomplish
this, we perform a simultaneous fit of the static screening pseudoscalar and
the axial-pseudoscalar correlator [see Eqs. (6)-(7)].

The spectral decomposition of the aforementioned correlators reads

$\displaystyle G_{P}^{s}(x_{3},T)$
$\displaystyle=\frac{\bra{0}\bar{u}\gamma_{5}d\ket{\pi}\bra{\pi}\bar{d}\gamma_{5}u\ket{0}}{2\,m_{\pi}}\left(e^{-m_{\pi}x_{3}}+e^{-m_{\pi}\cdot(L-x_{3})}\right)\,,$
(93) $\displaystyle G_{AP}^{s}(x_{3},T)$
$\displaystyle=\frac{\bra{0}\bar{u}\gamma_{3}\gamma_{5}d\ket{\pi}\bra{\pi}\bar{d}\gamma_{5}u\ket{0}}{2\,m_{\pi}}\left(e^{-m_{\pi}x_{3}}-e^{-m_{\pi}\cdot(L-x_{3})}\right)\,,$
(94)

which suggests the following one-state fit ansatz:

$\displaystyle G_{P}^{\text{fit}}(x_{3},T)$
$\displaystyle=\absolutevalue{A}^{2}\left(e^{-Cx_{3}}+e^{-C\cdot(L-x_{3})}\right)\,,$
(95) $\displaystyle G_{AP}^{\text{fit}}(x_{3},T)$
$\displaystyle=\absolutevalue{A}\cdot\absolutevalue{B}\left(e^{-Cx_{3}}-e^{-C\cdot(L-x_{3})}\right)\,,$
(96)

$A$, $B$ and $C$ being free fit parameters. The pion screening mass $m_{\pi}$
and the screening pion decay constant $f_{\pi}$ are obtained from the fit
parameters via

$m_{\pi}\overset{\circ}{=}C,\qquad
f_{\pi}\overset{\circ}{=}\sqrt{\frac{2\absolutevalue{B}^{2}}{C}}\;.$ (97)

The ‘cosh mass’ is defined as the solution of the algebraic equation,

$\displaystyle\frac{G_{P}^{s}(x_{3}+a,T)}{G_{P}^{s}(x_{3},T)}=\frac{\text{cosh}[m_{\text{cosh}}(x_{3}+a-L/2)]}{\text{cosh}[m_{\text{cosh}}(x_{3}-L/2)]}\,.$
(98)

It is visualized in Fig. 15. Note that there is a different equation and a
different solution for $m_{\text{cosh}}$ for each value of $x_{3}$. The fits
are performed using the Levenberg-Marquardt’s method press2007numerical  and
the results are shown in Fig. 16.

Figure 16: Renormalized screening correlation functions in temperature units
$G_{P}^{s}(x_{3},T)/T^{3}$ (l.h.s) and $G_{AP}^{s}(x_{3},T)/T^{3}$ (r.h.s.)
and the result of the simultaneous fit. Top panels: Correlators measured on
the hadronic E250Nt24 ensemble. Middle panels: Correlators measured on the
chiral crossover E250Nt20 ensemble. Bottom panels: Correlators measured on the
high temperature E250Nt16 ensemble.

We fit the screening correlators $G_{P}^{s}(x_{3},T)$ and
$G_{AP}^{s}(x_{3},T)$ simultaneously using Eqs. (95)-(96). To assign fit
qualities in view of different fit windows, we apply the Akaike information
criterion Akaike:1998zah  and the model averaging method from Ref. Jay:2020jkz
. We assign a weight

$w_{i}=N\cdot\text{exp}\left[-\frac{1}{2}\left(\chi_{i}^{2}+2p_{i}-2n_{i}\right)\right]\,,$
(99)

where $p_{i}$ is the number of fit parameters and $n_{i}$ the number of data
points in the fit666Since we are performing simultaneous one-state fits
according to Eqs. (95)-(96) the number of fit parameters is three and we vary
only the starting points and lengths of the fit intervals. with minimized
$\chi_{i}^{2}$. $N$ is chosen such that the weights $w_{i}$ are normalized
according to $\sum_{i}w_{i}=1$. The results of the fits are shown in Fig. 16.

###  E.1 Alternative determination of the screening pion decay constant

Figure 17:  Upper panel: Screening correlator
$G_{A}^{s}(x_{3},T=154\,\text{MeV})$ together with the corresponding fit. Left
panel: Integrated autocorrelation time $\tau_{int}$ for $f_{\pi}$ determined
as described in Sec. III.2 Right panel : Integrated autocorrelation time
$\tau_{int}$ for $f_{\pi}$ determined from $G_{A}^{s}(x_{3},T)$.

We determine the screening pion decay constant on all three ensembles in Sec.
III.2. The method described therein relies on a simultaneous fit ansatz of the
screening pseudoscalar and axial-pseudoscalar correlator. As the latter
correlator is an order parameter for chiral symmetry restoration, we are
confronted with huge integrated autocorrelation times $\tau_{int}$ on the
E250Nt20 ensemble which is generated at a temperature close to the
pseudocritical temperature $T_{pc}$ [see Fig. 4]. Alternatively, one could
proceed as in Ref. Ce:2022dax  and extract the pion decay constant
simultaneously from the pseudoscalar and axial correlator. However, this fit
ansatz relies on the PCAC mass, which itself needs to be extracted from the
axial-pseudoscalar correlator. Thus, as a crosscheck, we determine the
screening decay constant on the E250Nt20 ensemble solely from the noisy axial
correlator $G_{A}^{s}$ [see top panel of Fig. 17]. The asymptotic form of this
correlator defines the screening pion decay constant $f_{\pi}$ [see Eq. (5)].
Consequently, using a cosh fit ansatz

$G_{A}^{s}(x_{3},T)\overset{\circ}{=}A_{1}^{2}\,\text{cosh}[m_{1}(x_{3}-L/2)]\,,$
(100)

the screening pion decay constant $f_{\pi}$ is obtained as

$f_{\pi}=A_{1}\sqrt{e^{m_{1}L/2}\,m_{1}}\,.$ (101)

We extract $f_{\pi}/T=0.179(24)$, compatible with our result
$f_{\pi}/T=0.172(12)$ [see Table 3].

##  Appendix F Free theory expression for spectral function of the axial
charge

The spectral representation for the axial-charge correlator reads

$\int\mathrm{d}^{3}x\;\langle
A_{0}^{a}(x)A_{0}^{b}(0)\rangle=\delta^{ab}\int_{0}^{\infty}\mathrm{d}\omega\;\frac{\cosh(\omega(\beta/2-x_{0}))}{\sinh(\omega\beta/2)}\rho_{A}(\omega,\bm{0}).$
(102)

For non-interacting quarks of mass $m$, the spectral function reads

$\rho_{A}(\omega,\bm{0})=\chi_{A}\,\omega\,\delta(\omega)+\frac{m^{2}N_{c}}{4\pi^{2}}\,\theta(\omega^{2}-4m^{2})\sqrt{1-\frac{4m^{2}}{\omega^{2}}}\,\tanh\left(\frac{\beta\omega}{4}\right),$
(103)

with the axial-charge susceptibility reading

$\chi_{A}=2\beta
N_{c}\int\frac{d^{3}p}{(2\pi)^{3}}\,\frac{\bm{p}^{2}}{E_{\bm{p}}^{2}}\,f^{F}_{\bm{p}}(1-f^{F}_{\bm{p}}).$
(104)

It is worth noting that the analogue of Eq. (103) for the isospin charge,

$\int\mathrm{d}^{3}x\;\langle
V_{0}^{a}(x)V_{0}^{b}(0)\rangle=\delta^{ab}\int_{0}^{\infty}\mathrm{d}\omega\;\frac{\cosh(\omega(\beta/2-x_{0}))}{\sinh(\omega\beta/2)}\rho_{00}(\omega,\bm{0}),$
(105)

reads

$\rho_{00}(\omega,\bm{0})=\chi_{s}\,\omega\,\delta(\omega)$ (106)

with

$\chi_{s}-\chi_{A}=2\beta
N_{c}m^{2}\int\frac{\mathrm{d}^{3}p}{(2\pi)^{3}}\,\frac{1}{E_{\bm{p}}^{2}}\,f^{F}_{\bm{p}}(1-f^{F}_{\bm{p}}).$
(107)

##  Appendix G Numerical values for temporal correlators

In this appendix we list the means and errors of the (anti)symmetrized
temporal correlators used in this work. For convenience, the $PA_{0}$ and
$V-A$ correlators are displayed in Figs. 18 and 20, respectively. The temporal
$A_{0}A_{0}$ correlator at $T=128$ MeV and its effective mass are displayed in
Fig. 19.

Figure 18: Ratios of the temporal thermal $(PA_{0})$-correlators and the
reconstructed correlator $(PA_{0})$-correlator. The errors on the ratio have
been estimated using the $\Gamma$ method in the implementation of the pyerrors
package introduced in Ref. Joswig:2022qfe . Top and middle panel: Result for
the E250Nt24 and E250Nt20 ensemble respectively. The blue band shows the
result from a correlated fit to the plateau. In both cases the
$(PA_{0})-$correlators have been $\mathcal{O}(a)$-improved. Bottom panel:
Result for the E250Nt16 ensemble. Since there is no clear plateau visible we
quote the result for $x_{0}/a=7$. In the high-temperature phase the
$(PA_{0})-$correlators have not been $\mathcal{O}(a)$-improved.

Figure 19: Both plots show results from the hadronic E250Nt24 ensemble. Left
panel: Effective mass plot for the cosh mass $m_{\text{cosh}}(x_{0})/T$ as a
function of the $x_{0}$-coordinate in temperature units, obtained from the
temporal axial correlation function at zero spatial momentum
$G_{A_{0}}(x_{0},T)$. Right panel: The blue bars correspond to the lattice
data of the temporal $G_{A_{0}}(x_{0},T)$-correlator in the hadronic phase,
while the red bars represent the prediction of the aforementioned correlator
in terms of the screening quantities $f_{\pi}$ and $m_{\pi}$. The green band
is the result of a direct one-state fit to the temporal correlator. From this
fit we obtain $\omega_{\mathbf{0}}/T=0.91(9)$.

Figure 20: L.h.s.: The improved thermal correlator for the difference “$(V-A)$”. The top panel shows the result in the hadronic phase, the middle panel at the finite temperature crossover and the bottom panel in the high-temperature phase. R.h.s.: The corresponding reconstructed correlators [following Eq. (67)] for the $\mathcal{O}(a)-$improved difference “$(V-A)$”. All correlators are renormalized and the corresponding ratios are shown in Fig. 11.  Table 7: Temporal symmetrized $V-A$ correlators projected to zero momentum. All errors quoted are purely statistical. The correlators are renormalized and quoted in units of temperature (i.e. scaled by a factor $N_{\tau}^{3}$). |  $T=128$ MeV |  $T=154$ MeV |  $T=192$ MeV
---|---|---|---
$x_{0}/a$ | $\frac{G_{(V-A)}(x_{0})}{T^{3}}$ | $\frac{G^{\text{\,rec}}_{(V-A)}(x_{0})}{T^{3}}$ | $\frac{G_{(V-A)}(x_{0})}{T^{3}}$ | $\frac{G^{\text{\,rec}}_{(V-A)}(x_{0})}{T^{3}}$ | $\frac{G_{(V-A)}(x_{0})}{T^{3}}$ | $\frac{G^{\text{\,rec}}_{(V-A)}(x_{0})}{T^{3}}$
3 | 2.7840(56) | 2.971(67) | 1.4447(32) | 1.761(40) | 0.73012(30) | 0.917(21)
4 | 0.5508(54) | 0.725(17) | 0.1571(30) | 0.461(11) | 0.06793(26) | 0.2554(58)
5 | 0.3230(50) | 0.482(11) | 0.0360(29) | 0.3220(68) | 0.00468(25) | 0.1889(42)
6 | 0.2760(44) | 0.4170(95) | 0.0215(28) | 0.2868(62) | -0.00288(24) | 0.1768(39)
7 | 0.2504(36) | 0.3755(85) | 0.0197(25) | 0.2673(60) | -0.00373(21) | 0.1725(39)
8 | 0.2295(31) | 0.3415(81) | 0.0191(25) | 0.2532(59) | -0.00373(24) | 0.1713(40)
9 | 0.2127(28) | 0.3133(79) | 0.0188(23) | 0.2441(59) |  |
10 | 0.2001(26) | 0.2918(72) | 0.0186(23) | 0.2415(59) |  |
11 | 0.1914(26) | 0.2780(72) |  |  |  |
12 | 0.1888(26) | 0.2741(74) |  |  |  |

## References

  * (1) Y. Aoki, G. Endrodi, Z. Fodor, S. Katz and K. Szabo, _The Order of the quantum chromodynamics transition predicted by the standard model of particle physics_ , _Nature_ 443 (2006) 675 [hep-lat/0611014].
  * (2) HotQCD collaboration, _Chiral crossover in QCD at zero and non-zero chemical potentials_ , _Phys. Lett. B_ 795 (2019) 15 [1812.08235].
  * (3) Wuppertal-Budapest collaboration, _Is there still any Tc mystery in lattice QCD? Results with physical masses in the continuum limit III_ , _JHEP_ 09 (2010) 073 [1005.3508].
  * (4) HotQCD collaboration, _Chiral Phase Transition Temperature in ( 2+1 )-Flavor QCD_ , _Phys. Rev. Lett._ 123 (2019) 062002 [1903.04801].
  * (5) E.V. Shuryak, _Physics of the pion liquid_ , _Phys. Rev. D_ 42 (1990) 1764.
  * (6) J. Goity and H. Leutwyler, _On the Mean Free Path of Pions in Hot Matter_ , _Phys.Lett._ B228 (1989) 517.
  * (7) J. Gasser and H. Leutwyler, _Thermodynamics of Chiral Symmetry_ , _Phys. Lett. B_ 188 (1987) 477.
  * (8) P. Gerber and H. Leutwyler, _Hadrons Below the Chiral Phase Transition_ , _Nucl. Phys. B_ 321 (1989) 387.
  * (9) A. Schenk, _Pion propagation at finite temperature_ , _Phys.Rev._ D47 (1993) 5138.
  * (10) D. Toublan, _Pion dynamics at finite temperature_ , _Phys.Rev._ D56 (1997) 5629 [hep-ph/9706273].
  * (11) R.D. Pisarski and M. Tytgat, _Propagation of cool pions_ , _Phys.Rev._ D54 (1996) 2989 [hep-ph/9604404].
  * (12) D.T. Son and M.A. Stephanov, _Pion propagation near the QCD chiral phase transition_ , _Phys. Rev. Lett._ 88 (2002) 202302 [hep-ph/0111100].
  * (13) D.T. Son and M.A. Stephanov, _Real time pion propagation in finite temperature QCD_ , _Phys. Rev. D_ 66 (2002) 076011 [hep-ph/0204226].
  * (14) L.D. McLerran and T. Toimela, _Photon and Dilepton Emission from the Quark - Gluon Plasma: Some General Considerations_ , _Phys.Rev._ D31 (1985) 545.
  * (15) P. Braun-Munzinger, V. Koch, T. Schäfer and J. Stachel, _Properties of hot and dense matter from relativistic heavy ion collisions_ , _Phys. Rept._ 621 (2016) 76 [1510.00442].
  * (16) R. Rapp and J. Wambach, _Chiral symmetry restoration and dileptons in relativistic heavy ion collisions_ , _Adv. Nucl. Phys._ 25 (2000) 1 [hep-ph/9909229].
  * (17) M. Laine, _NLO thermal dilepton rate at non-zero momentum_ , _JHEP_ 11 (2013) 120 [1310.0164].
  * (18) G. Aarts, C. Allton, S. Hands, B. Jäger, C. Praki and J.-I. Skullerud, _Nucleons and parity doubling across the deconfinement transition_ , _Phys. Rev. D_ 92 (2015) 014503 [1502.03603].
  * (19) G. Aarts, C. Allton, D. De Boni, S. Hands, B. Jäger, C. Praki et al., _Light baryons below and above the deconfinement transition: medium effects and parity doubling_ , _JHEP_ 06 (2017) 034 [1703.09246].
  * (20) B.B. Brandt, A. Francis, H.B. Meyer and D. Robaina, _Chiral dynamics in the low-temperature phase of QCD_ , _Phys. Rev. D_ 90 (2014) 054509 [1406.5602].
  * (21) B.B. Brandt, A. Francis, H.B. Meyer and D. Robaina, _Pion quasiparticle in the low-temperature phase of QCD_ , _Phys. Rev. D_ 92 (2015) 094510 [1506.05732].
  * (22) M. Cè, T. Harris, A. Krasniqi, H.B. Meyer and C. Török, _Aspects of chiral symmetry in QCD at T=128 MeV_ , _Phys. Rev. D_ 107 (2023) 054509 [2211.15558].
  * (23) A. Krasniqi, M. Cè, T. Harris, H. B. Meyer, C. Török and S. Ruhl, _A $(2+1)$-flavor lattice study of the pion quasiparticle in the thermal hadronic phase at physical quark masses_, _PoS_ Lattice2022 (2022) 0181.
  * (24) M. Lüscher, _Advanced lattice QCD_ , in _Les Houches Summer School in Theoretical Physics, Session 68: Probing the Standard Model of Particle Interactions_ , pp. 229–280, 2, 1998 [hep-lat/9802029].
  * (25) H.B. Meyer, _Transport Properties of the Quark-Gluon Plasma: A Lattice QCD Perspective_ , _Eur. Phys. J. A_ 47 (2011) 86 [1104.3708].
  * (26) S.R. Sharpe, _Applications of Chiral Perturbation theory to lattice QCD_ , in _Workshop on Perspectives in Lattice QCD_ , 7, 2006 [hep-lat/0607016].
  * (27) ALPHA collaboration, _Non-perturbative quark mass renormalisation and running in $N_{f}=3$ QCD_, _Eur. Phys. J. C_ 78 (2018) 387 [1802.05243].
  * (28) J. Bulava and S. Schaefer, _Improvement of $N_{f}$ = 3 lattice QCD with Wilson fermions and tree-level improved gauge action_, _Nucl. Phys. B_ 874 (2013) 188 [1304.7093].
  * (29) M. Bruno et al., _Simulation of QCD with N ${}_{f}=$ 2 $+$ 1 flavors of non-perturbatively improved Wilson fermions_, _JHEP_ 02 (2015) 043 [1411.3982].
  * (30) D. Mohler, S. Schaefer and J. Simeth, _CLS 2+1 flavor simulations at physical light- and strange-quark masses_ , _EPJ Web Conf._ 175 (2018) 02010 [1712.04884].
  * (31) M. Bruno, T. Korzec and S. Schaefer, _Setting the scale for the CLS $2+1$ flavor ensembles_, _Phys. Rev._ D95 (2017) 074504 [1608.08900].
  * (32) M. Cè et al., _Window observable for the hadronic vacuum polarization contribution to the muon $g-2$ from lattice QCD_, 2206.06582.
  * (33) M. Lüscher and S. Schaefer, _Lattice QCD with open boundary conditions and twisted-mass reweighting_ , _Comput. Phys. Commun._ 184 (2013) 519 [1206.2809].
  * (34) M.A. Clark and A.D. Kennedy, _Accelerating dynamical fermion computations using the rational hybrid Monte Carlo (RHMC) algorithm with multiple pseudofermion fields_ , _Phys. Rev. Lett._ 98 (2007) 051601 [hep-lat/0608015].
  * (35) D. Mohler and S. Schaefer, _Remarks on strange-quark simulations with Wilson fermions_ , _Phys. Rev. D_ 102 (2020) 074506 [2003.13359].
  * (36) S. Kuberski, _Low-mode deflation for twisted-mass and RHMC reweighting in lattice QCD_ , _Comput. Phys. Commun._ 300 (2024) 109173 [2306.02385].
  * (37) N. Madras and A.D. Sokal, _The Pivot algorithm: a highly efficient Monte Carlo method for selfavoiding walk_ , _J. Statist. Phys._ 50 (1988) 109.
  * (38) ALPHA collaboration, _Monte Carlo errors with less errors_ , _Comput. Phys. Commun._ 156 (2004) 143 [hep-lat/0306017].
  * (39) A. Ramos, _Automatic differentiation for error analysis of Monte Carlo data_ , _Comput. Phys. Commun._ 238 (2019) 19 [1809.01289].
  * (40) F. Joswig, S. Kuberski, J.T. Kuhlmann and J. Neuendorf, _pyerrors: A python framework for error analysis of Monte Carlo data_ , _Comput. Phys. Commun._ 288 (2023) 108750 [2209.14371].
  * (41) UKQCD collaboration, _Quark mass dependence of hadron masses from lattice QCD_ , _Phys. Rev. D_ 59 (1999) 074503 [hep-lat/9810021].
  * (42) P.A. Boyle, A. Juttner, C. Kelly and R.D. Kenway, _Use of stochastic sources for the lattice determination of light quark physics_ , _JHEP_ 08 (2008) 086 [0804.1501].
  * (43) ETM collaboration, _Dynamical Twisted Mass Fermions with Light Quarks: Simulation and Analysis Details_ , _Comput. Phys. Commun._ 179 (2008) 695 [0803.0224].
  * (44) G.S. Bali, S. Collins and A. Schafer, _Effective noise reduction techniques for disconnected loops in Lattice QCD_ , _Comput. Phys. Commun._ 181 (2010) 1570 [0910.3970].
  * (45) G. Backus and F. Gilbert, _The Resolving Power of Gross Earth Data_ , _Geophysical Journal International_ 16 (1968) 169 [https://academic.oup.com/gji/article-pdf/16/2/169/5891044/16-2-169.pdf].
  * (46) N. Wen, X. Cao, J. Chao and H. Liu, _Neutral pion masses within a hot and magnetized medium in a lattice-improved soft-wall AdS/QCD model_ , 2402.06239.
  * (47) D. Goderidze, A. Friesen and Y. Kalinovsky, _Pion damping width and pion spectral function in hot pion gas_ , _Int. J. Mod. Phys. A_ 37 (2022) 2250135 [2205.11436].
  * (48) A. Bazavov et al., _Meson screening masses in (2+1)-flavor QCD_ , _Phys. Rev. D_ 100 (2019) 094510 [1908.09552].
  * (49) S. Borsanyi, Z. Fodor, S.D. Katz, S. Krieg, C. Ratti and K. Szabo, _Fluctuations of conserved charges at finite temperature from lattice QCD_ , _JHEP_ 01 (2012) 138 [1112.4416].
  * (50) R. Hagedorn, _How We Got to QCD Matter from the Hadron Side: 1984_ , _Lect. Notes Phys._ 221 (1985) 53.
  * (51) F. Karsch, K. Redlich and A. Tawfik, _Hadron resonance mass spectrum and lattice QCD thermodynamics_ , _Eur. Phys. J. C_ 29 (2003) 549 [hep-ph/0303108].
  * (52) B.B. Brandt, A. Francis, B. Jäger and H.B. Meyer, _Charge transport and vector meson dissociation across the thermal phase transition in lattice QCD with two light quark flavors_ , _Phys. Rev. D_ 93 (2016) 054510 [1512.07249].
  * (53) M. Davier, A. Hocker and Z. Zhang, _The Physics of hadronic tau decays_ , _Rev.Mod.Phys._ 78 (2006) 1043 [hep-ph/0507078].
  * (54) J.I. Kapusta and E.V. Shuryak, _Weinberg type sum rules at zero and finite temperature_ , _Phys. Rev._ D49 (1994) 4694 [hep-ph/9312245].
  * (55) P.M. Hohler and R. Rapp, _Is $\rho$-Meson Melting Compatible with Chiral Restoration?_, _Phys.Lett._ B731 (2014) 103 [1311.2921].
  * (56) M. Dey, V.L. Eletsky and B.L. Ioffe, _Mixing of vector and axial mesons at finite temperature: an Indication towards chiral symmetry restoration_ , _Phys. Lett. B_ 252 (1990) 620.
  * (57) V.L. Eletsky and B.L. Ioffe, _On the current correlators in QCD at finite temperature_ , _Phys. Rev. D_ 47 (1993) 3083 [hep-ph/9302298].
  * (58) V.L. Eletsky and B.L. Ioffe, _Next-to-leading order temperature corrections to correlators in QCD_ , _Phys. Rev. D_ 51 (1995) 2371 [hep-ph/9405371].
  * (59) H.B. Meyer, _The bulk channel in thermal gauge theories_ , _JHEP 04 (2010) 099_ (2010) [arXiv:1002.3343].
  * (60) J. Gasser and H. Leutwyler, _Light Quarks at Low Temperatures_ , _Phys. Lett. B_ 184 (1987) 83.
  * (61) G. Aarts and A. Nikolaev, _Electrical conductivity of the quark-gluon plasma: perspective from lattice QCD_ , _Eur. Phys. J. A_ 57 (2021) 118 [2008.12326].
  * (62) D. Bernecker and H.B. Meyer, _Vector Correlators in Lattice QCD: Methods and applications_ , _Eur.Phys.J._ A47 (2011) 148 [1107.4388].
  * (63) S. Caron-Huot, _Asymptotics of thermal spectral functions_ , _Phys. Rev. D_ 79 (2009) 125009 [0903.3958].
  * (64) S. Weinberg, _Precise relations between the spectra of vector and axial vector mesons_ , _Phys. Rev. Lett._ 18 (1967) 507.
  * (65) S. Kuberski, _Muon $g-2$: Lattice calculations of the hadronic vacuum polarization_, _PoS_ LATTICE2023 (2024) 125 [2312.13753].
  * (66) M. Galassi et al., “Gnu scientific library reference manual.” https://www.gnu.org/software/gsl/, 2019.
  * (67) L. Fousse, G. Hanrot, V. Lefèvre, P. Pélissier and P. Zimmermann, _Mpfr: A multiple-precision binary floating-point library with correct rounding_ , _ACM Trans. Math. Softw._ 33 (2007) 13–es.
  * (68) T. Granlund and the GMP development team. https://gmplib.org/, 2020.
  * (69) ALPHA collaboration, _Nonperturbative results for the coefficients b(m) and b(a) - b(P) in O(a) improved lattice QCD_ , _Nucl. Phys. B_ 595 (2001) 44 [hep-lat/0009021].
  * (70) ALPHA collaboration, _Non-perturbative improvement of the axial current in $N_{f}$=3 lattice QCD with Wilson fermions and tree-level improved gauge action_, _Nucl. Phys. B_ 896 (2015) 555 [1502.04999].
  * (71) ALPHA collaboration, _Ward identity determination of $Z_{\mathrm{S}}/Z_{\mathrm{P}}$ for $N_{\mathrm{f}}=3$ lattice QCD in a Schrödinger functional setup_, _Eur. Phys. J. C_ 80 (2020) 765 [2005.01352].
  * (72) P. Korcyl and G.S. Bali, _Non-perturbative determination of improvement coefficients using coordinate space correlators in $N_{f}=2+1$ lattice QCD_, _Phys. Rev. D_ 95 (2017) 014505 [1607.07090].
  * (73) M. Dalla Brida, T. Korzec, S. Sint and P. Vilaseca, _High precision renormalization of the flavour non-singlet Noether currents in lattice QCD with Wilson quarks_ , _Eur. Phys. J. C_ 79 (2019) 23 [1808.09236].
  * (74) RQCD collaboration, _Octet baryon isovector charges from Nf=2+1 lattice QCD_ , _Phys. Rev. D_ 108 (2023) 034512 [2305.04717].
  * (75) A. Gerardin, T. Harris and H.B. Meyer, _Nonperturbative renormalization and $O(a)$-improvement of the nonsinglet vector current with $N_{f}=2+1$ Wilson fermions and tree-level Symanzik improved gauge action_, _Phys. Rev. D_ 99 (2019) 014519 [1811.08209].
  * (76) M. Cè, A. Gérardin, G. von Hippel, H.B. Meyer, K. Miura, K. Ottnad et al., _The hadronic running of the electromagnetic coupling and the electroweak mixing angle from lattice QCD_ , _JHEP_ 08 (2022) 220 [2203.08676].
  * (77) A. Agadjanov, D. Djukanovic, G. von Hippel, H.B. Meyer, K. Ottnad and H. Wittig, _Nucleon Sigma Terms with Nf=2+1 Flavors of O(a)-Improved Wilson Fermions_ , _Phys. Rev. Lett._ 131 (2023) 261902 [2303.08741].
  * (78) A. Krasniqi, M. Cè, T. Harris, R.J. Hudspith, H.B. Meyer and C. Török, _The thermal photon emissivity at the QCD chiral crossover from imaginary momentum correlators_ , _PoS_ LATTICE2023 (2024) 180 [2401.05951].
  * (79) W. Press, S. Teukolsky, W. Vetterling and B. Flannery, _Numerical Recipes: The Art of Scientific Computing_ , Cambridge University Press, 3 ed. (2007).
  * (80) H. Akaike, _Information Theory and an Extension of the Maximum Likelihood Principle_ , (New York), Springer Science+Business Media (1998), DOI.
  * (81) W.I. Jay and E.T. Neil, _Bayesian model averaging for analysis of lattice field theory results_ , _Phys. Rev. D_ 103 (2021) 114502 [2008.01069].