HotQCD Collaboration

# Lattice QCD estimates of thermal photon production from the QGP

Sajid Alia,f, Dibyendu Balaa, Anthony Francisb, Greg Jacksonc,d, Olaf
Kaczmareka, Jonas Turnwalde, Tristan Uedinga, Nicolas Winke aFakultät für
Physik, Universität Bielefeld, D-33615 Bielefeld, Germany
bInstitute of Physics, National Yang-Ming Chiao Tung University, 30010
Hsinchu, Taiwan
cInstitute for Nuclear Theory, Box 351550, University of Washington, Seattle,
WA 98195-1550, United States
dSUBATECH (Nantes Université, IMT Atlantique, IN2P3/CNRS), 4 rue Alfred
Kastler, La Chantrerie BP 20722, 44307 Nantes, France
eInstitut für Kernphysik, Technische Universität Darmstadt, D-64289 Darmstadt,
Germany
fGovernment College University Lahore, Department of Physics, Lahore 54000,
Pakistan

###### Abstract

Thermal photons produced in heavy-ion collision experiments are an important
observable for understanding quark-gluon plasma (QGP). The thermal photon rate
from the QGP at a given temperature can be calculated from the spectral
function of the vector current correlator. Extraction of the spectral function
from the lattice correlator is known to be an ill-conditioned problem, as
there is no unique solution for a spectral function for a given lattice
correlator with statistical errors. The vector current correlator, on the
other hand, receives a large ultraviolet contribution from the vacuum, which
makes the extraction of the thermal photon rate difficult from this channel.
We therefore consider the difference between the transverse and longitudinal
part of the spectral function, only capturing the thermal contribution to the
current correlator, simplifying the reconstruction significantly. The lattice
correlator is calculated for light quarks in quenched QCD at $T=470~{}$MeV
($\sim 1.5\,T_{c}$), as well as in 2+1 flavor QCD at $T=220~{}$MeV ($\sim
1.2\,T_{pc}$) with $m_{\pi}=320$ MeV. In order to quantify the non-
perturbative effects, the lattice correlator is compared with the
corresponding $\text{NLO}+\text{LPM}^{\text{LO}}$ estimate of correlator. The
reconstruction of the spectral function is performed in several different
frameworks, ranging from physics-informed models of the spectral function to
more general models in the Backus-Gilbert method and Gaussian Process
regression. We find that the resulting photon rates agree within errors.

##  I Introduction

One of the important predictions of quantum chromodynamics (QCD), the theory
of strong interactions, is that the QCD coupling constant decreases with the
energy scale. This characteristic makes QCD strongly coupled (non-
perturbative) at low energies, and weakly coupled at very high energies. As a
result, non-perturbative studies of QCD are important at low energies. First-
principles lattice QCD equation-of-state calculations show that, as the
temperature increases, QCD matter undergoes a crossover transition from a
hadronic phase to a quark-gluon plasma phase (QGP) at a pseudo-critical
temperature of $T_{pc}={156.5\pm 1.5}~{}\textrm{MeV}$ [1, 2, 3].

Experimentally, this phase has been studied in ultrarelativistic heavy-ion
collisions at RHIC and LHC. Photons and dileptons produced in a heavy-ion
collision emerge directly from the point of creation without any further
interaction with their surroundings [4, 5, 6, 7, 8]. As a result, they carry
information about the local region around their place of origin. However, in
experiment, one measures the total number of photons produced from different
stages of the QCD matter evolution after the collision. Separating photons
produced at different stages of evolution is a challenging task, as the vast
majority ($\sim 80...90$%) of photons emitted because of hadronic decays (e.g.
$\pi^{0}\to\gamma\,\gamma$) during the hadronization stage [9]. These photons
are known as the hadronic background. The remaining photons, which are
produced directly from the QGP or before its formation, are called direct
photons. Extraction of these direct photons from the hadronic background
requires sophisticated techniques and careful analysis [10].

The Au-Au collisions at PHENIX [11] and Pb-Pb collisions at ALICE [12] have
revealed large yields of direct photons. The high-$p_{T}$ segment of the
photon spectra (prompt photons) can be related to the yield of photons
generated in proton-proton collisions, scaled according to the number of
participating binary collisions. The excess photons observed in the
low-$p_{T}$ region of the spectra may be tentatively attributed to in-medium
plasma effects. Measurements of direct photons and dileptons also reveal a
large elliptical flow, quantified by the coefficients $v_{2}$ [13, 14].
However, understanding the origin of this phenomenon from first principles
poses a significant theoretical challenge [15].

Any attempt to explain these observables must involve integrating the
(multistage) hydrodynamic expansion of the plasma with the differential photon
or dilepton production rate embedded at each stage of its evolution.
Simulations have revealed that the ‘slopes’ of the $p_{T}$ spectrum (for
photons) and the invariant mass spectrum (for dileptons) can be systematically
related to an average temperature of the fireball [16, 17, 18]. The
anisotropic flow of direct photons and dileptons has been shown to be
sensitive to the initial conditions and shear viscosity [19, 20, 21]. Since
electromagnetic probes are generated at every stage of plasma evolution, they
may supply useful constraints on the viscosity and other transport
coefficients. One of the important ingredients in such analyses is a precise
understanding of the photon production mechanism from each stage of the plasma
evolution. Since hydrodynamic expansion assumes local thermal equilibrium,
computation of the photon production rate at each stage requires the use of
equilibrium thermal field theory calculations.

In this paper, we aim to estimate the thermal photon rate from the QGP using
lattice QCD. A major challenge is posed by the fact that the thermal
production of photons requires real-time information, which is, in the context
of lattice QCD, only accessible via analytic continuation. This analytic
continuation is formally an ill-conditioned inverse operation [22, 23].
Therefore, we apply different methods to address this problem, ranging from
physics informed model fits and the Backus-Gilbert method to Gaussian process
regression.

Using these methods, we can calculate the thermal photon rate from the lattice
correlator for light quarks in two different scenarios: Quenched QCD at a
temperature of $470\text{\,}\mathrm{MeV}$ ($1.5\,T_{c}$), and (2+1)-flavor QCD
at a temperature of $220\text{\,}\mathrm{MeV}$ ($1.22\,T_{pc}$). For the
quenched QCD scenario, we have extrapolated the correlator to the continuum,
which means that our prediction for the photon production rate is also at the
continuum. For the full QCD scenario, our results are at a finite lattice
spacing.

The paper is organized into the following sections: In Section II, we will
briefly discuss the thermal photon rate and the correlator that we used to
calculate it. In Section III, we present the lattice details and perform the
continuum extrapolation of quenched lattice data. In Section IV, we compare
the lattice correlator with the correlator obtained from perturbative spectral
function, to quantify the non-perturbative effects. We undertake the spectral
reconstruction using multiple techniques in Section V. Finally, in Section
VII, we compare the photon production rates obtained from different techniques
and provide the final value for the photon production rate.

##  II Thermal photon rate and Spectral function

In a thermalized plasma at temperature $T$, the thermal photon production rate
is defined as the number of photons radiated from the plasma per unit time and
per unit volume. To leading order in the QED coupling $\alpha_{\rm em}$, this
rate can be calculated as

$\frac{{\rm d}\Gamma_{\gamma}}{{\rm d}^{3}{\vec{k}}}\;=\;-\,\frac{\alpha_{\rm
em}n_{b}(k)}{2\pi^{2}k}\Big{\\{}{\textstyle\sum_{i=1}^{N_{f}}Q_{i}^{2}}\Big{\\}}g^{\mu\nu}\rho_{\mu\nu}\big{(}\,\omega=|\vec{k}|,\,\vec{k}\,\big{)}\;,$
(1)

where $n_{b}(\omega)=1/(e^{\omega/T}-1)$ is the Bose distribution function,
and $Q_{i}$ is the quark charge in units of $e$ ($i$ denotes flavor)[24]. We
take the Minkowski space metric tensor as $g_{\mu\nu}={\rm diag}(1,-1,-1,-1)$.

The spectral function $\rho_{\mu\nu}(\omega,\vec{k})$ contains the full
information about QCD dynamics and can be calculated through analytic
continuation from the Matsubara frequency modes,111These discrete energies are
denoted $\omega_{n}=2\pi Tn\,$, where $n\in\mathbb{Z}$. i.e.,
$\rho_{\mu\nu}(\omega,\vec{k})=\mathrm{Im}[G^{E}_{\mu\nu}(\omega_{n}\to-i\omega,\vec{k})]$,
where the Euclidean current correlator is given by

$G^{E}_{\mu\nu}(\omega_{n},\vec{k})=\int{\rm d}^{3}\vec{x}\int_{0}^{\beta}{\rm
d}\tau\,e^{i(\omega_{n}\tau-\vec{k}.\vec{x})}\langle
J_{\mu}(\tau,\vec{x})J_{\nu}(0,\vec{0})\rangle\,,$ (2)

where $\beta=1/T$ is the temporal extent and the electromagnetic current
$J_{\mu}(\tau,\vec{x})=\bar{\psi}(\tau,\vec{x})\gamma_{\mu}\psi(\tau,\vec{x})\,$.

On the lattice we calculate the Euclidean correlation function
$G^{E}_{\mu\nu}(\tau,\vec{k})=T\,\sum_{n}\exp(i\omega_{n}\tau)G^{E}_{\mu\nu}(\omega_{n},\vec{k})$.
This correlator can be written in terms of the spectral function as

$G^{E}_{\mu\nu}(\tau,\vec{k})=\int_{0}^{\infty}\frac{{\rm
d}\omega}{\pi}\rho_{\mu\nu}(\omega,\vec{k})\frac{\cosh[\omega(\tau-\frac{1}{2T})]}{\sinh[\frac{\omega}{2T}]}\,.$
(3)

As mentioned before, extracting the spectral function from the Euclidean
correlator is an ill-conditioned problem. Therefore, when inverting Equation 3
naively, any errors in the lattice data lead to exponentially large errors in
the corresponding spectral function. There are, consequently, many possible
spectral functions that could reproduce the lattice data. The incorporation of
additional, physics motivated, information is needed in order to ‘regularize’
the problem and achieve reasonable reconstructions.

The contracted (or “vector channel”) spectral function
$g^{\mu\nu}\rho_{\mu\nu}$, appearing in Equation 1, is proportional to
$\omega^{2}$ in the ultraviolet (UV) domain. As a result, the lattice
correlator, determined by the integration in Equation 3, receives a
significant contribution from the UV part of the spectral function. This fact
complicates the reconstruction process in general, as the thermal photon rate
is dominated by the infrared (IR) part of the spectral function. Therefore,
following Ref. [25], we will estimate the thermal photon rate from the
difference between the transverse and longitudinal (T-L) spectral functions.
This approach has the advantage in that the UV part of the spectral function
is heavily suppressed, meaning that the corresponding Euclidean correlator is
governed by the IR part of the spectral function, resulting in more reliable
reconstructions.

In order to obtain the T-L part, the spectral function $\rho_{\mu\nu}$ is
decomposed in terms of $\rho_{T}$ (the transverse component) and $\rho_{L}$
(the longitudinal component) of the spectral function as

$\rho_{\mu\nu}(\omega,\vec{k})=P^{T}_{\mu\nu}\rho_{T}(\omega,\vec{k})+P^{L}_{\mu\nu}\rho_{L}(\omega,\vec{k})\,.$
(4)

