The rate of photon production in the quark-gluon plasma from lattice QCD

We calculate the thermal rate of real-photon production in the quark-gluon plasma at a temperature of $T=254$ MeV using lattice QCD. The calculation is based on the difference between the spatially transverse and longitudinal parts of the polarization tensor, which has the advantage of falling off rapidly at large frequencies. We obtain this linear combination in the time-momentum representation from lattice QCD with two flavors of quarks in the continuum limit with a precision of about two parts per mille. Applying a theoretically motivated fit ansatz for the associated spectral function, we obtain values for the photon rate that are in line with QCD weak-coupling calculations; for photon momenta $ 1.0\leq k[{\rm GeV}]\leq 1.4$, our non-perturbative results constrain the rate to be no larger than twice the weak-coupling prediction. We also provide a physics interpretation of the electromagnetic spectral functions valid for all frequencies and momenta.


I. INTRODUCTION
Strongly interacting matter undergoes a phase transition at a temperature of about 150 MeV [1][2][3]. Below the transition, the thermal medium is characterized by hadrons (nucleons, pions, kaons, . . . ) as primary degrees of freedom, while well above the transition it is characterized by quarks and gluons, the elementary degrees of freedom of quantum chromodynamics (QCD). The high-temperature phase, the quark-gluon plasma (QGP), is probed experimentally in high-energy heavy-ion collisions at T 500 MeV [4]. One of the remarkable properties of the medium is its ability to exhibit collective effects in spite of the rapid expansion occurring in heavyion collisions. The most prominent such effect is the large anisotropic flow observed in heavy-ion collisions at RHIC and the LHC, pointing to a small shear viscosity to entropy density ratio of the medium; see e.g. [5] and Refs. therein. In addition, probes of the medium that do not interact strongly are of great interest, since they escape largely unscathed once produced. In particular, the rate at which photons are emitted by the QGP is a classic -though challenging -observable in heavy-ion experiments. Direct photons with a transverse momentum below 2 GeV are found to admit an exponential spectrum, and models assuming the formation of the QGP are consistent with these measurements [6,7]. The production of weakly interacting particles by the QGP is also an important issue in early-universe cosmology, for instance in models which propose a keV-scale sterile neutrino as a dark matter candidate [8,9].
In this Letter we address the rate of photon emission from the QGP via lattice QCD simulations. One motivation to perform the calculation is that the rate vanishes in the limit of non-interacting quarks and gluons; therefore it is a measure of the strength of their interactions.
Secondly, direct photons emitted in heavy-ion collisions have been found to exhibit an unexpectedly large central value of elliptic flow [10,11] -albeit with significant uncertainty, therefore addressing their thermal production rate non-perturbatively can contribute to resolving the issue. Thirdly, a controlled calculation of the photon rate paves the way for calculating the production of other particles, such as lepton pairs -relevant in heavyion phenomenology -or sterile neutrinos -relevant for validating or ruling out a dark matter candidate.
The main computational difficulty stems from the production of weakly-interacting particles being a real-time process, which is accessible from the Matsubara path integral formalism implemented in lattice QCD only via an analytic continuation [12]. Numerically, the latter amounts to a poorly conditioned inverse problem discussed below.

