# Lattice QCD estimates of thermal photon production from the QGP

Sajid Ali, Dibyendu Bala, Anthony Francis, Greg Jackson, Olaf Kaczmarek, Jonas Turnwald, Tristan Ueding, Nicolas Wink

► To cite this version:

Sajid Ali, Dibyendu Bala, Anthony Francis, Greg Jackson, Olaf Kaczmarek, et al.. Lattice QCD estimates of thermal photon production from the QGP. Physical Review D, 2024, 110 (5), pp.054518. 10.1103/PhysRevD.110.054518. hal-04529120

HAL Id: hal-04529120

https://hal.science/hal-04529120v1

Submitted on 28 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Lattice QCD estimates of thermal photon production from the QGP

Sajid Ali,¹,⁶ Dibyendu Bala,¹ Anthony Francis,² Greg Jackson,³,⁴ Olaf Kaczmarek,¹ Jonas Turnwald,⁵ Tristan Ueding,¹ and Nicolas Wink⁵

(HotQCD Collaboration)

¹Fakultät für Physik, Universität Bielefeld, D-33615 Bielefeld, Germany

²Institute of Physics, National Yang-Ming Chiao Tung University, 30010 Hsinchu, Taiwan

³Institute for Nuclear Theory, University of Washington,
Box 351550, Seattle, Washington 98195-1550, USA

⁴SUBATECH (Nantes Université, IMT Atlantique, IN2P3/CNRS),
4 rue Alfred Kastler, La Chantrerie BP 20722, 44307 Nantes, France

⁵Institut für Kernphysik, Technische Universität Darmstadt, D-64289 Darmstadt, Germany

⁶Government College University Lahore, Department of Physics, Lahore 54000, Pakistan

<figure><img src="image_2.png" /></figure>

(Received 17 April 2024; accepted 9 September 2024; published 27 September 2024)

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 nonperturbative effects, the lattice correlator is compared with the corresponding NLO + LPM$^{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.

DOI: 10.1103/PhysRevD.110.054518

## I. INTRODUCTION

One of the important predictions of quantum chromo-dynamics (QCD), the theory of strong interactions, is that the QCD coupling constant decreases with the energy scale. This characteristic makes QCD strongly coupled (nonperturbative) at low energies, and weakly coupled at very high energies. As a result, nonperturbative 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$ MeV [1–3].

Experimentally, this phase has been studied in ultra-relativistic 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–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

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Funded by SCOAP³.
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–15]. However, understanding the origin of this
phenomenon from first principles poses a significant
theoretical challenge [16].

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 [17–19]. The anisotropic flow of direct photons
and dileptons has been shown to be sensitive to the initial
conditions and shear viscosity [20–22]. Since electro-
magnetic 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 under-
standing of the photon production mechanism from each
stage of the plasma evolution. Since hydrodynamic expan-
sion assumes local thermal equilibrium, computation of the
photon production rate at each stage requires the use of
equilibrium thermal field theory calculations.

Numerous studies have investigated the thermal photon
rate, both perturbatively and nonperturbatively, using
lattice techniques. The perturbative evaluation of the rate
has a long history, and was regarded as settled in the early
1990s when two groups independently arrived at the same
answer [23,24]. Those works implemented resummed
thermal fermion propagators to account for screening,
the first step needed to obtain a finite result. It was then
discovered that those treatments were incomplete and
additional processes (naively of higher order) were incor-
porated to establish the full $O(\alpha_{em}\alpha_s)$ result [25,26]. More
recently, it was understood how to take advantage of a
Euclideanization property for thermal light-cone correlators
to attain relative $O(\sqrt{\alpha_s})$ corrections [27].

On the lattice QCD side, a major challenge is the fact that
the thermal production of photons requires real-time
information. This can only be achieved after analytic
continuation of Euclidean data, which is formally an
ill-conditioned inverse operation [28,29]. Consequently,

various techniques and models have been employed for
spectral reconstruction. In Ref. [30], the photon rate was
estimated using the vector channel correlator, wherein the
UV part of the spectral function was calculated perturba-
tively, while the IR part was fixed using lattice data. A
similar method was employed in another study [31], where
the photon rate was estimated from the transverse channel
correlator. A new correlator was utilized in Ref. [32], based
on the difference between the transverse and longitudinal
channels. This approach is advantageous because it sup-
presses the UV part of the spectral function. Recently, a
novel idea has been proposed for calculating the photon
rate. In this method, the photon production rate is studied
using the imaginary momentum correlator, which does not
require spectral reconstruction [33,34].

In this paper, following [32], we estimate the thermal
photon rate from the QGP using the difference between the
transverse and longitudinal correlators. For the spectral
reconstruction, we used various techniques, 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 MeV (1.5*T<sub>c</sub>), and (2 + 1)-flavor QCD at a temperature
of 220 MeV (1.22*T<sub>pc</sub>). 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.

For the full QCD case, we use highly improved staggered quark (HISQ) configurations for the background gauge fields. This is advantageous because generating HISQ configurations is easier than generating dynamical Wilson fermion configurations, especially at the fine lattice spacing we are considering. One reason for this is that staggered fermions preserve a remnant of chiral symmetry, which protects the Dirac operator from having small real eigenvalues responsible for exceptional configurations. Additionally, chiral symmetry prevents the quark mass from having any additive renormalization constant, making quark mass tuning much easier. In the valence sector, we use clover-improved Wilson fermions. The reason for this choice is that HISQ correlators contain both oscillating and nonoscillating modes, which complicates the spectral reconstruction method at finite temperatures. The downside of this approach is the presence of mixed-action cutoff effects. However, the mixed-action effect, quantified in terms of the relative mass difference [35], turns out to be small at the lattice spacing we are working with. We therefore expect these mixed-action effects to vanish in the continuum limit.