Here, $P^{T}_{\mu\nu}$ and $P^{L}_{\mu\nu}$ are the projection operators,
given explicitly by

$\begin{split}&P^{T}_{ij}(\omega,\vec{k})\;=\;-g_{ij}-\frac{k_{i}k_{j}}{\vec{k}^{2}}\,,\\\\[4.30554pt]
&P^{T}_{0i}\;=\;P^{T}_{i0}\;=\;P^{T}_{00}\;=\;0\,,\\\\[4.30554pt]
&P^{L}_{\mu\nu}(\omega,\vec{k})\;=\;-g_{\mu\nu}+\frac{K_{\mu}K_{\nu}}{K^{2}}-P^{T}_{\mu\nu}(\omega,\vec{k})\,,\end{split}$
(5)

where $K=(\omega,\vec{k})$ is the associated four-momentum. The projection
operators satisfy the following relations: $P_{T}^{2}=P_{T}$,
$P_{L}^{2}=P_{L}$, and $P_{T}P_{L}=P_{L}P_{T}=0\,$. In terms of these
components, the vector channel spectral function is given by
$-\rho^{\mu}_{\mu}=2\rho_{T}+\rho_{L}$. The photon production rate is
therefore proportional to $\rho_{T}(\omega=|\vec{k}|,\vec{k})$, as on the
light cone longitudinal part of the spectral function vanishes,
$\rho_{L}(\omega=|\vec{k}|,\vec{k})=0\,$. This allows us to calculate the
photon production rate using the T-L spectral function as

$\rho_{H}(\omega,\vec{k})\;=\;2\big{\\{}\,\rho_{T}(\omega,\vec{k})-\rho_{L}(\omega,\vec{k})\,\big{\\}}\,.$
(6)

At zero temperature, because of the Ward identity and Lorentz invariance, the
T-L spectral function vanishes, $\rho_{H}(\omega,\vec{k})=0$. Therefore, non-
zero $\rho_{H}(\omega,\vec{k})$ displays purely thermal effects. Using the
Operator Product Expansion (OPE), it has also been shown that in the domain
where $\omega\gg k,\pi\,T$, this spectral function behaves asymptotically like
$1/\omega^{4}$ [26, 27]. One consequence of the OPE, is that this spectral
function satisfies a sum rule,

$\int_{0}^{\infty}\,{\rm d}\omega\,\omega\,\rho_{H}(\omega)=0.\,$ (7)

We should note that the photon production rate Equation 1 is often
(equivalently) expressed as

$\frac{{\rm d}\Gamma_{\gamma}}{{\rm d}^{3}{\vec{k}}}\;=\;\frac{\alpha_{\rm
em}n_{b}(k)\chi_{q}}{\pi^{2}}\Big{\\{}{\textstyle\sum_{i=1}^{N_{f}}Q_{i}^{2}}\Big{\\}}D_{\text{eff}}(k)\,,$
(8)

where the effective diffusion coefficient is defined by

$D_{\text{eff}}(k)\;\equiv\;\frac{\rho_{H}\big{(}|\vec{k}|,k\big{)}}{2\chi_{q}|\vec{k}|}.$
(9)

Here $\chi_{q}$ is the quark number susceptibility (defined later in Equation
12). In the hydrodynamic regime of large wavelengths, $D_{\text{eff}}(k)$
approaches the well-known diffusion coefficient $D$, i.e. $\lim_{k\to
0}D_{\text{eff}}(k)=D$.

##  III Lattice details

In lattice QCD calculations, space-time has finite spatial extent ($L$) and
discrete momenta, given by $\vec{k}=\left(\frac{2\pi n_{x}}{L},\frac{2\pi
n_{y}}{L},\frac{2\pi n_{z}}{L}\right)$, where
$(n_{x},n_{y},n_{z})\in\mathbb{Z}$. Correlation functions
$G_{xx}(\tau,\vec{k})$, $G_{yy}(\tau,\vec{k})$, $G_{zz}(\tau,\vec{k})$, and
$G_{\tau\tau}(\tau,\vec{k})$ are computed on the lattice. If the momentum is
in the the $x$ direction i.e. $\vec{k}=(k,0,0)$, then the correlator for
$\rho_{H}$ is given by [27]:

$G_{H}(\tau,k_{x}=k)\ =\
G_{yy}(\tau,k_{x}=k)\,+\,G_{zz}(\tau,k_{x}=k)\,-\,2\big{(}\,G_{xx}(\tau,k_{x}=k)-G_{\tau\tau}(\tau,k_{x}=k)\,\big{)}.$
(10)

The lattice correlator is symmetric under $\tau\to\beta-\tau$ and due to
rotational invariance,
$G_{H}(\tau,k_{x}=k)=G_{H}(\tau,k_{y}=k)=G_{H}(\tau,k_{z}=k)$. Therefore, we
considered the following quantity as an estimate of the T-L correlator:

$G_{H}(\tau,k)\ =\
\frac{1}{6}\Big{\\{}\,G_{H}(\tau,k_{x}=k)+G_{H}(\tau,k_{y}=k)+G_{H}(\tau,k_{z}=k)+(\tau\to\beta-\tau)\,\Big{\\}}.$
(11)

$N_{f}$ | $T/T_{c}$ | $\beta_{0}$ | $\kappa$ | $N_{\sigma}^{3}\times N_{\tau}$ | confs |  $a^{-1}$[GeV]
---|---|---|---|---|---|---
$0$ | $1.5$ | $7.394$ | $0.13407$ | $120^{3}\times 30$ | $1950$ | $14.1$
$7.192$ | $0.13440$ | $\ 96^{3}\times 24$ | $2000$ | $11.3$
$7.035$ | $0.13467$ | $\ 80^{3}\times 20$ | $1824$ | $9.39$
$2+1$ | $1.22$ | $8.249$ | $0.13515$ | $\ 96^{3}\times 32$ | $1750$ | $7.04$
Table 1:  Details of the lattices studied in this paper, where $\beta_{0}$ is
the coefficient of the Wilson plaquette and $\kappa$ is the hopping parameter.
The critical temperature $T_{c}=$313\text{\,}\mathrm{MeV}$$ for $N_{f}=0$ case
and the pseudo critical temperature is $T_{pc}=$180\text{\,}\mathrm{MeV}$$ for
$N_{f}=2+1$ case.

We calculated the $G_{H}$ correlator using Clover-improved Wilson fermions for
quenched configurations at a temperature of $\sim 1.5\,T_{c}$ and (2+1)-flavor
QCD configurations at a temperature of $\sim 1.2\,T_{pc}$. The scale for
quenched QCD has been set by $r_{0}=$0.47\text{\,}\mathrm{fm}$$ [28] and the
$T_{c}$ is obtained by $T_{c}=0.7457/r_{0}$ [29]. In the case of full QCD, the
scale is determined by $r_{1}=$0.3106\text{\,}\mathrm{fm}$$ and
$T_{pc}=$180\text{\,}\mathrm{MeV}$$ was obtained from disconnected chiral
susceptibility, cf. Ref. [30], as in Ref. [31].

We generated the quenched configurations using a standard Wilson gauge action
with heat-bath and overrelaxation updates. Each configuration was separated by
500 sweeps, with each sweep corresponding to one heat bath followed by four
overrelaxation steps. We generated pure gauge configurations in this way for
three lattice spacings, which allows us to perform a continuum extrapolation
of the correlators. For the valence sector, we selected a hopping parameter
$\kappa$ value very close to the $\kappa_{c}$ value, which was obtained by
cubic spline interpolation of $\kappa_{c}$ values given in Ref. [32]. The
Clover improvement constant, also known as the Sheikholeslami–Wohlert
coefficient ($c_{SW}$), was estimated from the non-perturbative
parametrization of Ref. [32].

For full QCD, we used gauge field configurations generated by the HotQCD
collaboration using Highly Improved Staggered Quark (HISQ) [33] and a tree-
level improved Lüscher-Weisz gauge action [34, 35] using the RHMC algorithm.
We tuned the $\kappa$ value to match the pion mass, which in this study was
unphysical and equal to $320\text{\,}\mathrm{MeV}$ with $m_{l}=m_{s}/5$. The
$c_{SW}$ coefficient was taken to be the tadpole-improved tree-level value,
i.e. $c_{SW}=1/u_{0}^{3}$, where $u_{0}$ is the fourth root of the plaquette
expectation value. In all cases considered, the quark mass was significantly
smaller than the temperature being studied. More detailed information
regarding the lattice parameters can be found in Table 1.

Figure 1:  Euclidean correlator from lattice simulations in the T-L channel,
cf. Equation 13, for quenched QCD. Continuum extrapolation for various values
of $\tau T$ at a fixed temperature appear on the left, and a comparison of the
continuum correlator with the correlator calculated on finite lattice spacing
for various spatial momenta is on the right.

Our implementation of the current operator on the lattice is not $O(a)$
improved, but as shown in Ref. [36] the effect of this improvement is
negligible in the coupling range we are working on. Therefore, we expect the
correction to the continuum result should be $O(a^{2})$ in both quenched and
full QCD lattice correlator.

The correlator $G_{H}$ calculated on the lattice requires renormalization
prior to performing the continuum limit. This correlator is multiplicatively
renormalizable, which allows us to extrapolate to the continuum by considering
the quantity $G_{H}/(2\chi_{q}T)$. Here, $\chi_{q}$ denotes the quark number
susceptibility and is defined as

$\chi_{q}=\int_{0}^{1/T}\,{\rm d}\tau\,G_{\tau\tau}(\tau,\vec{0})\,.$ (12)

Since $\chi_{q}$ and $G_{H}$ have the same renormalization coefficient, the
ratio $G_{H}/(2\chi_{q}T)$ does not require renormalization. To perform the
continuum extrapolation, we need lattice data for various lattice spacings at
a fixed value of $\tau T$ and which may not be available for coarser lattices.
To address this, we used cubic spline interpolation of the courser lattices to
estimate the correlator at the same value of $\tau T$ for different lattice
spacings. We determined the error on the interpolated data points through a
bootstrap procedure. This involved performing the interpolation on each
bootstrap sample. Once we estimated $G_{H}/(2\chi_{q}T)$ with the error at the
same $\tau T$ for different lattice spacing, we performed the continuum
extrapolation using the following ansatz:

$\frac{G_{H}}{2\chi_{q}T}(a)=\frac{G_{H}}{2\chi_{q}T}(a=0)+\frac{b}{N_{\tau}^{2}}\,.$
(13)

We illustrated the fitting of the above ansatz for a few values of $\tau T$ in
left panel of Figure 1. Since this is a linear fit, the error bar on the
continuum-extrapolated data points has been estimated using Gaussian error
propagation. The continuum extrapolated correlator along with the correlator
at finite lattice spacing is shown in the right panel of the same figure. We
observe that the cut-off effect on this correlator at the momenta $k/T=1.57$
and $3.14$ is barely visible; however, with increasing momentum, cut-off
effects start showing up when $\tau T\lesssim 0.25$. In contrast, the vector
channel shows a significant cut-off effect due to the presence of a large UV
component [37]. This also indicates that T-L correlators are dominated by the
IR part of the spectral function.

##  IV Comparison with perturbative estimate