II. THEORY BACKGROUND
We consider the full set of spectral functions of the electromagnetic current 1 (1) For any four-vector u µ , the form u † µ ρ µν (ω, k)u ν /ω is real and non-negative; for u µ real, it is also even in ω. Current conservation leads to ω 2 ρ 00 (ω, k) = k i k j ρ ij (ω, k), implying that 2 (k ikj ρ ij − ρ 00 )/ω has the same sign as K 2 ≡ ω 2 − k 2 , and that it vanishes at lightlike kinematics, K 2 = 0. It will be useful to consider the linear combination Defining the Euclidean correlator 3 the corresponding linear combination admits the spectral representation The production rate of dileptons with invariant masssquared equal to K 2 , which occurs via a timelike photon, is proportional to ρ(ω, k, 1) [13]. To leading order in the fine-structure constant α = e 2 /(4π), the differential photon rate per unit volume of plasma can be written as and does not depend on λ. The Euclidean correlator G(x 0 , k, λ) probes the spectral function for all virtualities K 2 ≥ −k 2 . It is therefore desirable to have an interpretation of the spectral function for negative virtualities. The cross-section per unit volume for an electron scattering on the medium through the exchange of a spacelike photon is given by µν ≡ 2(p µ p ν + p ν p µ − g µν (p · p )), with p and p respectively the initial and final electron momenta and k = p − p . Eq. (7) refers to the restframe of the thermal medium. More generally, the vector spectral functions can be interpreted as the ability of the medium to dissipate the energy stored in electromagnetic fields: consider coupling the plasma to a harmonic external vector potential A(t, x) = Re( A k e i( k· x−ωt) ), by adding the term ∆H = −e d 3 x j · A to the Hamiltonian. The energy of the external electromagnetic fields is given by E e.m. = 1 . A fraction of this energy gets transferred to the medium per unit time and turned 2 We use the notation k ≡ | k| into heat. We find, for the transverse and longitudinal cases, the following rates of energy transfer, dE e.m. dt = e 2k ikj ρ ij (ω, k) ω = e 2 ω k 2 ρ 00 (ω, k).
These equations provide an interpretation of the spectral functions for all virtualities. The positivity of the spectral functions on the right-hand side guarantees that the medium obeys the second law of thermodynamics. Given the goal of computing the photon rate, computationally its non-dependence on the value of the parameter λ can be exploited to one's advantage. We choose λ = −2, because as a combined consequence of current conservation and Lorentz invariance, ρ(ω, k, −2) vanishes identically in the vacuum (at zero temperature). Due to the latter property and because ρ(ω, k = 0, −2) vanishes exactly for ω = 0 due to charge conservation, we expect from the operator-product expansion This strong suppression in the ultraviolet implies a superconvergent sum rule for ρ(ω, k, −2), Spectral positivity implies that ρ(ω, k, −2)/ω is nonnegative for K 2 < 0, and it must become negative for K 2 > 0 in order to satisfy the sum rule (10). There are two regimes in which the functional form of the spectral function is known. In the infrared limit, the ρ 00 contribution parametrically dominates ρ(ω, k, −2) and the hydrodynamic prediction is where D is the diffusion coefficient and χ s ≡ β G 00 (x 0 , 0) the static susceptibility. Therefore, following [14], we define the effective diffusion coefficient which is proportional to the photon rate and tends to D in the limit k → 0. In the weak-coupling regime, results at order g 2 have recently become available for general (ω, k) [15,16]. The photon rate itself has been obtained at order g 3 in [17]. From here on we set λ = −2 and omit the last argument of ρ(ω, k, λ) and G(ω, k, λ).

III. THE LATTICE CALCULATION
We use lattice QCD with an isospin doublet of O(a) improved Wilson fermions at a temperature of T = 254 MeV; the details of the lattice action can be found in [20] and references therein. Table I lists our ensembles, which allow us to take the continuum limit at a fixed temperature. All but the finest ensemble have a renormalized quark mass of m MS 13 MeV in the MS scheme at a renormalization scale of µ = 2 GeV; on the finest ensemble, we have m MS 16 MeV. Quark-mass effects, which are suppressed by (m/T ) 2 in the chirally symmetric phase, are therefore expected to be negligible. The ensembles F7, O7 and X7 were generated using the MP-HMC algorithm [21] in the implementation described in Ref. [22] based on the DD-HMC package [23], while ensemble W7 was generated using twisted-mass Hasenbusch frequency splitting in the version 1.6 of open-QCD [24,25]. The ensembles labelled F7 and O7 have bare parameters identical to the zero-temperature F7 and O7 ensembles described in [20], for which the pion mass is 269 MeV.
We compute the correlator G(x 0 , k) of the isovector current 1 √ 2ψ γ µ τ 3 ψ, which consists of a single connected Wick contraction 4 . The corresponding static susceptibility amounts to G(x 0 , 0)/(2T 3 ) = χ s /T 2 = 0.880(9) stat (8) syst in the continuum limit, where the systematic error reflects the dependence on using different prescriptions for the renormalisation of the local vector current. We employ the local and the conserved vector currents, resulting in four discretizations of G(x 0 , k), and perform a constrained simultaneous continuum extrapolation. We have computed the leading-order perturbative lattice predictions, so that we are able to correct for the corresponding cutoff effects affecting our Monte-Carlo data. To avoid incurring large cutoff effects at short distances, we omit data points for x 0 < x min 0 , where x min 0 = β/4 is our default value. We thus have data points for 7,8,9,10,11,12}. Given 4 In order to keep the notation concise, we do not explicitly distinguish between the quantities derived from the isovector and from the electromagnetic current. To obtain the photon rate from our results for D eff (k), we recommend using Eq. the high accuracy of the data, we are led to leave out the ensemble with the coarsest lattice spacing from the continuum limit. Figure 1 illustrates the correlator obtained at different lattice spacings and its continuum limit. The relative statistical precision of the continuum correlator is one to two permille. It is well-known that the topological charge Q acquires a long autocorrelation time at small lattice spacings, and our simulations confirm this effect. However, we have found the dependence of the vector correlator of interest on |Q| to be at most at the 3% level. Therefore the vector correlator only suffers a modest increase in uncertainty from this algorithmic difficulty. We define the observable Expressed in terms of the spectral function, in the limit x 0 → β/2 it describes the ratio of the ω 2 moment to the ω 0 moment of ρ(ω, k)/ sinh(ωβ/2). The 1/k 2 factor allows for a finite k → 0 limit. It is instructive to compare the results from lattice QCD with the theory of noninteracting quarks as well as with an extreme opposite, namely the N = 4 super-Yang-Mills (SYM) theory in the limit of infinite 't Hooft coupling and infinite number of colors; the spectral functions of the latter are obtained via the AdS/CFT correspondence [26]. One qualitative difference between the spectral functions of the strongly coupled SYM theory and of free quarks is that in the former case the positive spectral weight of the spacelike region ω 2 < k 2 'leaks' into the timelike region; see especially the second panel of times larger in the former theory. It is thus interesting to ask how R(x 0 , k) behaves in QCD at the temperature of 254 MeV. The observable is displayed in Fig. 2. The QCD values lie less than 20% above the non-interacting values.