The paper is organized into the following sections: In
Sec. II, we will briefly discuss the thermal photon rate
and the correlator that we used to calculate it. In Sec. III,
we present the lattice details and perform the continuum extrapolation of quenched lattice data. In Sec. 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 Sec. V. Finally, in Sec. 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 photon production rate is defined as the number of photons radiated from the plasma per unit time and per unit volume. We consider the thermal photon production rate due to $N_f = 3$ degenerated quark flavors, whose total electromagnetic current that couples to the photon is given by $e \sum_f Q_f \bar{\psi}_f \gamma^\mu \psi_f$, where $Q_f$ is the quark charge fraction and $\psi_f$ is the quark field of flavor f. If we consider the case $N_f = 3$ where $f = u, d, s$, then $\sum_f Q_f = 0$ and $\sum_f Q_f^2 = \frac{2}{3}$. Under these conditions, the disconnected<fnref n="1" /> quark contribution vanishes exactly and therefore to leading order in the QED coupling $\alpha_{em}$, this rate can be calculated as

$$ \frac{d\Gamma_{\gamma}}{d^3\vec{k}} = -\frac{\alpha_{em}n_b(k)}{2\pi^2 k} \left\{ \sum_{i=1}^{N_f} Q_i^2 \right\} g^{\mu\nu} \rho_{\mu\nu}(\omega = |\vec{k}|, \vec{k}), \quad (1) $$

where $n_b(\omega) = 1/(e^{\omega/T} - 1)$ is the Bose distribution function [37]. We take the Minkowski metric as $g_{\mu\nu} = \text{diag}(1, -1, -1, -1)$. In this expression, the $\rho_{\mu\nu}$ is the connected part of the vector current spectral function for a single flavor.

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,<fnref n="2" /> i.e., $\rho_{\mu\nu}(\omega, \vec{k}) = \text{Im}[G_{\mu\nu}^E(\omega_n \rightarrow -i(\omega + i0^+), \vec{k})]$, where the Euclidean current correlator is given by