Figure 2:  Perturbative spectral function in the T-L channel computed using
the procedure from Ref. [38] at $\text{NLO}+\text{LPM}^{\text{LO}}$, as a
function of $\omega$ for several values of spatial momenta.
($\text{LPM}^{\text{LO}}$ means that near the light cone LPM resummation is
performed at leading order.) Here the renormalization scale was fixed at
$\mu=\mu_{\rm opt}$ as described in the main text. The left panel shows the
spectral function with no dynamical quarks at $T=$470\text{\,}\mathrm{MeV}$$,
and the right panel shows the spectral function in $N_{f}=2+1$ flavor QCD at
$T=$220\text{\,}\mathrm{MeV}$\,$.

Figure 3:  Euclidean correlator in T-L channel, calculated from the
$\text{NLO}+\text{LPM}^{\text{LO}}$perturbative spectral function and the
corresponding lattice measurements for several momenta. (As mentioned in the
main text, the actual lattice data are multiplied by $\chi_{q}=0.897T^{2}$ for
the case of $N_{f}=0$ and $\chi_{q}=0.842T^{2}$ for $N_{f}=2+1$ in order to
plot $G_{H}/T^{3}$.) Solid lines depict the perturbative results computed from
Equation 3, while data points represent lattice correlators for quenched on
the left and full QCD on the right.

Perturbation theory is a useful tool for calculating spectral functions, which
becomes rigorous at high temperature thanks to asymptotic freedom, and
significant progress has been made in using it for this purpose [39, 40, 38,
41, 42]. However, at the phenomenologically interesting temperature range, QCD
is strongly coupled and weak-coupling calculations may not be sufficient.
Therefore, non-perturbative calculations are necessary to access the photon
rate. In this section, we compare the lattice correlator with the one obtained
from the perturbative spectral function to judge the importance of non-
perturbative effects.

At leading order (LO), the spectral function is determined from the process
$q\bar{q}\to\gamma^{*}$. At the next-to-leading order (NLO), the spectral
function receives contributions from well-known processes such as Compton
scattering ($qg\to\gamma^{*}q$), pair annihilation ($q\bar{q}\to\gamma^{*}g$),
various interference terms between them, as well as virtual loop corrections
to the LO process [43, 44, 45]. On the light cone, the LO spectral function is
zero, but at NLO, it exhibits a well-known logarithmic singularity and a
finite discontinuous part. Far away from the light cone
($M^{2}=|\omega^{2}-k^{2}|\mathrel{\hbox to0.0pt{\lower
2.5pt\hbox{$\sim$}\hss}\raise 1.99997pt\hbox{$>$}}(\pi T)^{2}$), perturbative
calculations are well-behaved [39, 40]. However, the singularity of the NLO
spectral function at the light cone indicates the breakdown of the naive
perturbative calculation, necessitating resummation to all orders of a certain
class of diagrams [46, 47]. Therefore, a separate treatment is required to
calculate the spectral function close to the light cone.

When $M^{2}\mathrel{\hbox to0.0pt{\lower 2.5pt\hbox{$\sim$}\hss}\raise
1.99997pt\hbox{$<$}}g^{2}T^{2}$, the hard thermal loop (HTL) resummation
introduces an asymptotic (thermal) mass on the internal quark propagator in
$2\to 2$ NLO processes, which regulates the logarithmic singularity at the
light cone [48, 49]. Additionally, multiple scatterings due to soft gluons can
occur with the quark that emits the photon.222These interactions can happen
within the formation time of the photon and end up contributing significantly
to the rate. It is also mandatory to resum such scatterings beyond leading-
logarithmic accuracy. To include this contribution, a technique called Landau-
Pomeranchuk-Migdal (LPM) resummation is used in this regime [41, 42]. This
resummation removes the singularity and the discontinuity of the spectral
function at the light cone and therefore makes the spectral function smooth
across the light cone.

We computed the perturbative spectral function and its corresponding
correlator for the temperatures and momenta used in our lattice study,
following the calculation described in Ref. [38]. In that paper, leading order
LPM resummation near $M^{2}\sim g^{2}T^{2}$ was performed by solving a two-
dimensional Schrödinger equation with a LO light cone potential [50]. The
spectral function in the intermediate region is obtained by interpolating
between the $M^{2}\sim T^{2}$ and $M^{2}\sim g^{2}T^{2}$ region. The
renormalization scale chosen for the calculation is $\mu_{\rm
opt}=\sqrt{M^{2}+(\zeta\pi T)^{2}}$ with $\zeta^{(N_{f}=0)}=1$ and
$\zeta^{(N_{f}=3)}=2\,$. This choice is motivated by the fact that, far from
the light cone, the relevant scale is $K^{2}=\omega^{2}-k^{2}$ (the photon’s
virtuality), whereas near the light cone, the physics is that of dimensionally
reduced EQCD, and thus the relevant scale is temperature [51]. The running
coupling then has been calculated by using $\Lambda^{{\overline{\rm
MS}}}=0.8\,T_{c}$ for quenched [29] and $\Lambda^{{\overline{\rm MS}}}=339$
MeV $\simeq 1.88\,T_{pc}$ for full QCD [52] (these choices have uncertainties
$\sim 10...25$%).

The resulting spectral function is plotted for various momenta in Figure 2 for
both quenched and full QCD. As expected, the spectral function is UV
suppressed, and it has been shown [38] that it indeed behaves as
$1/\omega^{4}$ for $\omega\gg k,\pi\,T$ region, which is in agreement with the
OPE expectations and satisfies the sum rule [26].

As discussed in the previous section, the ratio $G_{H}/(2\chi_{q}T)$ on the
lattice is free of renormalization. In order to compare it with the
perturbative correlator, we used $\chi_{q}=0.897\,T^{2}$ for the case of
$N_{f}=0$ and $\chi_{q}=0.842\,T^{2}$ for $N_{f}=2+1$. The value for $N_{f}=0$
was estimated using the Lüscher non-perturbative method [37], while for the
case of $N_{f}=2+1$, we used the perturbative estimate of order $g^{6}\log(g)$
[53]. The comparison between the correlator obtained from the perturbative
spectral function and the lattice calculation is illustrated in Figure 3.

In the left panel, we plotted the quenched case correlator up to a large value
of momentum, since the quenched results are continuum extrapolated. However,
for full QCD, the results are only provided for the first three lattice
momenta, as the lattice finite spacing effects are larger for higher momenta.
In the quenched case, we observe that our lattice correlator data
significantly exceeds the perturbative correlator at the lowest available
momentum. However, at higher momentum, the situation improves, although a
small non-perturbative difference remains. On the other hand, in full QCD,
non-perturbative effects are already prominent for all momentum values we
considered. The temperatures for the $N_{f}=0$ and $N_{f}=2+1$ cases are
$1.5\,T_{c}$ and $\sim 1.2\,T_{pc}$ respectively. The $N_{f}=0$ case
temperature in this “scaled temperature” unit is $25\%$ higher than the
$N_{f}=2+1$ case. Therefore, Figure 3 also illustrates that by going $25\%$
higher in temperature, one sees a significant suppression in the non-
perturbative effects.

##  V Spectral reconstruction

Figure 4: The Euclidean correlator in the T-L channel, illustrating the result
of fitting the polynomial ansatz in Equation 14 for
$\omega_{0}=\sqrt{k^{2}+(\pi T)^{2}}$ (red band) to the measured lattice data.
The quenched case is shown on the left and full QCD is shown on the right,
both as a function of $\tau$ for just a few indicative momenta.

In this section, we will discuss some approaches to spectral reconstruction
from the lattice Euclidean correlator. Despite the ill-conditioned nature of
the inverse problem, one can make use of additional, physically motivated,
information to improve the reliability of extracted spectral functions—a
strategy which has been successful in previous works [54, 55, 27, 25]. As
discussed in Section II, the known physical constraints on the spectral
function contain the sum rule Equation 7 and the asymptotic behavior
$1/\omega^{4}$ at $\omega\gg k,\pi T$. In the following, we discuss several
methods and models for the spectral reconstruction that fold in these physical
constraints [56].

###  V.1 Models of the spectral function

Firstly, we consider a model of the T-L spectral function in which the
$\omega$-dependence arises from connecting two regions (delineated by
$\omega_{0}\mathrel{\hbox to0.0pt{\lower 2.5pt\hbox{$\sim$}\hss}\raise
1.99997pt\hbox{$>$}}k,\pi\,T$). The ‘ultraviolet’ region
$\omega\geq\omega_{0}$ is modeled by $\rho_{>}$ which is constructed from
inverse even powers of $\omega$ (starting with $1/\omega^{4}$) in order to
satisfy the OPE result. The ‘infrared’ region $\omega\leq\omega_{0}$ is
modeled by $\rho_{<}$, for which we adopt the same polynomial put forward in
Ref. [55]. The two regions are matched continuously,

$\rho_{H}^{\rm
poly}(\omega)\;=\;\rho_{<}(\omega)\,\Theta(\omega_{0}-\omega)+\rho_{>}(\omega)\,\Theta(\omega-\omega_{0})\,,$
(14)

where $\Theta$ is the Heaviside step function. Because the spectral function
is expected to be smooth and differentiable across the light cone [51], we
introduce the parameters $\beta$ and $\gamma$ such that

$\displaystyle\rho_{<}(\omega_{0})$
$\displaystyle=\;\rho_{>}(\omega_{0})\;=\;\beta\,,$ (15)
$\displaystyle\rho^{\prime}_{<}(\omega_{0})$
$\displaystyle=\;\rho^{\prime}_{>}(\omega_{0})\;=\;\gamma\,.$ (16)

These considerations lead to the following model of the spectral function

$\displaystyle\rho_{<}$
$\displaystyle\equiv\frac{\beta\omega^{3}}{2\omega_{0}^{3}}\left(5-3\frac{\omega^{2}}{\omega_{0}^{2}}\right)-\frac{\gamma\omega^{3}}{2\omega_{0}^{2}}\left(1-\frac{\omega^{2}}{\omega_{0}^{2}}\right)+\frac{\delta\,\omega}{\omega_{0}}\left(1-\frac{\omega^{2}}{\omega_{0}^{2}}\right)^{2}\,,$
(17) $\displaystyle\rho_{>}$
$\displaystyle\equiv-\frac{\beta\,\omega_{0}^{4}}{7\,\omega^{4}}\left(54\frac{\omega_{0}^{4}}{\omega^{4}}\,-\,94\frac{\omega_{0}^{2}}{\omega^{2}}\,+\,33\right)+\frac{\gamma\,\omega_{0}^{5}}{140\,\omega^{4}}\left(-81\frac{\omega_{0}^{4}}{\omega^{4}}\,+\,92\frac{\omega_{0}^{2}}{\omega^{2}}\,-\,11\right)-\frac{16\,\delta\,\omega_{0}^{4}}{35\,\omega^{4}}\left(1-\frac{\omega_{0}^{2}}{\omega^{2}}\right)^{2}\,,$
(18)

Figure 5:  The spectral function $\rho_{H}$ for various momenta obtained from
the polynomial fitting (solid lines), and the dependence on the matching point
$\omega_{0}$ is also illustrated. The results for quenched QCD are shown on
the left, while the results for full QCD are shown on the right. The
perturbative version of the same spectral function is also shown for
$\mu=\mu_{\rm opt}$ (red points).