IV. ANALYSIS OF THE SPECTRAL FUNCTION
To obtain a global picture of the spectral function without committing to any specific functional form, in [27] we applied the Backus-Gilbert method to our data. The results confirm the theoretical expectation that most of the spectral weight is contained in the spacelike region ω 2 < k 2 .
A second method [27], which we now pursue further, consists in applying an explicit fit ansatz for the spectral function, The ansatz satisfies the expected large-ω behavior (9). We always determine the parameter B in terms of (ω 0 , a, b) by imposing the sum rule (10) and require B ≥ −1/k 2 to satisfy the spectral positivity condition for ω 2 < k 2 . Thus, for a single momentum k, Eq. (14) amounts to a four-parameter fit. The Euclidean correlator resulting from the spectral function (14) can be expressed as a linear combination of Lerch transcendents Φ(e ±2πix0 , 1, 1 2 + i ωp 2π ), where ω p are the frequency poles of ρ(ω, k)/ tanh(ωβ/2).
We impose the following physically motivated constraints on the parameters. Spectral positivity implies  6) and (12). The color-coded vertical bars represent those values of D eff for which a spectral function of the form (14) exists that has a p-value above 0.32.
The colors indicate the smallest χ 2 /d.o.f. found for a given value of D eff . Shaded areas identify the momentum groups that are fitted simultaneously; for each momentum, results are shown both for the γ = 1 and γ = 2 parametrizations of the k-dependence of the nonlinear parameters. Analytical results from perturbative QCD [28] and from the strong-coupling limit of N = 4 super-Yang-Mills theory [26] are shown for comparison.
A ≥ 0 and B ≥ −1/k 2 . Furthermore, since there cannot be arbitrarily long relaxation times in the system, we impose the condition on the poles, where D strong = 1 2πT is the diffusion coefficient of the strongly coupled SYM theory and D −1 weak the inverse QCD diffusion coefficient at leading-order in the perturbative expansion, which we set to 0.46T based on the results of [29]. This condition reflects the fact that Dk 2 is the rate of dissipation of a perturbation in the charge density, while D −1 provides an estimate of the relaxation rate of a homogeneous current.
In order to increase the discriminative power of our fits, we simultaneously fit data at different momenta. The correlators have been computed for all spatial momenta k = π T 2 ν for ν ∈ Z 3 and n ≡ | ν| 2 ≤ 16. We found it convenient to split the set of available momenta into three groups, 1 ≤ n ≤ 3, 3 ≤ n ≤ 8 and 8 ≤ n ≤ 14, which contain respectively N k = 3, 5 and 7 momentum values. The number of data points entering a fit is thus given by N k N t , the number N t of Euclidean times being seven in our data set. Within a momentum group, we parameterize the momentum dependence of ansatz (14) by expressing the nonlinear parameters a, b and ω 0 as functions of the momentum. We consider two polynomial forms in our analysis, with γ = 1 or 2 and where k min is the smallest momentum in the group.
Since the covariance matrix C of the data points is sizeable, we have used the regularized matrixC, constructed according to 5 We have studied the stability of our results with respect to the regularization parameters x, y and found little dependence on them around the values we chose [30]. For instance, both x and y were set to 0.95 for the second momentum group. The importance of preserving correlations among the input data points when addressing the inverse problem has been emphasized previously [31].
For each momentum group, we performed a scan in the six-dimensional space of non-linear parameters (a 0 ,a 2 ,b 0 ,b 2 ,W 0 ,W 2 ), while, at each momentum, the parameter B is determined by imposing the sum rule (10) and the linear parameter A by minimizing the χ 2 . The number of fit parameters is thus given by 6 + N k , and the number of degrees of freedom for each of the three momentum groups is 12, 24 and 36 respectively. We calculate the p-value of each set of parameter values and consider that it provides a satisfactory description of the correlator whenever p > 0.32. If the condition is satisfied, the corresponding D eff (k) are marked as being compatible with the lattice data, and the associated p-value is recorded.
Before describing our results for D eff (k), we briefly present the outcome of our procedure when applied to mock Euclidean data generated from known spectral functions. For these tests, we have used the spectral functions of non-interacting quarks as well as those of the strongly coupled SYM theory. In order to be realistic, we re-use the covariance matrix of our lattice QCD data, rescaled so as to achieve the same relative error on the correlator. In both cases, we find that the correct value of D eff (k) is one of those having a p-value above 0.32. The output spectral functions yielding the highest p-value tend to have a somewhat larger value of D eff (k).
Our final results for the D eff (k) values yielding a pvalue above 0.32 for the QCD correlator at T = 254 MeV are displayed in Fig. 3. We show results for both the linear and the quadratic dependence on k, γ = 1 and 2. We observe that for the third momentum group, containing momenta above 1.0 GeV, the values of D eff (k)·GeV cover the interval [0, 0.7] and are thus compatible both with the leading-order weak-coupling prediction [28] and the strongly-coupled SYM prediction [26], which lie between 0.3 and 0.5. Moreover, the weak-coupling prediction is among those values with the highest p-value. In the second momentum group, the range of acceptable D eff (k) values covers a range up to about twice the stronglycoupled SYM value (for the ansatz quadratic in k), while again the weak-coupling prediction has one of the highest p-values. In the group of smallest momenta, the lattice data loses sensitivity to the photon rate. Particularly, the data does not exclude large values of D eff (k). Finally, we remark that our fits yield a strong correlation between the values of D eff (k) at successive k [30].
It is instructive to look at the full frequency dependence of the spectral functions which describe the QCD correlators. In Fig. 4, as representative examples for the three spatial momenta k = (0.40, 0.98, 1.49) GeV, we show spectral functions that correspond to the upper and lower end of the D eff (k) ranges shown in Fig. 3. We also display the spectral function leading to the smallest χ 2 , and for comparison, the spectral functions of noninteracting quarks as well as those of the strongly coupled SYM theory. For the second and third momenta, we observe that all three spectral functions describing the QCD correlators exhibit a smooth behaviour for ω 2 < k 2 and admit a maximum near the point ω = k, its precise location being tightly linked to the value of D eff (k) and hence to the photon emission rate.