$$ G_{\mu\nu}^{E}(\omega_n, \vec{k}) = \int d^3\vec{x} \int_0^\beta d\tau e^{i(\omega_n \tau - \vec{k} \cdot \vec{x})} \langle J_\mu(\tau, \vec{x}) J_\nu(0, \vec{0}) \rangle, \quad (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_{\mu\nu}^E(\tau, \vec{k}) = T \sum_n \exp(i\omega_n\tau) G_{\mu\nu}^E(\omega_n, \vec{k})$. This correlator can be written in terms of the spectral function as

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

As mentioned before, extracting the spectral function from the Euclidean correlator is an ill-conditioned problem. Therefore, when inverting Eq. (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 Eq. (1), is proportional to $\omega^2$ in the ultraviolet (UV) domain. As a result, the lattice correlator, determined by the integration in Eq. (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. [32], 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_{\mu\nu}^{T} \rho_{T}(\omega, \vec{k}) + P_{\mu\nu}^{L} \rho_{L}(\omega, \vec{k}). \quad (4) $$

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

$$ \begin{aligned} P_{ij}^{T}(\omega, \vec{k}) &= -g_{ij} - \frac{k_i k_j}{\vec{k}^2}, \\ P_{0i}^{T} &= P_{i0}^{T} = P_{00}^{T} = 0, \\ P_{\mu\nu}^{L}(\omega, \vec{k}) &= -g_{\mu\nu} + \frac{K_{\mu} K_{\nu}}{K^2} - P_{\mu\nu}^{T}(\omega, \vec{k}), \end{aligned} \quad (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\{\rho_T(\omega, \vec{k}) - \rho_L(\omega, \vec{k})\}. \quad (6) $$

<footnote n="1">¹In terms of the $\alpha_s$-expansion, the disconnected contribution only starts at $O(\alpha_s^3)$ [36].</footnote>

<footnote n="2">²These discrete energies are denoted $\omega_n = 2\pi T n$, where $n \in \mathbb{Z}$.</footnote>
At zero temperature, because of the Ward identity and Lorentz invariance, the T-L spectral function vanishes, $\rho_H(\omega, \vec{k}) = 0$. Therefore, nonzero $\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$ [38,39]. One consequence of the OPE, is that this spectral function satisfies a sum rule,

$$ \int_0^\infty d\omega \, \omega \rho_H(\omega) = 0. \qquad (7) $$

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

$$ \frac{d\Gamma_{\gamma}}{d^3\vec{k}} = \frac{\alpha_{\text{em}}n_b(k)\chi_q}{\pi^2} \left\{ \sum_{i=1}^{N_f} Q_i^2 \right\} D_{\text{eff}}(k), \qquad (8) $$

where the effective diffusion coefficient is defined by

$$ D_{\text{eff}}(k) = \frac{\rho_H(|\vec{k}|, k)}{2\chi_q|\vec{k}|}. \qquad (9) $$

Here $\chi_q$ is the quark number susceptibility [defined later in Eq. (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} = (\frac{2\pi n_x}{L}, \frac{2\pi n_y}{L}, \frac{2\pi n_z}{L})$, 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 $x$ direction i.e., $\vec{k} = (k, 0, 0)$, then the correlator for $\rho_H$ is given by [39]

$$ G_H(\tau, k_x = k) = G_{yy}(\tau, k_x = k) + G_{zz}(\tau, k_x = k) - 2(G_{xx}(\tau, k_x = k) - G_{\tau\tau}(\tau, k_x = k)). \qquad (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} \{ G_H(\tau, k_x = k) + G_H(\tau, k_y = k) + G_H(\tau, k_z = k) + (\tau \to \beta - \tau) \}. \qquad (11) $$

We calculated the $G_H$ correlator using Clover-improved Wilson fermions for quenched configurations at a

temperature of ~1.5$T_c$ and (2 + 1)-flavor QCD configurations at a temperature of ~1.2$T_{pc}$. The scale for quenched QCD has been set by $r_0 = 0.47$ fm [40] and the $T_c$ is obtained by $T_c = 0.7457/r_0$ [41]. In the case of full QCD, the scale is determined by $r_1 = 0.3106$ fm and $T_{pc} = 180$ MeV was obtained from disconnected chiral susceptibility, cf. Ref. [42], as in Ref. [43].

We generated the quenched configurations using a standard Wilson gauge action with heat-bath and over-relaxation 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. [44]. The clover improvement constant, also known as the Sheikholeslami-Wohlert coefficient ($c_{SW}$), was estimated from the nonperturbative parametrization of Ref. [44].

For full QCD, we used gauge field configurations generated by the HotQCD collaboration using highly improved staggered quark (HISQ) [45] and a tree-level improved Lüscher-Weisz gauge action [46,47]. We used the RHMC algorithm [48] for the production of configurations and every fifth configuration has been saved for measurement. These configurations have also been used in [43], and the auto-correlation time is estimated to be 10–30 configurations from a gradient flowed Polyakov loop. We have also looked at the gradient flowed topological charge, which is known to exhibit a larger auto-correlation time. We find that our gauge configurations are mostly in the trivial topological sector. This is expected, as at high temperatures, the topological susceptibility is small. However, the correlation functions we consider here are not too sensitive to topology, and therefore, we expect the correlator to exhibit only a mild dependence on the topological charge. This behavior has been observed in Ref. [32]. We tuned the $\kappa$ value to match the pion mass, which in this study was unphysical and equal to 320 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 I. Our implementation of the current operator on the lattice is not $O(a)$ improved, but as discussed in Ref. [39], the effect of this improvement is quark mass suppressed in the chirally symmetric phase. For the quenched case, the effect of the improvement coefficient is small in the coupling range we are working on [49]. Therefore, we expect the correction to the continuum result should be $O(a^2)$ in both quenched and full QCD lattice correlator.
<table><caption>TABLE I. 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$ MeV for $N_f = 0$ case and the pseudo critical temperature is $T_{pc} = 180$ MeV for $N_f = 2 + 1$ case.</caption><thead><tr><th>N<sub>f</sub></th><th>T</th><th>&beta;<sub>0</sub></th><th>&kappa;</th><th>N<sub>&sigma;</sub><sup>3</sup> &times; N<sub>&tau;</sub></th><th>Confs</th><th>a<sup>-1</sup> [GeV]</th></tr></thead><tbody><tr><td rowspan="3">0</td><td rowspan="3">1.5T<sub>c</sub></td><td>7.394</td><td>0.13407</td><td>120<sup>3</sup> &times; 30</td><td>1950</td><td>14.1</td></tr><tr><td>7.192</td><td>0.13440</td><td>96<sup>3</sup> &times; 24</td><td>2000</td><td>11.3</td></tr><tr><td>7.035</td><td>0.13467</td><td>80<sup>3</sup> &times; 20</td><td>1824</td><td>9.39</td></tr><tr><td>2 + 1</td><td>1.2T<sub>pc</sub></td><td>8.249</td><td>0.13515</td><td>96<sup>3</sup> &times; 32</td><td>1750</td><td>7.04</td></tr></tbody></table>

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 in the continuum is defined as

$$ \chi_q = \int_0^{1/T} d\tau G_{\tau\tau}(\tau, \vec{0}). \qquad (12) $$

Since $\chi_q$ and $G_H$ have the same renormalization coefficient, the ratio $G_H/(2\chi_q T)$ does not require renormalization. On the lattice, however, due to charge conservation, $G_{00}(\tau, \vec{0})$ is essentially constant in $\tau$, except for the short-distance part of $G_{00}(\tau, \vec{0})$ related to cutoff effects. As a result, we calculate $G_H/(2G_{00}(\tau = 1/2T, \vec{0}))$ for an estimate of $G_H/(2\chi_q T)$.

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 coarser 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}. \qquad (13) $$

We illustrated the fitting of the above ansatz for a few values of $\tau T$ in left panel of Fig. 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 cutoff effect on this correlator at the momenta $k/T = 1.57$ and 3.14 is barely visible; however, with increasing momentum, cutoff effects start showing up when $\tau T \lesssim 0.25$. In contrast, the vector channel shows a significant cutoff effect due to the presence of a large UV component [50]. This also indicates that T-L correlators are dominated by the IR part of the spectral function.

## IV. COMPARISON WITH PERTURBATIVE ESTIMATE

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 [25,51–54]. However, at the phenomenologically interesting temperature range, QCD is strongly coupled and weak-coupling calculations may not be sufficient. Therefore, nonperturbative 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 nonperturbative effects.

<figure><img src="image_3.png" /><figcaption>FIG. 1. Euclidean correlator from lattice simulations in the T-L channel, cf. Eq. (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.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 2. Perturbative spectral function in the T-L channel computed using the procedure from Ref. [51] at NLO + LPM$^{LO}$, as a function of $\omega$ for several values of spatial momenta. (LPM$^{LO}$ means that near the light cone LPM resummation is performed at leading order.) Here the renormalization scale was fixed at $\mu = \mu_{\text{opt}}$ as described in the main text. The left panel shows the spectral function with no dynamical quarks at $T = 470$ MeV, and the right panel shows the spectral function in $N_f = 2 + 1$ flavor QCD at $T = 220$ MeV.</figcaption></figure>

At leading-order (LO), the spectral function is determined from the process $q\bar{q} \to \gamma^*$. At 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 [55–57]. 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| \gtrsim (\pi T)^2$), perturbative calculations are well-behaved [52,53]. 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 [58,59]. Therefore, a separate treatment is required to calculate the spectral function close to the light cone.

When $M^2 \lesssim 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 [23,24]. Additionally, multiple scatterings due to soft gluons can occur with the quark that emits the photon.<fnref n="3" /> To include this contribution, a technique called Landau-Pomeranchuk-Migdal (LPM) resummation is used in this regime [25,54]. 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. [51]. 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 [60]. 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_{\text{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 [61]. The running coupling then has been calculated by using $\Lambda^{\text{MS}} = 0.8 T_c$ for quenched [41] and $\Lambda^{\text{MS}} = 339$ MeV $\simeq 1.88 T_{\text{pc}}$ for full QCD [62] (these choices have uncertainties $\sim 10...25\%$).

The resulting spectral function is plotted for various momenta in Fig. 2 for both quenched and full QCD. As expected, the spectral function is UV suppressed, and it has been shown [51] 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 [38].

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 \pm 0.012)T^2$ for the case of $N_f = 0$ and $\chi_q = 0.842_{-0.03}^{+0.05}T^2$ for $N_f = 2 + 1$. The value for $N_f = 0$ was estimated using the Lüscher nonperturbative method [50], while for the case of $N_f = 2 + 1$, we used the perturbative EQCD estimate of order $g^6 \log(g)$ [63]. The errorbar on the nonperturbative estimate is statistical in nature [50], whereas for the perturbative estimate, the error bar is obtained by varying the renormalization scale between $\mu = 0.5\mu_{\text{min}}$ to $\mu = 2.0\mu_{\text{min}}$, where $\mu_{\text{min}} = 9.082T$ is the

<footnote n="3">³These 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.</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 3. Euclidean correlator in T-L channel, calculated from the NLO + LPMLO 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 Eq. (3), while data points represent lattice correlators for quenched on the left and full QCD on the right.</figcaption></figure>

minimum sensitivity scale. For $N_f = 3$, $\chi_q$ for the physical pion mass was calculated in Ref. [64], which provides an estimate $\chi_q = (0.772 \pm 0.015)T^2$ at the temperature $1.2T_c$. In Ref. [65], a hard thermal loop (HTL) resummed perturbative estimate gives a susceptibility $\chi_q = 0.888_{-0.063}^{+0.023}T^2$. However, using these values for the $\chi_q$, the qualitative conclusion remains unchanged. We should also like to mention we need $\chi_q$ only to compare the perturbative correlator. The $D_{eff}$ can be directly calculated from the correlator $G_H/(2\chi_q T)$. The comparison between the correlator obtained from the perturbative spectral function and the lattice calculation is illustrated in Fig. 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 nonperturbative difference remains. On the other hand, in full QCD, nonperturbative 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.5T_c$ and $\sim 1.2T_{pc}$ respectively. The $N_f = 0$ case temperature in this “scaled temperature” unit is 25% higher than the $N_f = 2 + 1$ case. Therefore, Fig. 3 also illustrates that by going 25% higher in temperature, one sees a significant suppression in the nonperturbative effects.

V. SPECTRAL RECONSTRUCTION

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 [30,32,39,66]. As discussed in Sec. II, the known physical constraints on the spectral function contain the sum rule Eq. (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 [67].

A. 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 \gtrsim k, \pi T$). The ultraviolet region $\omega \ge \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 \le \omega_0$ is modeled by $\rho_<$, for which we adopt the same polynomial put forward in Ref. [30]. The two regions are matched continuously,

$$ \rho_H^{\text{poly}}(\omega) = \rho_<(\omega)\Theta(\omega_0 - \omega) + \rho_>(\omega)\Theta(\omega - \omega_0), \quad (14) $$

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

$$ \rho_<(\omega_0) = \rho_>(\omega_0) = \beta, \quad (15) $$

$$ \rho'_<(\omega_0) = \rho'_>(\omega_0) = \gamma. \quad (16) $$

These considerations lead to the following model of the spectral function
$$ \rho_{<} = \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, \quad (17) $$

$$ \rho_{>} = -\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, \quad (18) $$

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} \le 0$. The resulting fits of the lattice data with this model of the spectral function is shown in Fig. 4 for quenched and full QCD, while the spectral function obtained from this fit is illustrated in Fig. 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 Fig. 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_{\text{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 II for quenched and in Table III for full QCD. The results shown are from continuum extrapolated dataset 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. [32] and is given by

$$ \rho_H^{\text{Padé}}(\omega) = \frac{A \tanh\left(\frac{\omega}{2T}\right) (1 + B\omega^2)}{(\omega^2 + a^2)[(\omega + \omega_0)^2 + b^2][(\omega - \omega_0)^2 + b^2]} \quad (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 [29]. The remaining part of the spectral function is inspired by qualities from AdS/CFT (see, e.g., Ref. [68]) 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 Eq. (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. [39]. Therefore, we quote the variation of the minimum and maximum $D_{\text{eff}}T$ of

<figure><img src="image_6.png" /><figcaption>FIG. 4. The Euclidean correlator in the T-L channel, illustrating the result of fitting the polynomial ansatz in Eq. (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.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIG. 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_{opt}$ (red points).</figcaption></figure>

the distribution within the error bars (rather than the standard deviation from the mean values). The estimated $D_{eff}$ can be found in Table II for quenched and Table III for full QCD. Sample spectral functions for the two extreme cases are plotted in Fig. 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 compared to the perturbative spectral function.

## B. Backus-Gilbert estimate

The Backus-Gilbert (BG) [69] 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^{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^{BG}(\omega, \vec{k}) = \int d\omega^* \delta(\omega, \omega^*) \rho(\omega^*, \vec{k}). \quad (20) $$
<table><caption>TABLE II. Extracted values of $D_{\text{eff}}T$ from different methods for $N_f = 0$.</caption><thead><tr><th>k/T</th><th>Polynomial</th><th>Padé</th><th>Backus-Gilbert</th><th>GP regression</th></tr></thead><tbody><tr><td>1.57</td><td>0.549 ± 0.035</td><td>0.469 ± 0.368</td><td>0.302 ± 0.175</td><td>0.552 ± 0.129</td></tr><tr><td>3.14</td><td>0.161 ± 0.092</td><td>0.129 ± 0.123</td><td>0.119 ± 0.057</td><td>0.140 ± 0.039</td></tr><tr><td>4.71</td><td>0.087 ± 0.048</td><td>0.089 ± 0.067</td><td>0.054 ± 0.031</td><td>0.061 ± 0.016</td></tr><tr><td>6.28</td><td>0.054 ± 0.018</td><td>0.058 ± 0.044</td><td>0.022 ± 0.018</td><td>0.048 ± 0.012</td></tr><tr><td>7.85</td><td>0.041 ± 0.007</td><td>0.066 ± 0.066</td><td>0.016 ± 0.016</td><td>0.033 ± 0.010</td></tr><tr><td>9.42</td><td>0.031 ± 0.007</td><td>0.074 ± 0.045</td><td>0.013 ± 0.013</td><td>0.020 ± 0.009</td></tr></tbody></table>

<table><caption>TABLE III. Extracted values of $D_{\text{eff}}T$ from different methods for $N_f = 2 + 1$.</caption><thead><tr><th>k/T</th><th>Polynomial</th><th>Padé</th><th>Backus-Gilbert</th><th>GP regression</th></tr></thead><tbody><tr><td>2.09</td><td>0.369 ± 0.040</td><td>0.206 ± 0.085</td><td>0.197 ± 0.127</td><td>0.404 ± 0.219</td></tr><tr><td>4.19</td><td>0.102 ± 0.034</td><td>0.062 ± 0.058</td><td>0.063 ± 0.047</td><td>0.055 ± 0.048</td></tr><tr><td>6.28</td><td>0.031 ± 0.013</td><td>0.010 ± 0.010</td><td>-0.011 ± 0.047</td><td>0.004 ± 0.027</td></tr></tbody></table>

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 Eq. (3) is rewritten as

$$G_H(\tau, \vec{k}) = \int_0^\infty d\omega \frac{\rho(\omega, \vec{k})}{f(\omega, \vec{k})} \tilde{K}(\omega, \tau), \quad (21)$$

where $\tilde{K}(\omega, \tau) = \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

$$\begin{aligned} \frac{\rho_H^{\text{BG}}(\omega, \vec{k})}{f(\omega, \vec{k})} &= \sum_i q_i(\omega, \vec{k}) G_H(\tau_i, \vec{k}) \\ &= \int_0^\infty d\omega^* \delta(\omega^*, \omega) \frac{\rho_H(\omega^*, \vec{k})}{f(\omega^*, \vec{k})}, \end{aligned} \quad (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)$ (with respect to $q_i$), a measure of the width of the resolution function

$$A(\omega) \equiv \int d\omega^* (\omega - \omega^*)^2 \delta(\omega, \omega^*)^2 = q^\dagger(\omega) W(\omega) q(\omega), \quad (23)$$

where $[W]_{ij}(\omega) \equiv \int_0^\infty d\omega^* (\omega - \omega^*)^2 \tilde{K}(\omega^*, \tau_i) \tilde{K}(\omega^*, \tau_j)$.

<figure><img src="image_8.png" /><figcaption>FIG. 6. Spectral function $\rho_H$ obtained from the Padé-like ansatz in Eq. (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.</figcaption></figure>
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 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 Eq. (26).</figcaption></figure>

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 $[q]_i$ is given by

$$q(\omega) = \frac{W^{-1}(\omega)r}{r^T W^{-1}(\omega)r}, \qquad (24)$$

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

One of the popular regularization schemes is to take

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

where $S$ is the covariance matrix of the lattice correlator. Other regularization schemes also exist in literature, e.g., the Tikhonov regularization [70] where $S$ is replaced by the identity matrix. In this paper, we have used a diagonal covariance matrix for $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, \qquad (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 Fig. 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 Eq. (9), which is plotted on the left panel of Fig. 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 Fig. 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 \lesssim 2$.

The reason for this can be understood by looking at the resolution function, which is depicted in Fig. 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
<figure><img src="image_11.png" /><figcaption>FIG. 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.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 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$.</figcaption></figure>

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 II for quenched and Table III 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.

### C. Gaussian process regression

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 applica-
tions, primarily in modeling and classification problems;
see Ref. [71] for a review and Ref. [72] for a textbook
introduction. Because no explicit functional basis is
required, GPR offers a very general approach to interpolat-
ing 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 nonparametric
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. [73,74] 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 distribu-
tion 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 trans-
formation, cf. Eq. (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. [73] and applied to
a number of numerical reconstructions in the context of
QCD, see, e.g., Refs. [74–76].

GPR has two main appeals, when considering inverse
problems. First, 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.,
Eq. (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 (C4). 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 [77,78]. However, when having additional infor-
mation 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.

1. GPR reconstruction details

Since we have data for several spatial momenta, the
reconstruction of the spectral function is not only per-
formed in the ω, 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 direc-
tion. For the reconstruction, we include the lattice data
for all the available spatial momenta as well as the sum
rule Eq. (7).

As a very common GP-kernel choice, the RBF kernel
equation (C4) 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, ρ_H^UV(ω) ∝ 1/ω^4, the GP-kernel
is modified in the UV regime to restrict the functional basis
to the asymptotic behavior resulting in Eq. (C5). 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_ω, l_k, σ, l_{UV}, v. The lenghtscale in k and
ω direction is controlled by l_k and l_ω respectively, while σ
gives a prior estimate on the variance of the GP. The
position of the transition from the RBF kernel toward the
UV kernel is controlled by ω₀ = √k² + v(πT)², where v is
varied, in analogy to the polynomial fit ansatz, and l_{UV} is
the smoothness of the transition between the two kernels.

These parameters are optimized by minimizing the
associated negative log-likelihood (NLL) Eq. (C8). Since
this optimization is performed in a five dimensional
parameter space and the parameters are not fully indepen-
dent, this optimization generally does not converge con-
sistently. Additionally, open directions in the parameter
space, e.g., toward vanishing length scales are possible.
<figure><img src="image_13.png" /><figcaption>FIG. 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.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIG. 11. Effective diffusion coefficient $D_{\text{eff}}$ obtained, as in Eq. (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.</figcaption></figure>

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 Appendix 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 Fig. 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 Fig. 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 toward smaller or higher momenta comes with large uncertainties. Ultimately, when extrapolating toward $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 II for quenched and in Table III for full QCD.

## VI. COMPARISON OF $D_{\text{eff}}T$

In the previous sections we discussed various models and methods for the spectral reconstruction from lattice data.
<figure><img src="image_15.png" /><figcaption>FIG. 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.</figcaption></figure>

Now, we compare the photon production rate obtained from these different methods.

In Fig. 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 from$^4\mu = 0.5\mu_{\text{opt}}$ to $\mu = 2\mu_{\text{opt}}$ (where the optimal scale is defined in Sec. IV). The error band on the perturbative estimate is slightly smaller for quenched QCD at $T = 470$ MeV compared to full QCD at $T = 220$ MeV, which may indicate that the perturbative framework is better behaved in this case.

In the case of quenched QCD, our nonperturbative 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 nonperturbative 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 nonperturbative 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$.

Qualitatively, we can compare our dynamical results with other lattice studies of $D_{\text{eff}}$. Our results are in qualitative agreement with the study in Ref. [32] at a temperature of 254 MeV (1.2 $T_{\text{pc}}$). In both studies, at the lowest momentum ($k = 0.4$ GeV), $D_{\text{eff}}$ ranges from 0.5 to $2.7 \text{ GeV}^{-1}$. At the highest available momentum, the results are consistent with zero within the error bars. However, the error on our $D_{\text{eff}}$ is smaller for the highest available momentum compared to Ref. [32]. This could be due to the different systematics involved in the spectral reconstruction methods, as well as the fact that the results in Ref. [32] are from continuum-extrapolated 2-flavor QCD correlators, whereas the dynamical results in this paper are from (2 + 1)-flavor QCD at a fixed lattice spacing.

In Ref [31], the photon production rate is estimated from the transverse channel of the correlator. In this case, the results are consistent at smaller momentum. However, $D_{\text{eff}}$ estimated here shows a smaller value at the highest

<footnote>$^4$The smaller value of $\mu$ corresponds to an upper estimate of $D_{\text{eff}}$, due to the larger QCD coupling being probed. In our case, the extreme choice of $\mu = 0.5\mu_{\text{opt}}$ leads to $\alpha_s \simeq 0.5$ which is where most of the sensitivity to scale variation originates.</footnote>
momentum, which could again be due to different system-atics involved, similar to the reason mentioned above.

## VII. CONCLUSION

We estimated the thermal photon production rate from the Euclidean lattice correlators in quenched QCD ($T = 470$ MeV) and full QCD ($T = 220$ MeV) with an unphysical pion mass of 320 MeV. Motivated by Ref. [32], 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 signifi-cant 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 nonperturbative effects in this corre-lator by comparing it with its perturbative estimate, shown in Fig. 3. 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 [51].

We observe that in quenched QCD, the nonperturbative 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 pertur-bative correlator. As we move toward larger momenta, the lattice data also begins to exhibit curvature similar to the perturbative correlator. However, it appears that nonper-turbative 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 recon-struction. 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 Eq. (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 Fig. 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.<fnref n="5" /> Here we also observed that except the smallest momen-tum, for the quenched case the resulting nonpertur-bative 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 [32]. A sample spectral function is shown in Fig. 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 Eq. (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 with respect to 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 regres-sion 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 continu-ous prediction for the effective diffusion coefficient, see Fig. 11. By considering the UV asymptotic behavior from the OPE and the sum rule, we can reduce the vari ance on the GP posterior in order to make meaningful predictions.

The photon production rate in terms of $D_{\text{eff}}$ is plotted in Fig. 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 esti-mate. In the full QCD case, the photon production rate becomes smaller than the perturbative estimate and becomes consistent with zero at a momentum of approx-imately 0.8 GeV...1.6 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

<footnote n="5">⁵The uncertainty with respect to the variation of $\omega_0$ is taken to be the systematic uncertainty on our result.</footnote>
<figure><img src="image_16.png" /><figcaption>FIG. 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 NLO + LPM$^{LO}$ perturbative correlator. In the right panel, the reconstructed spectral function is displayed along with the exact NLO + LPM$^{LO}$ spectral function.</figcaption></figure>

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.

All data from our calculations, presented in the figures of this paper, can be found in [79].

## 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 [80]. 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 No. 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 No. 111-2112-M-A49-018-MY2.

## APPENDIX A: MOCK-ANALYSIS OF PERTURBATIVE DATA WITH POLYNOMIAL ANSATZ

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 NLO + LPM$^{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 \approx 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 Fig. 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 Fig. 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 approximately predict the 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 FINITED 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 Fig. 14, while the same comparison for the GPR reconstruction can be found in
<figure><img src="image_17.png" /><figcaption>FIG. 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$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIG. 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$.</figcaption></figure>

Fig. 15. We observe that for both methods the results for small spatial momenta agree remarkably well, as different values remain 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.

## APPENDIX C: INVERSION WITH GPR

Here, we will give a short introduction to the reconstruction of spectral functions with GPR. For additional details and a more general discussion, we refer to Refs. [72–74] for an introduction to GPs in general. Note that in the following part, we closely follow the introduction to GP reconstruction found in Ref. [76] and modify the notation to match the current reconstruction in order to ensure a self-contained and reproducible presentation.

To start, we define a GP prior for the spectral function that is reconstructed. This assumes some Gaussian distribution over different realizations of the spectral function.

The GP prior is generally written as

$$ \rho(\omega) \sim \mathcal{GP}(\mu(\omega), C(\omega, \omega)), \quad (C1) $$

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.

For the reconstruction, the GP prior is then conditioned on the available data, e.g., observations of the propagator $G_{H,i}$ at discrete Euclidean times $\tau_i \equiv [\tau]_i$. This means that we compute the conditional posterior distribution of the GP, which has a closed analytic form and is given by [74]

$$ \begin{aligned} \rho(\omega) | G(\tau) &\sim \mathcal{GP}(\mathbf{w}^{\top}(\omega)(\mathbf{W} + \sigma_n^2 \mathbf{1})^{-1} G_H(\tau), \\ C(\omega, \omega) - \mathbf{w}^{\top}(\omega)(\mathbf{W} + \sigma_n^2 \mathbf{1})^{-1} \mathbf{w}(\omega)), \end{aligned} \quad (C2) $$

where

$$ \begin{aligned} [\mathbf{w}]_i(\omega) &= \int d\omega' K(\tau_i, \omega') C(\omega', \omega), \\ [\mathbf{W}]_{ij} &= \int d\omega' d\omega'' K(\tau_i, \omega') K(\tau_j, \omega'') C(\omega', \omega''). \end{aligned} \quad (C3) $$
and $K(\tau, \omega) \equiv \frac{\cosh(\omega(1/2T-\tau))}{\pi \sinh(\omega/2T)}$. Equation (C3) 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 of the underlying function 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

$$C_{\text{RBF}}(\omega, \omega') = \sigma^2 \exp\left(-\frac{(\omega - \omega')^2}{2l^2}\right), \quad (C4)$$

where $\sigma$ controls the overall magnitude and $l$ is a generic lengthscale within the GP. The RBF kernel falls under the class of so-called universal kernels [77,78]. 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. [76]. 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{UV}}(\omega) \propto 1/\omega^4$ of the spectral function, this can be achieved by taking the asymptotic kernel

$$C_{\text{UV}}(\omega, \omega') = (\omega\omega')^{-4}. \quad (C5)$$

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

$$\theta^{\pm}(\omega; l_{\text{UV}}, \omega_0) = \frac{1}{1 + \exp(\pm(\omega - \omega_0))/l_{\text{UV}}}, \quad (C6)$$

where $\omega_0$ is the position and $l_{\text{UV}}$ controls the steepness of the transition from the universal to the UV kernel. This results in the full kernel

$$C(\omega, \omega') = \theta^{-}(\omega)\theta^{-}(\omega')C_{\text{RBF}}(\omega, \omega') + \theta^{+}(\omega)\theta^{+}(\omega')C_{\text{UV}}(\omega, \omega'). \quad (C7)$$

The parameters of the soft step function, i.e., the position and the lengthscale of the transition, and of the RBF kernel

are subject to likelihood optimization, described in the following.

# 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 Eq. (C3), is given by

$$-\log p(G_H(\boldsymbol{\tau})|\boldsymbol{\sigma}) = \frac{1}{2} G_H(\boldsymbol{\tau})^\mathrm{T} (\mathbf{W}_\sigma + \sigma_n^2 \mathbf{1})^{-1} G_H(\boldsymbol{\tau}) \\ + \frac{1}{2} \log \det(\mathbf{W}_\sigma + \sigma_n^2 \mathbf{1}) + \frac{N}{2} \log 2\pi, \quad (C8)$$

where the dependence of the kernel on the hyperparameters is indicated by $\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{UV}}, \nu$, presented in Sec. V C 1. $\sigma$ is the overall magnitude of the kernel, while $l_\omega$ and $l_k$ are the RBF lengthscales in $\omega$ and $k$ direction respectively. $l_{\text{UV}}$ is the steepness of the asymptotic transition, while $\nu$ controls the position of the transition in Eq. (C6) with $\omega_0 = \sqrt{k^2 + \nu(\pi T)^2}$, in analogy to the polynomial fit ansatz.

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 \to \sigma/k^2$. In Fig. 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 \ge 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{UV}} \ge 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.
<figure><img src="image_19.png" /><figcaption>FIG. 16. Metropolis-Hastings sampled GP kernel parameters in the negative log-likelihood plane.</figcaption></figure>

[1] A. Bazavov *et al.* (HotQCD Collaboration), Phys. Lett. B **795**, 15 (2019).

[2] A. Bazavov *et al.* (HotQCD Collaboration), Phys. Rev. D **90**, 094503 (2014).

[3] S. Borsányi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, and K. K. Szabó, Phys. Lett. B **730**, 99 (2014).

[4] R. D. Pisarski, Phys. Lett. **110B**, 155 (1982).

[5] L. D. McLerran and T. Toimela, Phys. Rev. D **31**, 545 (1985).

[6] K. Kajantie and H. I. Miettinen, Z. Phys. C **9**, 341 (1981).

[7] K. Kajantie, J. Kapusta, L. McLerran, and A. Mekjian, Phys. Rev. D **34**, 2746 (1986).

[8] C. Gale and J. Kapusta, Phys. Rev. C **35**, 2107 (1987).

[9] S. Turbide, R. Rapp, and C. Gale, Phys. Rev. C **69**, 014903 (2004).

[10] A. Adare *et al.* (PHENIX Collaboration), Phys. Rev. C **91**, 064904 (2015).

[11] A. Adare *et al.* (PHENIX Collaboration), Phys. Rev. Lett. **104**, 132301 (2010).

[12] M. Wilde (ALICE Collaboration), Nucl. Phys. **A904–905**, 573c (2013).

[13] A. Adare *et al.* (PHENIX Collaboration), Phys. Rev. C **94**, 064901 (2016).

[14] D. Lohner and the ALICE Collaboration, J. Phys. Conf. Ser. **446**, 012028 (2013).

[15] R. Chatterjee, E. S. Frodermann, U. W. Heinz, and D. K. Srivastava, Phys. Rev. Lett. **96**, 202302 (2006).

[16] H. van Hees, C. Gale, and R. Rapp, Phys. Rev. C **84**, 054906 (2011).

[17] C. Shen, U. W. Heinz, J.-F. Paquet, and C. Gale, Phys. Rev. C **89**, 044910 (2014).

[18] J. Churchill, L. Du, C. Gale, G. Jackson, and S. Jeon, Phys. Rev. C **109**, 044915 (2024).

[19] J. Churchill, L. Du, C. Gale, G. Jackson, and S. Jeon, Phys. Rev. Lett. **132**, 172301 (2024).

[20] A. K. Chaudhuri and B. Sinha, Phys. Rev. C **83**, 034905 (2011).

[21] J.-F. Paquet, C. Shen, G. S. Denicol, M. Luzum, B. Schenke, S. Jeon, and C. Gale, Phys. Rev. C **93**, 044906 (2016).

[22] G. Vujanovic, J.-F. Paquet, G. S. Denicol, M. Luzum, S. Jeon, and C. Gale, Phys. Rev. C **94**, 014904 (2016).

[23] R. Baier, H. Nakkagawa, A. Niégawa, and K. Redlich, Z. Phys. C **53**, 433 (1992).

[24] J. Kapusta, P. Lichard, and D. Seibert, Phys. Rev. D **44**, 2774 (1991); **47**, 4171(E) (1993).

[25] P. Arnold, G. D. Moore, and L. G. Yaffe, J. High Energy Phys. **11** (2001) 057.

[26] P. Arnold, G. D. Moore, and L. G. Yaffe, J. High Energy Phys. **12** (2001) 009.

[27] J. Ghiglieri, J. Hong, A. Kurkela, E. Lu, G. D. Moore, and D. Teaney, J. High Energy Phys. **05** (2013) 010.

[28] G. Cuniberti, E. De Micheli, and G. A. Viano, Commun. Math. Phys. **216**, 59 (2001).

[29] H. B. Meyer, Eur. Phys. J. A **47**, 86 (2011).

[30] J. Ghiglieri, O. Kaczmarek, M. Laine, and F. Meyer, Phys. Rev. D **94**, 016005 (2016).

[31] M. Cè, T. Harris, A. Krasniqi, H. B. Meyer, and C. Török, Phys. Rev. D **106**, 054501 (2022).

[32] M. Cè, T. Harris, H. B. Meyer, A. Steinberg, and A. Toniato, Phys. Rev. D **102**, 091501 (2020).

[33] H. B. Meyer, Eur. Phys. J. A **54**, 192 (2018).

[34] M. Cè, T. Harris, A. Krasniqi, H. B. Meyer, and C. Török, Phys. Rev. D **109**, 014507 (2024).

[35] D.-J. Zhao, G. Wang, F. He, L. Jin, P. Sun, Y.-B. Yang, and K. Zhang (χQCD Collaboration), Phys. Rev. D **107**, L091501 (2023).

[36] P. A. Baikov, K. G. Chetyrkin, J. H. Kühn, and J. Rittinger, Phys. Lett. B **714**, 62 (2012).

[37] L. D. McLerran and T. Toimela, Phys. Rev. D **31**, 545 (1985).

[38] S. Caron-Huot, Phys. Rev. D **79**, 125009 (2009).

[39] B. B. Brandt, A. Francis, T. Harris, H. B. Meyer, and A. Steinberg, EPJ Web Conf. **175**, 07044 (2018).

[40] R. Sommer, Proc. Sci. LATTICE2013 (2014) 015 [arXiv: 1401.3270].

[41] A. Francis, O. Kaczmarek, M. Laine, T. Neuhaus, and H. Ohno, Phys. Rev. D **91**, 096002 (2015).
[42] A. Bazavov et al. (MILC Collaboration), Proc. Sci. LATTICE2010 (2010) 074 [arXiv:1012.0868].

[43] L. Altenkort, D. de la Cruz, O. Kaczmarek, R. Larsen, G. D. Moore, S. Mukherjee, P. Petreczky, H.-T. Shu, and S. Stendebach (HotQCD Collaboration), Phys. Rev. Lett. 132, 051902 (2024).

[44] M. Lüscher, S. Sint, R. Sommer, P. Weisz, and U. Wolff, Nucl. Phys. B491, 323 (1997).

[45] E. Follana, Q. Mason, C. Davies, K. Hornbostel, G. P. Lepage, J. Shigemitsu, H. Trottier, and K. Wong (HPQCD and UKQCD Collaborations), Phys. Rev. D 75, 054502 (2007).

[46] M. Lüscher and P. Weisz, Commun. Math. Phys. 97, 59 (1985); 98, 433(E) (1985).

[47] M. Lüscher and P. Weisz, Phys. Lett. 158B, 250 (1985).

[48] M. A. Clark and A. D. Kennedy, Nucl. Phys. B, Proc. Suppl. 129, 850 (2004).

[49] M. Guagnelli and R. Sommer, Nucl. Phys. B, Proc. Suppl. 63, 886 (1998).

[50] H.-T. Ding, O. Kaczmarek, and F. Meyer, Phys. Rev. D 94, 034504 (2016).

[51] G. Jackson and M. Laine, J. High Energy Phys. 11 (2019) 144.

[52] M. Laine, J. High Energy Phys. 05 (2013) 083.

[53] G. Jackson, Phys. Rev. D 100, 116019 (2019).

[54] P. Aurenche, F. Gelis, G. D. Moore, and H. Zaraket, J. High Energy Phys. 12 (2002) 006.

[55] R. Baier, B. Pire, and D. Schiff, Phys. Rev. D 38, 2814 (1988).

[56] Y. Gabellini, T. Grandou, and D. Poizat, Ann. Phys. (N.Y.) 202, 436 (1990).

[57] T. Altherr and P. Aurenche, Z. Phys. C 45, 99 (1989).

[58] E. Braaten, R. D. Pisarski, and T.-C. Yuan, Phys. Rev. Lett. 64, 2242 (1990).

[59] G. D. Moore and J.-M. Robert, arXiv:hep-ph/0607172.

[60] M. Laine and A. Rothkopf, J. High Energy Phys. 07 (2013) 082.

[61] S. Caron-Huot, Phys. Rev. D 79, 065039 (2009).

[62] S. Aoki et al. (Flavour Lattice Averaging Group), Eur. Phys. J. C 80, 113 (2020).

[63] A. Vuorinen, Phys. Rev. D 67, 074032 (2003).

[64] S. Borsányi, Z. Fodor, S. D. Katz, S. Krieg, C. Ratti, and K. Szabó, J. High Energy Phys. 01 (2012) 138.

[65] S. Mogliacci, J. O. Andersen, M. Strickland, N. Su, and A. Vuorinen, J. High Energy Phys. 12 (2013) 055.

[66] S. Gupta, Phys. Lett. B 597, 57 (2004).

[67] D. Bala, S. Ali, A. Francis, G. Jackson, O. Kaczmarek, and T. Ueding, Proc. Sci. LATTICE2022 (2023) 169 [arXiv:2212.11509].

[68] P. K. Kovtun and A. O. Starinets, Phys. Rev. D 72, 086009 (2005).

[69] G. Backus and F. Gilbert, Geophys. J. Int. 16, 169 (1968).

[70] N. Yu. Astrakhantsev, V. V. Braguta, and A. Yu. Kotov, Phys. Rev. D 98, 054515 (2018).

[71] M. Kanagawa, P. Hennig, D. Sejdinovic, and B. K. Sriperumbudur, arXiv:1807.02582.

[72] C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning (Springer, New York, 2006), Vol. 1.

[73] A. P. Valentine and M. Sambridge, Geophys. J. Int. 220, 1632 (2020).

[74] 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).

[75] J. M. Pawlowski, C. S. Schneider, J. Turnwald, J. M. Urban, and N. Wink, Phys. Rev. D 108, 076018 (2023).

[76] J. Horak, J. M. Pawlowski, J. Turnwald, J. M. Urban, N. Wink, and S. Zafeiropoulos, Phys. Rev. D 107, 076019 (2023).

[77] I. Steinwart, J. Mach. Learn. Res. 2, 67 (2001).

[78] C. A. Micchelli, Y. Xu, and H. Zhang, J. Mach. Learn. Res. 7, 2651 (2006).

[79] S. Ali, D. Bala, A. Francis, G. Jackson, O. Kaczmarek, J. Turnwald, T. Ueding, and N. Wink, Data publication for “Lattice QCD estimates of thermal photon production from the QGP”, 10.4119/unibi/2992848.

[80] L. Mazur et al. (HotQCD Collaboration), Comput. Phys. Commun. 300, 109164 (2024).