where the coefficients have been chosen to ensure the matching. The quantity
$\delta$ is another parameter which controls the slope of the spectral
function at $\omega=0$. The parameters $\beta$, $\gamma$, and $\delta$ are
chosen such that the sum rule is satisfied. The fitting parameters are
constrained to values that satisfy $\delta>0$, $\rho_{H}(k,k)>0$ and
$\frac{\partial G_{H}}{\partial\tau}\leq 0$. The resulting fits of the lattice
data with this model of the spectral function is shown in Figure 4 for
quenched and full QCD, while the spectral function obtained from this fit is
illustrated in Figure 5.

We estimate the error of the reconstruction by comparing different values of
$\omega_{0}=\sqrt{k^{2}+\nu\,(\pi T)^{2}}$ for $\nu=\\{1,2,5\\}$. In Figure 5
we show the variation of the spectral function for a these values of
$\omega_{0}$. At small momentum, we see that for both quenched and full QCD,
the variation with respect to $\omega_{0}$ becomes significant only at large
$\omega_{0}$. The polynomial spectral function also predicts a large negative
part at large $\omega$ compared to the perturbative spectral function. At
higher momentum, however, we see that the variation of the spectral function
with respect to $\omega_{0}$ becomes small for quenched QCD. The large
$\omega$ part of the spectral function is also approximately consistent with
the perturbative spectral function. At the same time, for full QCD, we see
that the variation with respect to $\omega_{0}$ is comparatively larger than
the quenched case, and the polynomial model spectral function shows a large
deviation from the corresponding perturbative spectral function. We pursue a
mock analysis of the perturbative correlator in Appendix A, showing that this
range of $\omega_{0}$ approximately captures the perturbative spectral
function.

The photon production rate, which is proportional to $D_{\text{eff}}$, can
then be obtained from the spectral function. The systematic error on $D_{\rm
eff}$ was obtained by combining the values of $\omega_{0}$ prescribed by
$\nu=\\{1,2,5\\}\,$. We provide the resulting numerical values for
$D_{\text{eff}}$ in Table 2 for quenched and in Table 3 for full QCD. The
results shown are from continuum extrapolated data set for quenched QCD. In
Appendix B, we discuss the dependence of the $D_{\text{eff}}$ on the lattice
spacing.

The second model of the T-L spectral function is a Padé-like ansatz that has
already been applied to the reconstruction problem in Ref. [25] and is given
by