V. CONCLUSION
Using lattice simulations in the quark-gluon plasma phase of QCD with two dynamical quark flavors at a temperature of 254 MeV, we have computed one particularly ultraviolet-soft component of the polarization tensor in the continuum limit. This component determines the photon emission rate from the medium via analytic continuation, in practice however one is faced with the inverse problem (Eq. (5)) for the spectral function. We explored exhaustively the parameter space of the Padéform spectral functions in Eq. (14). The photon rate is given, up to kinematical factors, by the spectral function at photon kinematics, ω = k, and normalizing this quantity by the well-determined static charge susceptibility, one obtains the effective (momentum-dependent) diffusion coefficient. Within the explored family of spectral functions, we determined which values of this coefficient are compatible with the Euclidean data; our result is displayed in Fig. 3. We have validated our handling of the inverse problem by applying the same procedure to two field theories that represent extreme opposite caricatures of the quark-gluon plasma. Our results imply non-perturbative constraints on the possible rate of pho-ton emission from the QGP at a temperature typical for the strongly interacting system created in heavy-ion collision experiments. We largely confirm the weak-coupling predictions, in spite of the relatively low temperature of 254 MeV. Our results are also in good agreement with those of a previous lattice calculation performed in the quenched approximation [14].
As a study based on correlators in the theory of noninteracting quarks shows, adding data points at shorter Euclidean time significantly enhances the ability of the data to exclude large values of D eff , particularly at low photon momenta. This calls for even finer lattices to be used, which represents a challenge in view of the large lattice sizes required and the long associated autocorrelation times.
As described in [27], an analogous combination of correlators can be applied to energy-momentum tensor correlators to extract an effective shear viscosity η eff (k). We finally remark that a different strategy has also recently been proposed to compute the photon rate using a dispersion relation at fixed, vanishing virtuality [32]. The ultraviolet-soft channel employed here also plays an important role in the implementation of this alternative method.