$\displaystyle\rho^{\text{Pad\'{e}}}_{H}(\omega)=\frac{A\tanh(\frac{\omega}{2T})\big{(}1+B\omega^{2}\big{)}}{(\omega^{2}+a^{2})\big{[}(\omega+\omega_{0})^{2}+b^{2}\big{]}\big{[}(\omega-\omega_{0})^{2}+b^{2}\big{]}}\,.$
(19)

At small $\omega$ this spectral function reproduces the hydrodynamic
prediction $\rho_{H}(\omega)\sim\frac{A\omega}{\omega^{2}+a^{2}}$ for the IR
limit [23]. The remaining part of the spectral function is inspired by
qualities from AdS/CFT (see, e.g. Ref. [57]) and at large $\omega$ is
consistent with the OPE. In order to satisfy the sum rule, $B$ becomes a
function of $\omega_{0},\,a$ and $b$ (here the parameter $\omega_{0}$ is not
the same as the one from Equation 14). The fit then was performed with respect
to $A,\,a,\,b,\,\omega_{0}$ by minimizing the uncorrelated $\chi^{2}$. We
perform this fit in each of the bootstrap samples. The resulting distribution
of $D_{\text{eff}}T$ over the bootstrap samples is broad, indicating a near-
degenerate minimum $\chi^{2}$, as also identified in Ref. [27]. Therefore, we
quote the variation of the minimum and maximum $D_{\text{eff}}T$ of the
distribution within the error bars (rather than the standard deviation from
the mean values). The estimated $D_{\text{eff}}$ can be found in Table 2 for
quenched and Table 3 for full QCD. Sample spectral functions for the two
extreme cases are plotted in Figure 6, along with the perturbative spectral
functions. We see that the qualitative behavior of these spectral functions
are similar to the spectral function obtained from the polynomial ansatz.
However at large $\omega$ the Padé-like ansatz shows a large negative
contribution than the perturbative spectral function.

Figure 6:  Spectral function $\rho_{H}$ obtained from the Padé-like ansatz in
Equation 19, after fitting to the corresponding lattice Euclidean correlator.
On the left, the spectral functions are shown for quenched QCD at momentum
$k/T=4.71$, and on the right, for full QCD at momentum $k/T=4.19$, with the
minimum $D_{\text{eff}}T$ and maximum $D_{\text{eff}}T$ samples of the
bootstrap distribution.

###  V.2 Backus-Gilbert Estimate

Figure 7:  Backus-Gilbert estimated spectral function, with different colored
points representing different choices of $\lambda$, the regularizing
parameter. The same $\nu=1$ is used for both quenched (left) and full QCD
(right), corresponding to $\omega_{0}=\sqrt{k^{2}+(\pi T)^{2}}$ in the prior
function Equation 26.

Figure 8:  Dependence of the diffusion coefficient $D_{\text{eff}}T$, from the
BG-estimated spectral function, on the regularization parameter $\lambda$ for
$\nu=1$ is shown on the left and of the bootstrap average (indicated by black
horizontal lines). The same dependence on $\nu$ appears on the right. The top
two panels are for $N_{f}=0$ and the lower two panels are for $N_{f}=2+1\,$,
in each case several momenta are shown.

Figure 9:  Dependence of the resolution function
$\delta(\omega,\omega^{*}=|\vec{k}|)$ on $\omega\,$, for parameter values of
$\nu=\\{0,8\\}$ and at $\lambda=0.7\,$. Both panels are for quenched QCD, the
left for $k/T=1.57$ and the right for $k/T=4.71\,$.

The Backus-Gilbert (BG) [58] method is a technique commonly used for spectral
reconstruction in the context of QCD. It involves computing a smeared
estimator for the spectral function, denoted as $\rho^{\rm
BG}(\omega,\vec{k})$. Consequently, it represents the actual spectral function
$\rho(\omega,\vec{k})$ convoluted with a so-called resolution function
$\delta(\omega,\omega^{*})$, namely

$\rho^{\rm BG}(\omega,\vec{k})=\int{\rm
d}\omega^{*}\,\delta(\omega,\omega^{*})\,\rho(\omega^{*},\vec{k})\,.$ (20)

When the resolution function $\delta(\omega,\omega^{*})$ is a Dirac delta
function $\delta(\omega-\omega^{*})$, the Backus-Gilbert estimate provides an
exact reconstruction of the spectral function. However, due to the problem at
hand being ill-conditioned, this is not possible in practice. Therefore, the
objective of this method is to minimize the width of the resolution function
in order to achieve better accuracy.

In order to identify the resolution function in this method, the correlator in
Equation 3 is rewritten as

$G_{H}(\tau,\vec{k})\;=\;\int_{0}^{\infty}{\rm
d}\omega\,\frac{\rho(\omega,\vec{k})}{f(\omega,\vec{k})}\,\tilde{K}(\omega,\tau)\,,$
(21)

where
$\tilde{K}(\omega,\tau)\equiv\frac{f(\omega,\vec{k})\cosh[\omega(1/2T-\tau)]}{\sinh(\omega\beta/2)\pi}$.
Here, some arbitrary function $f(\omega,\vec{k})$ is introduced by hand. It
can be chosen to incorporate the known physics information about the spectral
function asymptotics and should also remove the singularity at $\omega=0$
originating from the kernel.

The BG-estimated spectral function is expressed as a linear combination of the
measured lattice correlator points

$\frac{\rho_{H}^{\rm
BG}(\omega,\vec{k})}{f(\omega,\vec{k})}\;=\;\sum_{i}q_{i}(\omega,\vec{k})\,G_{H}(\tau_{i},\vec{k})\\\
\;=\;\int_{0}^{\infty}{\rm
d}\omega^{*}\,\delta(\omega^{*},\omega)\,\frac{\rho_{H}(\omega^{*},\vec{k})}{f(\omega^{*},\vec{k})}\,,$
(22)

where we identify
$\delta(\omega,\omega^{*})=\sum_{i}q_{i}(\omega,\vec{k})\tilde{K}(\omega^{*},\tau_{i})$.

For any $\omega$, the resolution function is therefore fully characterized by
the coefficients $q_{i}$. One way of constructing these coefficients, is by
minimizing $A(\omega)$ (w.r.t. $q_{i}$), a measure of the width of the
resolution function

$A(\omega)\;\equiv\;\int{\rm
d}\omega^{*}\,(\omega-\omega^{*})^{2}\delta(\omega,\omega^{*})^{2}\;=\;{\bm{q}^{\intercal}(\omega)\bm{W}(\omega)\bm{q}(\omega)}\,,$
(23)

where $[\bm{W}]_{ij}(\omega)\equiv\int_{0}^{\infty}{\rm
d}\omega^{*}\,(\omega-\omega^{*})^{2}\tilde{K}(\omega^{*},\tau_{i})\tilde{K}(\omega^{*},\tau_{j})$.

In the BG method, the minimization is performed under the constraint that the
integral of the resolution function $\delta(\omega,\omega^{*})$ over $\omega$
is equal to 1. With this constraint, the solution for the coefficients
$[\bm{q}]_{i}$ is given by

$\bm{q}(\omega)\;=\;\frac{\bm{W}^{-1}(\omega)\,\bm{r}}{\bm{r}^{\intercal}\,\bm{W}^{-1}(\omega)\,\bm{r}}\,,$
(24)

where $[\bm{r}]_{i}\equiv\int{\rm d}\omega\,\tilde{K}(\omega,\tau_{i})\,$. In
reality, the matrix $\bm{W}$ is close to singular—due to its large condition
number—and therefore requires regularization.

One of the popular regularization schemes is to take

$\bm{W}^{\lambda}\;\equiv\;\lambda\,\bm{W}+(1-\lambda)\,\bm{S}\,,$ (25)

where $\bm{S}$ is the covariance matrix of the lattice correlator. Other
regularization schemes also exist in literature, e.g. the Tikhonov
regularization [59] where $\bm{S}$ is replaced by the identity matrix. In this
paper, we have used a diagonal covariance matrix for $\bm{S}$.

The chosen regularization introduces an additional dependence of the spectral
function on the parameter $\lambda$. The spectral function also depends on the
prior function. We choose the form of the prior function as

$f(\omega,\vec{k})=\left(\frac{\omega_{0}}{\omega}\right)^{4}\tanh\left(\frac{\omega}{\omega_{0}}\right)^{5}\,,$
(26)

where $\omega_{0}=\sqrt{k^{2}+\nu(\pi T)^{2}}$ as before. In this case, we
will vary $\nu=\\{0,\frac{1}{2},1,2,5,8\\}\,$. The functional form is based on
the information that we already have about the spectral function, i.e it
behaves as $1/\omega^{4}$ for large $\omega\gg k$ values consistent with the
OPE expectation, and linearly with $\omega$ for small $\omega$ values.

The resulting BG-estimated spectral function, and thus the photon production
rate, depend on the parameters $\nu$ and $\lambda$. The BG-estimated spectral
functions for the momentum $k/T=1.57$ (for $N_{f}=0$) and $k/T=2.09$ (for
$N_{f}=2+1$) are shown in Figure 7. In this figure, we showed the dependence
of spectral function on $\lambda\in[0.1,0.8]$ for $\nu=1\,$. (The errors on
the spectral functions are statistical.) Here we also observed the general
feature, that the UV part of the spectral function is suppressed compared to
the IR part. The photon production rate is proportional to the diffusion
coefficient $D_{\text{eff}}$, defined from $\rho_{H}$ in Equation 9, which is
plotted on the left panel of Figure 8 as a function of $\lambda$ for various
momenta, but for a given value of $\nu$. We observe that for both $N_{f}=0$
and $N_{f}=2+1$ case the estimated $D_{\text{eff}}$ is stable with respect to
the variation of the $\lambda$ parameter, although the error starts growing
larger as $\lambda\to 1\,$. The values of $D_{\text{eff}}$ for various
$\lambda$ are combined using a bootstrap average over the plateau in region
$\lambda=[0.1,0.7]$ in steps of $\Delta\lambda=0.1$. This $\lambda$-averaged
version of $D_{\text{eff}}T$ is plotted as a function of $\nu$ on the right
panel of Figure 8. For both $N_{f}=0$ and $N_{f}=2+1\,$, the coefficient
$D_{\text{eff}}T$ exhibits small variations with respect to $\nu\,$, except
for the smallest cases displayed, where deviations emerge for
$\nu\mathrel{\hbox to0.0pt{\lower 2.5pt\hbox{$\sim$}\hss}\raise
1.99997pt\hbox{$<$}}2\,$.

The reason for this can be understood by looking at the resolution function,
which is depicted in Figure 9. In this plot, we show the resolution function
for the momentum $k/T=1.57$ and $k/T=4.71$ for the minimum and maximum $\nu$
values in the case of $N_{f}=0$ and $\lambda=0.7\,$. We find that for
$k/T=4.71$, the resolution function is peaked around the light cone and
exhibits variations of $\sim 9\%$ from $\nu=0$ to $\nu=8$. As a result, the
BG-estimated spectral function in this case gets most of its contribution from
the light cone region. Apart from the lowest momentum, this pattern is
repeated for the other momenta studied here. For $k/T=1.57$ and $\nu=0$, we
find that the resolution function is peaked at the light cone. However, for
$k/T=1.57$ and $\nu=8$, the resolution function instead shows a peak at
$\omega=0\,$. Therefore, the BG-estimated photon production rate from these
values of $\nu$ gets a significant contribution from the $\omega=0$ region.
Nevertheless, when quoting the value of the photon production rate, we
included these variations with respect to $\nu$ as a source of systematic
error. In the case of full QCD, we find the results for the momentum
$k/T=6.28$ are negative up to $\nu=2$. The tabulated values of
$D_{\text{eff}}T$, determined from the BG method, can be found in Table 2 for
quenched and Table 3 for full QCD. This implies a lower bound for
$D_{\text{eff}}$ equal to zero for the photon production rate because a
negative rate would be unphysical.

###  V.3 Gaussian Process Regression

Figure 10:  Spectral function from GPR: Quenched on the left, and full QCD on
the right. The error band includes both the error propagation of the lattice
data and the systematic error from the GP model.

Gaussian Processes Regression (GPR), also known as ‘Kriging’ or Wiener-
Kolmogorov prediction, is a widely employed probabilistic interpolation
method, especially well-suited for noisy and irregularly distributed data.
Originating from early applications in geostatistics, today Gaussian Processes
(GPs) have a wide range of applications, primarily in modeling and
classification problems; see Ref. [60] for a review and Ref. [61] for a
textbook introduction. Because no explicit functional basis is required, GPR
offers a very general approach to interpolating noisy data. Instead, the
characteristics of the underlying function are defined implicitly through the
so-called kernel of the GP, making GPR a very flexible non-parametric model
when considering interpolations of data with an unknown functional basis.
Here, we will summarize the main concepts and features of the inversion
procedure with GPs, while we refer to Appendix C for the technical discussion
and Refs. [62, 63] for additional details on the method.

In essence, GPs can be viewed as a way of defining a Gaussian distribution
over functions, with some prior mean and covariance. When describing some
function as a GP, we can then make prediction about function values by
conditioning the GP on the available data. This is a standard procedure in
Bayesian statistics and leads to the posterior predictive distribution. Due to
the probability density being inherently Gaussian, this posterior probability
can be computed analytically, and the resulting posterior is again a
multivariate normal distribution with a posterior mean and covariance, cf.
Appendix C. This predictive distribution is a distribution of all possible
functions fitting the given data and we use the mean of this distribution in
order to make predictions.

The important property that makes GPs practicable for the inversion of
integral transformations in particular, is the fact that linear
transformations preserve Gaussian statistics. This means that we can condition
the predictive distribution of the GP not only on direct data, but also, in
the case of the spectral function, on the correlator data, since the spectral
function is connected to the correlator by a linear transformation, cf.
Equation 3. The resulting predictive distribution is then a distribution over
the family of functions which fulfill the integral constraint given by the
correlator data. This method has been proposed in Ref. [62] and applied to a
number of numerical reconstructions in the context of QCD, see e.g. Refs. [63,
64, 65].

GPR has two main appeals, when considering inverse problems. Firstly, we do
not choose a finite basis for the interpolation. Instead, we generally focus
on GP models that perform the interpolation in an infinite dimensional
functional basis. This means that we can model every continuous function and
do not have to assume a prior functional basis. Secondly, the inclusion of
additional data, either direct knowledge of the spectral function or other
linearly connected constraints such as the sum rule, e.g. Equation 7, is
straightforward and can further restrict the range of possible functions.
Otherwise, when only considering the correlator data, the uncertainty in the
spectral function prediction can be rather large, as expected from an ill-
conditioned inversion.

As mentioned before, the GP prior covariance is fully characterized by the so-
called GP-kernel. The GP-kernel constitutes the prior assumption about the
function that is interpolated. Typically, these assumptions encompass very
general features, such as continuity and a general length scale. A popular
choice is the radial basis function (RBF) kernel Equation 30. As a so-called
universal GP-kernel, the associated GP is a universal approximator, meaning
that this choice does not restrict the possible spectral functions [66, 67].
However, when having additional information about the functional structure of
the spectral function, such as UV asymptotics, the GP-kernel can be extended
to incorporate this information in certain regions, as described in Appendix
C.

####  V.3.1 GPR Reconstruction Details

Figure 11:  Effective diffusion coefficent $D_{\text{eff}}$ obtained, as in
Equation 9, from the GPR reconstruction of $\rho_{H}$ for quenched on the left
and for full QCD on the right. The error band includes both the error
propagation of the lattice data and the systematic error from the GP model.

Since we have data for several spatial momenta, the reconstruction of the
spectral function is not only performed in the $\omega$, but also
simultaneously in the $k$-direction. This ensures continuity in both
directions and consequently increases the stability of the reconstruction.
However, apart from continuity, we do not assume any additional structure of
the spectral function in the momentum direction. For the reconstruction, we
include the lattice data for all the available spatial momenta as well as the
sum rule Equation 7.

As a very common GP-kernel choice, the RBF kernel Equation 30 is used. In
general, we observe that choosing different smooth, stationary and universal
GP-kernels have very little impact on the reconstruction. In order to include
the known UV asymptotics, $\rho_{H}^{\text{\tiny{UV}}}(\omega)\propto
1/\omega^{4}$, the GP-kernel is modified in the UV regime to restrict the
functional basis to the asymptotic behavior resulting in Equation 31. The
transition from the RBF kernel to the UV kernel is controlled by smooth step
functions, for more details on the exact form of the resulting GP-kernel see
Appendix C. Combining the RBF kernel in two dimension with the asymptotic
kernel results in a total number of five GP-kernel parameters:
$l_{\omega},\,l_{k},\,\sigma,\,l_{\text{\tiny{UV}}},\,\nu$. The lenghtscale in
$k$ and $\omega$ direction is controlled by $l_{k}$ and $l_{\omega}$
respectively, while $\sigma$ gives a prior estimate on the variance of the GP.
The position of the transition from the RBF kernel towards the UV kernel is
controlled by $\omega_{0}=\sqrt{k^{2}+\nu(\pi T)^{2}}$, where $\nu$ is varied,
in analogy to the polynomial fit ansatz, and $l_{\text{\tiny{UV}}}$ is the
smoothness of the transition between the two kernels.

These parameters are optimized by minimizing the associated negative log-
likelihood (NLL) Section C.1. Since this optimization is performed in a five
dimensional parameter space and the parameters are not fully independent, this
optimization generally does not converge consistently. Additionally, open
directions in the parameter space, e.g. towards vanishing length scales are
possible. This can be viewed as a manifestation of the inherent ill-
conditioning of this inverse problem. Therefore, and in order to capture the
systematic error from the introduction of the UV kernel, we perform a
Metropolis-Hastings sampling of the NLL in the parameter space, see Section
C.1 for details. While the error on the lattice data is propagated through the
GP and is given by the covariance of the posterior predictive distribution of
the GP, the systematic error is estimated by the uncertainty in the parameters
while minimizing the associated NLL. This systematic error is mostly the
result of introducing an explicit functional basis in the UV, while a change
in the RBF kernel parameters results in a much smaller deviation in the
spectral function. The final error estimate on the spectral function from the
GP reconstruction is given by a combination of the statistical and systematic
error.

The resulting spectral functions corresponding to the quenched and full QCD
lattice correlator are presented in Figure 10. In order to compare with the
reconstructions from the fits and the BG method, we have evaluated the
spectral function at the spatial momenta $k$, where the lattice data is
available. In general however, the GP gives a statistical estimate for all
momenta in between the available data. This gives a continuous estimate on the
thermal photon rate not only at the given spatial momentum values, but also in
between, see Figure 11. The validity of the interpolation between data points
has been additionally confirmed by reconstructing the effective diffusion
coefficient while disregarding lattice data for a single spatial momentum.
This gives, as expected, very similar results with an increased error around
this momentum value. Although the interpolation in momentum-direction give
reliable estimates, the extrapolation towards smaller or higher momenta comes
with large uncertainties. Ultimately, when extrapolating towards $k=0$, we
recover the GP prior as the current assumptions about the spectral function do
not allow for a systematic extrapolation. The resulting values for
$D_{\text{eff}}T$ from GP reconstruction are given in Table 2 for quenched and
in Table 3 for full QCD.

##  VI Comparison of $D_{\text{eff}}T$

Figure 12:  Effective diffusion coefficient from all spectral reconstruction
methods studied in this paper. The $N_{f}=0$ case appears on the left panel,
$N_{f}=2+1$ on the right. Here we have converted the results in physical
units. For better visibility, error bars at the same momentum value are
slightly shifted on the horizontal axis.

In the previous sections we discussed various models and methods for the
spectral reconstruction from lattice data. Now, we compare the photon
production rate obtained from these different methods.

In Figure 12, we show the effective diffusion coefficient as a function of
momentum for all the reconstruction methods we performed. We observe that the
uncertainty in the estimation of $D_{\text{eff}}$ is larger at small momenta.
This is expected because at small momenta the lattice correlators are flat
compared to the larger momenta, corresponding to a narrow peak in
$\rho_{H}/\omega$ at $\omega=0$ for small momenta. On the other hand, at
larger momenta, $D_{\text{eff}}$ shows better agreement, as different methods
provide closer estimates, while the relative error on the reconstructions
remains approximately constant.

On the same plot, we also display the corresponding perturbative estimate of
$D_{\text{eff}}(k)$. The convergence of the perturbative estimate can be
checked by varying the renormalization scale $\mu$. Here, we estimate the
error band by varying the scale from333 The smaller value of $\mu$ corresponds
to an upper estimate of $D_{\rm eff}\,$, due to the larger QCD coupling being
probed. In our case, the extreme choice of $\mu=0.5\,\mu_{\text{\tiny{opt}}}$
leads to $\alpha_{s}\simeq 0.5$ which is where most of the sensitivity to
scale variation originates.  $\mu=0.5\,\mu_{\text{\tiny{opt}}}$ to
$\mu=2\,\mu_{\text{\tiny{opt}}}$ (where the ‘optimal’ scale is defined in
Section IV). The error band on the perturbative estimate is slightly smaller
for quenched QCD at $T=$470\text{\,}\mathrm{MeV}$$ compared to full QCD at
$T=$220\text{\,}\mathrm{MeV}$$, which may indicate that the perturbative
framework is better behaved in this case.

In the case of quenched QCD, our non-perturbative estimate of $D_{\text{eff}}$
from the Padé-like ansatz is consistent with the perturbative estimate within
the error bar. However, at the lowest momentum, the Polynomial-ansatz and GP
regression favor a higher photon rate compared to the perturbative estimate,
whereas at large momentum, the polynomial-ansatz, BG, and GP estimates give a
lower photon rate than the perturbative estimate. When it comes to full QCD,
we observe in all the non-perturbative estimates a significantly more rapid
decrease in the photon production rate compared to the perturbative estimate .

This aligns with the fact that, in units of the respective transition
temperature, quenched QCD is at approximately $25\%$ higher temperature than
full QCD. As a result, one would expect that the quenched QCD result will be
closer to the perturbative result than the full QCD case.

It is important to note that the perturbative calculation of the spectral
function on the light cone relies on a coupling controlled entirely by the
temperature. Consequently, one should expect the non-perturbative estimate to
approach the perturbative estimate only at high temperatures and for momenta
$k\sim T$.

We combined the $D_{\text{eff}}$ from different models to obtain an overall
error bar arising from model estimates. This is done by spanning the highest
and lowest value of $D_{\text{eff}}$ from models for a given momentum $k$.

$k/T$ | polynomial | Padé | Backus-Gilbert | GP regression
---|---|---|---|---
$1.57$ |  $0.549\pm 0.035$ |  $0.469\pm 0.368$ |  $0.302\pm 0.175$ | $0.616\pm 0.168$
$3.14$ |  $0.161\pm 0.092$ |  $0.129\pm 0.123$ |  $0.119\pm 0.057$ | $0.138\pm 0.046$
$4.71$ |  $0.087\pm 0.048$ |  $0.089\pm 0.067$ |  $0.054\pm 0.031$ | $0.042\pm 0.021$
$6.28$ |  $0.054\pm 0.018$ |  $0.058\pm 0.044$ |  $0.022\pm 0.018$ | $0.039\pm 0.015$
$7.85$ |  $0.041\pm 0.007$ |  $0.066\pm 0.066$ |  $0.016\pm 0.016$ | $0.029\pm 0.013$
$9.42$ |  $0.031\pm 0.007$ |  $0.074\pm 0.045$ |  $0.013\pm 0.013$ | $0.019\pm 0.010$
Table 2:  Extracted values of $D_{\text{eff}}T$ from different methods for $N_{f}=0\,$.  $k/T$ | polynomial | Padé | Backus-Gilbert | GP regression
---|---|---|---|---
$2.09$ |  $0.369\pm 0.040$ |  $0.206\pm 0.085$ |  $0.197\pm 0.127$ | $0.412\pm 0.179$
$4.19$ |  $0.102\pm 0.034$ |  $0.062\pm 0.058$ |  $0.063\pm 0.047$ | $0.046\pm 0.047$
$6.28$ |  $0.031\pm 0.013$ |  $0.010\pm 0.010$ |  $-0.011\pm 0.047$ | $0.014\pm 0.022$
Table 3:  Extracted values of $D_{\text{eff}}T$ from different methods for
$N_{f}=2+1\,$.

##  VII Conclusion

We estimated the thermal photon production rate from the Euclidean lattice
correlators in quenched QCD ($T=$470\text{\,}\mathrm{MeV}$$) and full QCD
($T=$220\text{\,}\mathrm{MeV}$$) with an unphysical pion mass of
$320\text{\,}\mathrm{MeV}$. Motivated by Ref. [25], we have utilized the T-L
correlator for calculating the photon rate. This choice offers several
advantages over the vector correlator, primarily due to the absence of a
significant UV component in the T-L correlator. The spectral function of this
correlator satisfies a sum rule along with a known asymptotic form in the UV.
In the quenched case this T-L correlator has been extrapolated to the
continuum whereas in full QCD the results are from a single lattice spacing.
We studied non-perturbative effects in this correlator by comparing it with
its perturbative estimate, shown in Figure 3. Figure 2 displays the
perturbative spectrum, which was calculated at NLO away from the light cone,
while leading-order LPM resummation has been done near the light cone [38].

We observe that in quenched QCD, the non-perturbative effects are much more
pronounced at lower momenta compared to higher momenta. This is evident as the
lattice correlators are significantly flatter in contrast to the perturbative
correlator. As we move towards larger momenta, the lattice data also begins to
exhibit curvature similar to the perturbative correlator. However, it appears
that non-perturbative effects persist even at the highest available momentum
we have.

In the case of full QCD, the lattice data displays a substantially flatter
behavior (as a function of $\tau$) compared to the perturbative correlator,
indicating significant non-perturbative modification.

To obtain photon rate, we performed the spectral reconstruction using the
lattice correlators we determined. Since the inverse problem is an ill-
conditioned problem, we used multiple models and methods for the spectral
reconstruction. In all models, we constrained the spectral function according
to the sum rule and incorporated the known asymptotic form. The first method
we used is a simple polynomial ansatz of the spectral function in the IR
region, which is based on the smoothness of the spectral function across the
light cone. We smoothly connected the IR part with a UV part of the spectral
function consistent with the OPE expansion, as given in Equation 14 . This
spectral function also satisfies the sum rule. Using this form we performed
the fitting with the lattice data with respect to the free parameters of the
model. The resulting spectral function is plotted in Figure 5 along with the
perturbative spectral function. We see that at small momentum, the spectral
function has a larger dependence on the matching point $\omega_{0}$ compared
to the case of large momentum. As the correlator is very flat at small
momentum, it signifies a very narrow peak at the origin, which is difficult to
determine.444The uncertainty with respect to the variation of $\omega_{0}$ is
taken to be the systematic uncertainty on our result. Here we also observed
that except the smallest momentum, for the quenched case the resulting non-
perturbative spectral function is much closer to the perturbative spectral
functions compared to the full QCD case. This is expected because effectively
the temperature for the quenched case is larger than the full QCD case. As a
second model we used a Padé-like ansatz of the spectral function as put
forward in [25]. A sample spectral function is shown in Figure 6.

We also used the well-known BG method, to obtain a smeared spectral function,
which is a convolution of the actual spectral function with a resolution
function. The additional prior function we used is given in Equation 26,
whereby it contains the known asymptotic behavior of the spectral function.
This method has the stability parameter $\lambda$ and the parameter $\nu$ from
the prior function. We found that the photon production rate is stable w.r.t.
the variation of $\lambda$ and $\nu$ for all the momenta. For the lowest
available momentum we found that although the result is stable with respect to
$\lambda$, it shows some dependence on $\nu$, while the photon production rate
nevertheless saturates at large $\nu$. This variation with respect to $\nu$
has been taken as a systematic error on the photon production rate.

As a third method we applied Gaussian Process regression in order to obtain a
prediction for the spectral function. Instead of giving a singular prediction,
this method returns a distribution over functions that agree with the given
correlator data. The T-L correlator data is reconstructed continuously in both
$\omega$ and $k$ directions, giving a continuous prediction for the effective
diffusion coefficient, see Figure 11. By considering the UV asymptotic
behavior from the OPE and the sum rule, we can reduce the variance on the GP
posterior in order to make meaningful predictions.

The photon production rate in terms of $D_{\text{eff}}$ is plotted in Figure
12. The error bar from the models have been combined to provide a single error
bar. As momentum increases, the results from different methods begin to
converge within their systematic error bars. In the quenched case, we observe
that the lattice estimate of the photon production rate approaches the
perturbative estimate. In the full QCD case, the photon production rate
becomes smaller than the perturbative estimate and becomes consistent with
zero at a momentum of approximately
$0.8\text{\,}\mathrm{GeV}1.6\text{\,}\mathrm{GeV}$.

There are many ways one can proceed for future studies. Currently, the full
QCD results are obtained at finite lattice spacing. For larger momentum, the
cutoff effects could be significant, and a continuum extrapolation is needed.
The full QCD results in this paper are based on unphysical pion mass
configurations. To make predictions aligned with experimental results, it is
important to calculate the photon rate using physical pion mass gauge
configurations. Although it is expected that at high temperatures, the effect
of quark mass on the photon rate should be small. It is also important to
calculate the photon rate for different temperatures, which will be relevant
as input to hydrodynamics.

##  VIII Acknowledgments

We would like to thank Luis Altenkort and Hai-Tao Shu for their work on
integrating the mesonic correlator into the QUDA code and generating gauge
field configurations using SIMULATeQCD [68]. D. B. would like to thank Guy D.
Moore for various discussions on the photon production rate. D. B. would also
like to appreciate the discussion with Harvey B. Meyer during the Lattice 2022
conference for the discussion on sum rules and asymptotic behavior of the
spectral function. S. A., D. B., O. K., J. T., T. U. and N. W. are supported
by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -
Project number 315477589-TRR 211. N. W. acknowledges support by the State of
Hesse within the Research Cluster ELEMENTS (Project ID 500/10.006). For the
computational work We used the Bielefeld GPU cluster and the JUWELS
supercomputer. The authors gratefully acknowledge the Gauss Centre for
Supercomputing e.V. for funding this project by providing computing time
through the John von Neumann Institute for Computing (NIC) on the GCS
Supercomputer JUWELS at Jülich Supercomputing Centre (JSC). G. J. was funded
by the U.S. Department of Energy (DOE), under grant No. DE-FG02-00ER41132, and
now by the ANR under grant No. ANR-22-CE31-0018 (AUTOTHERM). A. F.
acknowledges support by the National Science and Technology Council of Taiwan
under grant 111-2112-M-A49-018-MY2.

##  Appendix A Mock-analysis of perturbative data with Polynomial ansatz

Figure 13:  Mock analysis of the data generated from the perturbative
correlator. In the left panel, the red band represents the fitting of the data
points with the polynomial ansatz, while the black continuous line represents
the $\text{NLO}+\text{LPM}^{\text{LO}}$ perturbative correlator. In the right
panel, the reconstructed spectral function is displayed along with the exact
$\text{NLO}+\text{LPM}^{\text{LO}}$ spectral function.

In this section, we conduct a mock analysis to demonstrate the effectiveness
of the polynomial ansatz fitting. To generate the mock data set, we use the
perturbative $\text{NLO}+\text{LPM}^{\text{LO}}$ correlator, which is obtained
from the known perturbative spectral function $N_{f}=0$ case. In the mock
data, we include an error bar equivalent to the corresponding lattice
correlator. We used 11 equally spaced data points in the interval from $\tau
T\simeq 0.17$ to $\tau T=0.5\,$, matching the number of continuum-extrapolated
lattice data points for $N_{f}=0$.

We perform the fitting with the polynomial ansatz for the spectral function
and obtain the fit, which is shown in the left panel of Figure 13 for the
momentum $k/T=3.14$ and $k/T=6.28$. The black line corresponds to the exact
correlator, while the red band corresponds to the fitting of the spectral
function for $\omega_{0}=\sqrt{k^{2}+(\pi T)^{2}}$.

The right panel of Figure 13 shows the resulting spectral function for various
values of $\omega_{0}\,$ at $k/T=4.71$. We can see that this ansatz is indeed
able to predict the exact perturbative spectral function within the range of
$\omega_{0}=\sqrt{k^{2}+\nu(\pi T)^{2}}$ where $\nu=1\,...\,5\,$.

##  Appendix B $D_{\text{eff}}$ from finite lattice spacing for $N_{f}=0$

In this section, we will compare the estimate of the effective diffusion
coefficient from different lattice spacings in order to estimate the cutoff
effects on the reconstruction. This will be illustrated for the polynomial
ansatz of the spectral function and the GPR reconstruction. A comparison of
the different $D_{\text{eff}}$ values from the polynomial ansatz can be found
in Figure 14, while the same comparison for the GPR reconstruction can be
found in Figure 15. We observe that for both methods the results for small
spatial momenta agree remarkably well, as different values remains consistent
in the margin of error. For the polynomial ansatz and spatial momenta
$k/T>4.71$ we observe mild cutoff effects on $D_{\text{eff}}$, while in the GP
reconstruction, the error generally grows for higher momenta. This outcome
aligns with our expectations, considering that cutoff effects are already
minimal in the correlator.

Figure 14:  Effective diffusion coefficient, $D_{\text{eff}}$ obtained from
the polynomial estimate for different lattice spacing along with continuum
extrapolated value. Here $N_{f}=0\,$.

Figure 15:  Effective diffusion coefficient, $D_{\text{eff}}$ obtained from
the GP regression model for different lattice spacing along with continuum
extrapolated value. Here $N_{f}=0\,$.

##  Appendix C Inversion with GPR

Figure 16: Metropolis-Hastings sampled GP kernel parameters in the negative
log-likelihood plane.

Here, we will give a short introduction to the inversion of spectral functions
with GPR, while we refer to Refs. [62, 63] for a more extensive introduction.
To start, let us assume that the knowledge of the spectral function before
making observations of the correlator can be described by a GP prior. This is
written as

$\displaystyle\rho(\omega)\sim\mathcal{GP}\big{(}\mu(\omega),C(\omega,\omega)\big{)}\,,$
(27)

where $\mu$ and $C$ denote the mean and the covariance of the GP prior. The
prior mean can be generally set to zero, as its contribution can be fully
absorbed into the covariance.

The conditional posterior distribution for $\rho(\omega)$ given observations
of the propagator $G_{H,i}$ at $N_{G}$ discrete Euclidean times
$\tau_{i}\equiv[\bm{\tau}]_{i}$ can be derived in closed form,

$\displaystyle\rho(\omega)|G(\bm{\tau})\sim$
$\displaystyle\,\mathcal{GP}\Bigl{(}\bm{w}^{\intercal}(\omega)(\bm{W}+\sigma_{n}^{2}1)^{-1}G_{H}(\bm{\tau}),$
$\displaystyle
k(\omega,\omega)-\bm{w}^{\intercal}(\omega)(\bm{W}+\sigma_{n}^{2}1)^{-1}\bm{w}(\omega)\Bigr{)}\,,$
(28)

where

$\displaystyle[\bm{w}]_{i}(\omega)=\int\text{d}\omega^{\prime}\,K(\tau_{i},\omega^{\prime})C(\omega^{\prime},\omega)\,,$
$\displaystyle[\bm{W}]_{ij}=\int\text{d}\omega^{\prime}\text{d}\omega^{\prime\prime}\,K(\tau_{i},\omega^{\prime})K(\tau_{j},\omega^{\prime\prime})C(\omega^{\prime},\omega^{\prime\prime})\,.$
(29)

and
$K(\tau,\omega)\equiv\frac{\cosh(\omega(1/2T-\tau))}{\pi\sinh(\omega/2T)}$.
Appendix C is a standard result in multivariate statistics and is essentially
equivalent to GPR with direct observation only with the addition of the
integral transformation, that is supposed to be inverted. This posterior
encodes the knowledge of the spectral function under the constraint of the
correlator observation and directly accounts for the error estimations on the
correlator, here denoted as $\sigma_{n}$.

The prior covariance of the GP is characterized by the so-called kernel. This
kernel is generally a function with a small number of hyperparameters and
fully characterizes the GP. The kernel controls very general features of the
interpolation, such as differentiability of the function, a characteristic
lengthscale or an estimate on the variance of the interpolation. A commonly
used kernel choice is the radial basis function (RBF) or Gaussian kernel,
defined as

$\displaystyle
C_{\text{\tiny{RBF}}}(\omega,\omega^{\prime})=\sigma^{2}\exp\left(-\frac{(\omega-\omega^{\prime})^{2}}{2l^{2}}\right)\,,$
(30)

where $\sigma$ controls the overall magnitude and $l$ is a generic
lengthscale. The RBF kernel falls under the class of so-called universal
kernels [66, 67]. Such kernels have the universal approximation property; any
continuous function can be approximated by the corresponding GP.

The kernel can be extended in order to capture knowledge about the asymptotics
of the spectral function, as proposed in Ref. [65]. The main idea is to reduce
the basis function of the GP kernel in a certain region to just represent the
known asymptotics of the spectral function. For the known UV asymptotics
$\rho_{H}^{\text{\tiny{UV}}}(\omega)\propto 1/\omega^{4}$ of the spectral
function, this can be achieved by taking the asymptotic kernel

$\displaystyle
C_{\text{\tiny{UV}}}(\omega,\omega^{\prime})\;=\;\rho_{H}^{\text{\tiny{UV}}}(\omega)\rho_{H}^{\text{\tiny{UV}}}(\omega^{\prime})\;=\;\big{(}\omega\,\omega^{\prime}\big{)}^{-4}\,.$
(31)

In order to have a smooth transition between the universal and the asymptotic
kernel, the transition is controlled by a soft step function

$\displaystyle\theta^{\pm}(\omega;l_{\text{\tiny{UV}}},\omega_{0})=\frac{1}{1+\exp(\pm(\omega-\omega_{0})/l_{\text{\tiny{UV}}}}\,.$
(32)

This results in the full kernel

$\displaystyle C(\omega,\omega^{\prime})$
$\displaystyle=\theta^{-}(\omega)\theta^{-}(\omega^{\prime})C_{\text{\tiny{RBF}}}(\omega,\omega^{\prime})$
$\displaystyle+\theta^{+}(\omega)\theta^{+}(\omega^{\prime})C_{\text{\tiny{UV}}}(\omega,\omega^{\prime})\,.$
(33)

The parameters of the soft step function, i.e. the position and the
lenghtscale of the transition, and of the RBF kernel are subject to likelihood
optimization, described in the following.

###  C.1 GPR parameter optimization

The GP kernel parameter optimization is a pivotal part of the reconstruction,
since the kernel parameters control the implicit structure of the
reconstructed function. In order to optimize these parameters, the negative
log-likelihood (NLL) of the associated GP is minimized. The NLL of a general
GP, with matrices as defined in Appendix C, is given by

$\displaystyle-\log
p(G_{H}(\bm{\tau})|\bm{\sigma})=\frac{1}{2}G_{H}(\bm{\tau})^{\intercal}\left(\bm{W_{\sigma}}+\sigma_{n}^{2}1\right)^{-1}G_{H}(\bm{\tau})$
$\displaystyle+\frac{1}{2}\log\det(\bm{W_{\sigma}}+\sigma_{n}^{2}1)+\frac{N}{2}\log
2\pi\,,$ (34)

where the dependence of the kernel on the hyperparameters is indicated by
$\bm{\sigma}$. These parameters are optimized by performing a Metropolis-
Hastings scan of the NLL in the five dimensional parameters space spanned by
$l_{\omega},\,l_{k},\,\sigma,\,l_{\text{\tiny{UV}}},\,\nu$. When optimizing
the parameters for the one dimensional GP, e.g. for one fixed spatial
momentum, we notice that the overall magnitude $\sigma$ of the kernel
decreases with $1/k^{2}$, while the other parameters fluctuate, but do not
have a qualitatively different behavior. We therefore rescale the magnitude
parameter for the two-dimensional reconstruction as
$\sigma\rightarrow\sigma/k^{2}$. In Figure 16, the parameters scan of the NLL
is shown. We find a mild dependence of the spectral function on the
lengthscales and the magnitude of the RBF kernel, e.g.
$l_{\omega},\,l_{k},\,\sigma$, while the dependence on the asymptotic
parameters, as expected, introduces the majority of the systematic error. The
systematic error, captured by varying these hyperparameters has similar values
to the error propagation on the error of the underlying lattice data. The
parameter space is limited to values of $\nu\geq 1$, since the perturbative
behavior does not reach this far into the IR. However, such solutions are
still represented by the universal RBF kernel, we merely do not want to prefer
the min this region of the spectral function. Additionally, the lengthscale of
the transition is also restricted from below, $l_{\text{\tiny{UV}}}\geq 0.1$,
in order to avoid edges in the reconstructed spectral function. These open
directions in the parameter space can be attributed to inversion being ill-
conditioned. The error on the reconstructed spectral function is finally
computed as a combination of the intrinsic error propagation of the GP, which
takes the lattice data error into account, and the variation of the spectral
function with the parameters given by the NLL scan.

## References

  * Bazavov  _et al._ [2019] A. Bazavov _et al._ (HotQCD), Phys. Lett. B 795, 15 (2019), arXiv:1812.08235 [hep-lat] .
  * Bazavov  _et al._ [2014] A. Bazavov _et al._ (HotQCD), Phys. Rev. D 90, 094503 (2014), arXiv:1407.6387 [hep-lat] .
  * Borsanyi  _et al._ [2014] S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, and K. K. Szabo, Phys. Lett. B 730, 99 (2014), arXiv:1309.5258 [hep-lat] .
  * Pisarski [1982] R. D. Pisarski, Phys. Lett. B 110, 155 (1982).
  * McLerran and Toimela [1985a] L. D. McLerran and T. Toimela, Phys. Rev. D 31, 545 (1985a).
  * Kajantie and Miettinen [1981] K. Kajantie and H. I. Miettinen, Z. Phys. C 9, 341 (1981).
  * Kajantie  _et al._ [1986] K. Kajantie, J. I. Kapusta, L. D. McLerran, and A. Mekjian, Phys. Rev. D 34, 2746 (1986).
  * Gale and Kapusta [1987] C. Gale and J. I. Kapusta, Phys. Rev. C 35, 2107 (1987).
  * Turbide  _et al._ [2004] S. Turbide, R. Rapp, and C. Gale, Phys. Rev. C 69, 014903 (2004), arXiv:hep-ph/0308085 .
  * Adare  _et al._ [2015] A. Adare _et al._ (PHENIX), Phys. Rev. C 91, 064904 (2015), arXiv:1405.3940 [nucl-ex] .
  * Adare  _et al._ [2010] A. Adare _et al._ (PHENIX), Phys. Rev. Lett. 104, 132301 (2010), arXiv:0804.4168 [nucl-ex] .
  * Wilde [2013] M. Wilde (ALICE), Nucl. Phys. A 904-905, 573c (2013), arXiv:1210.5958 [hep-ex] .
  * Adare  _et al._ [2016] A. Adare _et al._ (PHENIX), Phys. Rev. C 94, 064901 (2016), arXiv:1509.07758 [nucl-ex] .
  * Lohner and the ALICE Collaboration [2013] D. Lohner and the ALICE Collaboration, Journal of Physics: Conference Series 446, 012028 (2013).
  * van Hees  _et al._ [2011] H. van Hees, C. Gale, and R. Rapp, Phys. Rev. C 84, 054906 (2011).
  * Shen  _et al._ [2014] C. Shen, U. W. Heinz, J.-F. Paquet, and C. Gale, Phys. Rev. C 89, 044910 (2014), arXiv:1308.2440 [nucl-th] .
  * Churchill  _et al._ [2023a] J. Churchill, L. Du, C. Gale, G. Jackson, and S. Jeon, (2023a), arXiv:2311.06675 [nucl-th] .
  * Churchill  _et al._ [2023b] J. Churchill, L. Du, C. Gale, G. Jackson, and S. Jeon, (2023b), arXiv:2311.06951 [nucl-th] .
  * Chaudhuri and Sinha [2011] A. K. Chaudhuri and B. Sinha, Phys. Rev. C 83, 034905 (2011), arXiv:1101.3823 [nucl-th] .
  * Paquet  _et al._ [2016] J.-F. m. c. Paquet, C. Shen, G. S. Denicol, M. Luzum, B. Schenke, S. Jeon, and C. Gale, Phys. Rev. C 93, 044906 (2016).
  * Vujanovic  _et al._ [2016] G. Vujanovic, J.-F. Paquet, G. S. Denicol, M. Luzum, S. Jeon, and C. Gale, Phys. Rev. C 94, 014904 (2016), arXiv:1602.01455 [nucl-th] .
  * Cuniberti  _et al._ [2001] G. Cuniberti, E. De Micheli, and G. A. Viano, Commun. Math. Phys. 216, 59 (2001), arXiv:cond-mat/0109175 .
  * Meyer [2011] H. B. Meyer, Eur. Phys. J. A 47, 86 (2011), arXiv:1104.3708 [hep-lat] .
  * McLerran and Toimela [1985b] L. D. McLerran and T. Toimela, Phys. Rev. D 31, 545 (1985b).
  * Cè  _et al._ [2020] M. Cè, T. Harris, H. B. Meyer, A. Steinberg, and A. Toniato, Phys. Rev. D 102, 091501 (2020), arXiv:2001.03368 [hep-lat] .
  * Caron-Huot [2009a] S. Caron-Huot, Phys. Rev. D 79, 125009 (2009a), arXiv:0903.3958 [hep-ph] .
  * Brandt  _et al._ [2018] B. B. Brandt, A. Francis, T. Harris, H. B. Meyer, and A. Steinberg, EPJ Web Conf. 175, 07044 (2018), arXiv:1710.07050 [hep-lat] .
  * Sommer [2014] R. Sommer, PoS LATTICE2013, 015 (2014), arXiv:1401.3270 [hep-lat] .
  * Francis  _et al._ [2015] A. Francis, O. Kaczmarek, M. Laine, T. Neuhaus, and H. Ohno, Phys. Rev. D 91, 096002 (2015), arXiv:1503.05652 [hep-lat] .
  * Bazavov  _et al._ [2010] A. Bazavov _et al._ (MILC), PoS LATTICE2010, 074 (2010), arXiv:1012.0868 [hep-lat] .
  * Altenkort  _et al._ [2024] L. Altenkort, D. de la Cruz, O. Kaczmarek, R. Larsen, G. D. Moore, S. Mukherjee, P. Petreczky, H.-T. Shu, and S. Stendebach (HotQCD), Phys. Rev. Lett. 132, 051902 (2024), arXiv:2311.01525 [hep-lat] .
  * Luscher  _et al._ [1997] M. Luscher, S. Sint, R. Sommer, P. Weisz, and U. Wolff, Nucl. Phys. B 491, 323 (1997), arXiv:hep-lat/9609035 .
  * Follana  _et al._ [2007] E. Follana, Q. Mason, C. Davies, K. Hornbostel, G. P. Lepage, J. Shigemitsu, H. Trottier, and K. Wong (HPQCD, UKQCD), Phys. Rev. D 75, 054502 (2007), arXiv:hep-lat/0610092 .
  * Luscher and Weisz [1985] M. Luscher and P. Weisz, Commun. Math. Phys. 97, 59 (1985), [Erratum: Commun.Math.Phys. 98, 433 (1985)].
  * Lüscher and Weisz [1985] M. Lüscher and P. Weisz, Physics Letters B 158, 250 (1985).
  * Heitger and Joswig [2021] J. Heitger and F. Joswig (ALPHA), Eur. Phys. J. C 81, 254 (2021), arXiv:2010.09539 [hep-lat] .
  * Ding  _et al._ [2016] H.-T. Ding, O. Kaczmarek, and F. Meyer, Phys. Rev. D 94, 034504 (2016), arXiv:1604.06712 [hep-lat] .
  * Jackson and Laine [2019] G. Jackson and M. Laine, JHEP 11, 144 (2019), arXiv:1910.09567 [hep-ph] .
  * Laine [2013] M. Laine, JHEP 05, 083 (2013), arXiv:1304.0202 [hep-ph] .
  * Jackson [2019] G. Jackson, Phys. Rev. D 100, 116019 (2019), arXiv:1910.07552 [hep-ph] .
  * Arnold  _et al._ [2001] P. B. Arnold, G. D. Moore, and L. G. Yaffe, JHEP 11, 057 (2001), arXiv:hep-ph/0109064 .
  * Aurenche  _et al._ [2002] P. Aurenche, F. Gelis, G. D. Moore, and H. Zaraket, JHEP 12, 006 (2002), arXiv:hep-ph/0211036 .
  * Baier  _et al._ [1988] R. Baier, B. Pire, and D. Schiff, Phys. Rev. D 38, 2814 (1988).
  * Gabellini  _et al._ [1990] Y. Gabellini, T. Grandou, and D. Poizat, Annals Phys. 202, 436 (1990).
  * Altherr and Aurenche [1989] T. Altherr and P. Aurenche, Z. Phys. C 45, 99 (1989).
  * Braaten  _et al._ [1990] E. Braaten, R. D. Pisarski, and T.-C. Yuan, Phys. Rev. Lett. 64, 2242 (1990).
  * Moore and Robert [2006] G. D. Moore and J.-M. Robert, (2006), arXiv:hep-ph/0607172 .
  * Kapusta  _et al._ [1991] J. I. Kapusta, P. Lichard, and D. Seibert, Phys. Rev. D 44, 2774 (1991), [Erratum: Phys.Rev.D 47, 4171 (1993)].
  * Baier  _et al._ [1992] R. Baier, H. Nakkagawa, A. Niegawa, and K. Redlich, Z. Phys. C 53, 433 (1992).
  * Laine and Rothkopf [2013] M. Laine and A. Rothkopf, JHEP 07, 082 (2013), arXiv:1304.4443 [hep-ph] .
  * Caron-Huot [2009b] S. Caron-Huot, Phys. Rev. D 79, 065039 (2009b), arXiv:0811.1603 [hep-ph] .
  * Aoki  _et al._ [2020] S. Aoki _et al._ (Flavour Lattice Averaging Group), Eur. Phys. J. C 80, 113 (2020), arXiv:1902.08191 [hep-lat] .
  * Vuorinen [2003] A. Vuorinen, Phys. Rev. D 67, 074032 (2003), arXiv:hep-ph/0212283 .
  * Gupta [2004] S. Gupta, Phys. Lett. B 597, 57 (2004), arXiv:hep-lat/0301006 .
  * Ghiglieri  _et al._ [2016] J. Ghiglieri, O. Kaczmarek, M. Laine, and F. Meyer, Phys. Rev. D 94, 016005 (2016), arXiv:1604.07544 [hep-lat] .
  * Bala  _et al._ [2023] D. Bala, S. Ali, A. Francis, G. Jackson, O. Kaczmarek, and T. Ueding, PoS LATTICE2022, 169 (2023), arXiv:2212.11509 [hep-lat] .
  * Kovtun and Starinets [2005] P. K. Kovtun and A. O. Starinets, Phys. Rev. D 72, 086009 (2005), arXiv:hep-th/0506184 .
  * Backus and Gilbert [1968] G. Backus and F. Gilbert, Geophys. J. Int. 16, 169 (1968).
  * Astrakhantsev  _et al._ [2018] N. Y. Astrakhantsev, V. V. Braguta, and A. Y. Kotov, Phys. Rev. D 98, 054515 (2018), arXiv:1804.02382 [hep-lat] .
  * Kanagawa  _et al._ [2018] M. Kanagawa, P. Hennig, D. Sejdinovic, and B. K. Sriperumbudur, (2018), arXiv:1807.02582 [stat.ML] .
  * Rasmussen and Williams [2006] C. E. Rasmussen and C. K. I. Williams,  _Gaussian processes for machine learning_ , Vol. 1 (Springer, 2006).
  * Valentine and Sambridge [2020] A. P. Valentine and M. Sambridge, Geophysical Journal International 220, 1632 (2020).
  * Horak  _et al._ [2022] J. Horak, J. M. Pawlowski, J. Rodríguez-Quintero, J. Turnwald, J. M. Urban, N. Wink, and S. Zafeiropoulos, Phys. Rev. D 105, 036014 (2022), arXiv:2107.13464 [hep-ph] .
  * Pawlowski  _et al._ [2023] J. M. Pawlowski, C. S. Schneider, J. Turnwald, J. M. Urban, and N. Wink, Phys. Rev. D 108, 076018 (2023), arXiv:2212.01113 [hep-ph] .
  * Horak  _et al._ [2023] J. Horak, J. M. Pawlowski, J. Turnwald, J. M. Urban, N. Wink, and S. Zafeiropoulos, Phys. Rev. D 107, 076019 (2023), arXiv:2301.07785 [hep-ph] .
  * Steinwart [2001] I. Steinwart, Journal of Machine Learning Research 2, 67 (2001).
  * Micchelli  _et al._ [2006] C. A. Micchelli, Y. Xu, and H. Zhang, Journal of Machine Learning Research 7, 2651 (2006).
  * Mazur  _et al._ [2023] L. Mazur _et al._ (HotQCD), (2023), arXiv:2306.01098 [hep-lat] .