An estimate for the thermal photon rate from lattice QCD

We estimate the production rate of photons by the quark-gluon plasma in lattice QCD. We propose a new correlation function which provides better control over the systematic uncertainty in estimating the photon production rate at photon momenta in the range {\pi}T/2 to 2{\pi}T. The relevant Euclidean vector current correlation functions are computed with $N_{\mathrm f}$ = 2 Wilson clover fermions in the chirally-symmetric phase. In order to estimate the photon rate, an ill-posed problem for the vector-channel spectral function must be regularized. We use both a direct model for the spectral function and a model-independent estimate from the Backus-Gilbert method to give an estimate for the photon rate.


Introduction
Lattice QCD has had a tremendous impact in computing the equilibrium properties of strongly interacting matter at non-zero temperature. Its real-time properties, however, are difficult to compute from numerical lattice calculations because they require the analytic continuation of incomplete and noisy Euclidean two-point functions. Yet, the dynamical parameters of thermal systems, like transport coefficients and production rates of weakly-coupled external currents, are needed to interpret, for example, the equilibration of the thermal medium formed after a heavy-ion collision [1] or the abundance of particles in the early universe [2]. In addition, they are interesting in their own right to characterize the medium.
Speaker, e-mail: harris@him.uni-mainz.de Speaker, e-mail: amsteinb@uni-mainz.de 1 Here we assume that the number of temporal indices among µ, ν is not exactly one, in which case the kernel is different.
It represents an ill-posed inverse problem when the correlator is finitely sampled, even though definite algorithms for this task exist [3]. Nevertheless, there may exist particular observables which can be constrained by the Euclidean observables.
Here we propose a new correlator whose spectral function has no UV growth, from which the production rate of real photons from a thermal medium of quarks and gluons can be computed. Furthermore, we present preliminary results of the photon rate using a continuum-extrapolated correlator computed with N f = 2 flavours of O(a)-improved Wilson fermions.

The photon production rate
The differential production rate of weakly-coupled massless vector bosons from a thermal medium of gluons and quarks in the flavour-symmetric limit is given, to leading order in the electromagnetic coupling α em , by [4] where n B (ω) = 1/[exp(βω) − 1] is the Bose-Einstein distribution. The rate is therefore determined by the spectral function evaluated on the photon mass-shell. This rate can be directly compared to the measurement of direct non-prompt real photons, i.e. those not originating from hadronic decays or initial partonic scatterings, which have been observed at RHIC [5] and the LHC [6,7]. Theoretically, this rate has been computed in QCD at weak-coupling [8,9], and QCD-like theories at strong-coupling [10]. These two scenarios show qualitatively different behaviour in the soft limit k ≡ |k| → 0, where the weak-coupling prediction for ρ µ µ (k, k)/k is of order 1/α 2 s , while the strong-coupling prediction is independent of the coupling. It would therefore be desirable to have a lattice prediction in the region of photon momenta of a few times the temperature to distinguish these possibilities, and further elucidate the nature of the medium itself. Pioneering numerical calculations have been performed in pure gauge theory [11], which utilized advanced perturbative results on the UV behaviour of the Euclidean correlator G µ µ (τ, k), closely-related to the dilepton production rate, to constrain the photon rate.

UV-finite correlator
As a consequence of the continuum vector Ward identity ω 2 ρ 00 = k i k j ρ i j , on the light cone It is therefore helpful to consider the linear combination The spectral function with λ = 1 corresponds to ρ 1 (ω, k) = ρ µ µ (ω, k), directly relevant to the dilepton rate. However, this channel contains the large ultraviolet contribution of order ω 2 at large ω; at small k, it also contains the diffusion pole. To avoid confronting these complications simultaneously, one can analyze separately 1. the purely transverse case λ = 0, which does not couple to the diffusion pole; 2. the case λ = −2, corresponding to the difference of the transverse and longitudinal channels.
The latter linear combination vanishes identically in the vacuum and is highly suppressed in the ultraviolet. Here we concentrate on the case λ = −2; in the future, we plan to also analyze the case λ = 0, which should yield consistent photon rates, thus providing a powerful cross-check. At the end, the spectral functions with λ = 0 and λ = −2 can be recombined in order to predict the dilepton rate. The importance of removing UV divergences from Euclidean correlators to estimate thermal real-time observables has also been discussed in ref. [12].   Figure 1 illustrates the effect of the cancellation on the tree-level spectral function in the solid lines. The standard correlation function (λ = 1) is shown for the lowest momenta in the dashed line, which diverges as ω 2 at large frequencies. The spectral function with λ = −2 on the other hand is very suppressed for ω > k, thus making this channel very sensitive to the infrared physics of interest. Note that the spectral function evaluated on the photon mass-shell (at the kink), and thus the photon rate, vanishes at this order in perturbation theory. If one thinks of the inverse problem as resulting in a 'smearing' of the actual spectral function, as is explicitly the case in the Backus-Gilbert method, then this represents a difficulty, since the spectral weight is of order unity for ω k.

Continuum limit
We have generated a series of ensembles to take the continuum limit at a single temperature, approximately T = 250 MeV, above the crossover to the chirally symmetric phase, and an additional ensemble at a single lattice spacing deep in the deconfined phase, approximately T = 500 MeV; see table 1. We use the non-perturbatively O(a)-improved Wilson action [13] with N f = 2 Wilson fermions and the Wilson gauge action. The parameters were chosen using the running of the coupling and quark masses as determined by the CLS collaboration [14]. The lattice with N τ ≡ β/a = 16 at T ≈ 250 MeV, where β = T −1 is the inverse temperature, has been used for our previous studies, see refs. [15][16][17].
In order to control the continuum limit we measured the two-point correlation functions of the vector current using both local and exactly-conserved discretizations of the current. Furthermore, in the case of the local-conserved correlation function, there are two discretizations of the λ = −2 linear combination, where in first case the conserved current is constructed on the lattice site and in the second case on the midpoint of the link. Therefore, we can define the λ = −2 correlator at half-integer or integer units of the lattice spacing, through the following replacements This gives in total four independent discretizations of the correlation function associated with the spectral function of eq. (4): local-local (LL), local-conserved site-centered (LC site), local-conserved link-centered (LC link) and conserved-conserved (CC). The correlation functions were measured with all on-and off-axis spatial momenta up to kβ 2π. The continuum limit was performed by making a piecewise spline interpolation of the correlation functions and performing a quadratic regression in the lattice spacing. This is motivated by the fact that in the chirally-symmetric phase, the matrix elements of the improvement operators are all suppressed by the quark mass in units of the temperature, and O(a)-effects are almost absent. This behaviour is also clearly observed in the data, as shown in figure 2 (left panel) for a given spatial momentum and Euclidean time separation. The continuum limit of each of the discretizations agrees, which shows that discretization effects should be under control, thus making a simultaneous extrapolation viable.   Figure 2. The simultaneous continuum limit of all four discretizations of the correlator at a fixed Euclidean distance, momentum and temperature. The continuum limit without and with tree-level improvement is shown in the left and right panels respectively. Furthermore, we employ an alternative improvement prescription by multiplicatively removing the tree-level lattice artifacts, via The effect of this improvement is shown in figure 2 (right panel). Although this provides a substantial correction to the finite lattice spacing data, reassuringly it does not alter the continuum limit. The continuum-extrapolated correlator is displayed as the open black symbols in figure 3, which also shows the data away from the continuum with the closed coloured symbols. The extrapolation is constrained by the data in both panels and also the other two discretizations not shown. Except at short Euclidean distances τ β/4, the β/a = 24 correlator lies close to the continuum-extrapolated values. The addition of this large and fine lattice is therefore very valuable for the analysis of cutoff effects.

Analysis of the photon rate
We employ two methods to regularize the inverse problem in the finite-dimensional analogue of eq. (1) in order to analyze the continuum-extrapolated correlators in terms of the spectral function. The first is a linear method which gives a model-independent estimator for a smeared spectral function, and the second is to fit the parameters of an explicit model to the data. These are complementary approaches because exact constraints that the spectral function must satisfy can be built into the model and checked a posteriori for the linear method.

Backus-Gilbert method
The Backus-Gilbert (BG) method [18] is defined by constructing a linear map from the space of functions in the time domain, G, to the space of functions on the frequency domain, ρ BG , by choosing the factor, n q n (ω)K(na, ω) ≡δ(ω, ω), called the resolution function, to be an optimally localized kernel in the frequency domain. In the case of a constant spectral function the BG estimator coincides with the spectral function, which explains the advantage of a slowly varying spectral function. Other kernels can be constructed with a variety of properties, at the loss of resolution in the frequency domain.
In figure 4 (left), the resolution functions are shown at a selection of values of the first argument, ω. Due to the finite number of data, in practice the resolution functions have support on the whole frequency domain, and for this temperature, are far from being close to Dirac distributions. This means that the estimator ρ BG (ω) picks up contributions from the true spectral function at nearby frequencies, and so is smeared out. Generally, the resolution decreases at larger values ofω.
In the right-hand panel of figure 4, an extra linear constraint is introduced,δ(ω, 0) = 0, which means that the estimator obtained from this method does not contain contributions from the true spectral function at ω = 0. At small k, hydrodynamics predicts a zero mode in the charge-charge correlator which could give a large contribution to the BG estimator in the region of the origin. We use this method to suppress this contribution to the photon rate, and to check the influence of this diffusion pole on our results. For larger frequencies,ω, this constraint is irrelevant, but for the lowest frequency shown,ωβ = 2.4, the support of the resolution function is pushed away from the origin.  in ref. [11], where χ s = d 4 x V em 0 (x)V em 0 (0) , which we know in the continuum as well. The light-cone frequency is indicated with the vertical solid bar. The estimator from the BG method with the extra constraint δ(ω, 0) = 0 is shown with the dashed blue curves. At the lower photon momentum, there is a systematic difference between the two estimators for the effective diffusion constant. This might be expected from increasing contributions from the diffusion pole in the hydrodynamic limit. On the other hand, for larger photon momenta, this effect is negligible, as seen in the right-hand panel of the figure.
In order to obtain an estimate for the systematic uncertainty on the effective diffusion constant from the BG method, we perfomed many variations of the method which give rise to systematic errors, including variations of the minimum time separation of the correlator and using the extra constraint. All the variations are listed in table 2. In figure 6, the distributions are shown for two photon momenta, kβ = π, 3π/2 in the left and right panels respectively. For the lower momenta (left), implementing the extra constraint results in a bimodal distribution, and the two estimators are clearly incomptabile. In the following the results from this variation are kept separate. At larger momenta (right), the results become compatible which indicates that contributions from the diffusion pole are absent in the estimator for the effective diffusion constant, and therefore combined to obtain the final estimate of the sytematic uncertainty. As a final estimator, we take the median of the distribution of the effective diffusion constant under all variations, and as a systematic uncertainty the 68% interval.

Model Fit
In the previous subsection 3.1, a model-independent approach is discussed, the Backus-Gilbert method. To complement this approach we study a physics-inspired ansatz for the tanh-regulated spectral function in this subsection. Thus, we are able to perform a maximum likelihood estimation by a nonlinear fit. Before discussing the ansatz we first prove some analytic properties of the spectral function which will motivate and strongly constrain the choice of simple models we can use for a maximum likelihood estimator.

Padé ansatz
The tanh-regulated spectral function can be described by a combination of two Padé approximants as (14) with two linear parameters A and B as well as three nonlinear parameters a, ω 0 , and b. Part I is inspired by the diffusion pole as it arises in the hydrodynamics prediction for the infrared limit; when one identifies a ↔ Dk 2 for small k, part I resembles the known expression [10] Part II models the pole structure of the AdS/CFT current correlator, see ref. [19] for details. In order to satisfy the superconvergent sum rule (13), a second linearly independent parameter B has to be introduced in part II of eq. (14). By imposing eq. (13), the second linear parameter B becomes a function of (a, ω 0 , b). The two poles at (±ω 0 , −b) in the lower half-plane do not only model the quasinormal modes of the retarded correlator, discussed in ref. [19], they also match the 1/ω 4 behaviour at large ω as dictated by the operator product expansion (11), see also figure 7.

Uncorrelated fit
In the fit ansatz (14), we are left with four independent fit parameters after determining B for given (a, ω 0 , b) via the sum rule (13). Then there are three degrees of freedom when we use seven data points between τ min /β = 0.25 and τ max /β = 0.5 on the continuum-extrapolated correlators. There is a huge subspace in the parameter space (A, a, ω 0 , b) for which the uncorrelated χ 2 is smaller than one. Figure 8 depicts the χ 2 -valley in the (ω 0 , b)-plane at photon momentum k/T ≈ 4.97 where χ 2 < 1 at two fixed values of a.
Because there is no local minimum in the uncorrelated χ 2 -landscape, this results in a plethora of acceptable solutions. In other words, there are many shapes of the spectral function, none of which can be strongly ruled out by our data as they all satisfy χ 2 (A, a, ω 0 , b) < 1. So rather than minimizing the uncorrelated χ 2 , we try and find bounds to the effective diffusion constant by taking the min and max values of all photon rates with χ 2 (A, a, ω 0 We exclude all solutions that result in a negative photon rate from this procedure. This is due to the known positivity of the spectral function below the light cone, When one allows the second pole at (ω 0 , b) to get too close to the real axis, the result is a very pronounced peak in the spectral function. This can happen when the nonlinear parameter a becomes large can approach the real axis. In figure 9, the result is seen as a sharp peak below the light cone (dashed blue curve). For a/T = 5 (right panel of figure 8), however, the second pole at (ω 0 , b) is pushed into the complex plane and does not dominate the shape of the spectral function.
Such an additional peak from a pole close to the real axis corresponds to a very long-lived excitation in the medium which is unphysical if it survives longer than the largest possible relaxation times in the system. So we constrain the imaginary parts of both poles, i.e. the fit parameters a and b, to fulfil the exclusion criterion [10,20] where D AdS/CFT · k 2 describes the diffusion of an electric charge with D AdS/CFT = 1/(2πT ) and D −1 PT accounts for the damping of a static current with D −1 PT = O(α 2 s ) · T , α s = 0.25. We claim that the exclusion criterion (17) amounts to a conservative constraint based on physics considerations. Figure 5 compares the BG estimators with and without the constraint on the resolution function with the model ansatz (14) for two different photon momenta kβ = π and 3π/2 in the left and right panel, respectively. For kβ = π the biggest discrepancy between the BG method and Padé ansatz occurs below the light cone for the BG estimator with the constraint. On the light cone and above, the discrepancy is small regardless of whether the resolution function is constrained or not. For kβ = 3π/2 there is almost no systematic difference between the BG estimators with or without a constraint and the fit ansatz. An important point to note is that on the light cone there is very good agreement of the BG estimators and the median of the model distribution in all cases.

Results
Having discussed and compared the two approaches in the previous section 3, we can now plot the effective diffusion constant, eq. (9), as a function of photon momentum k. Figure 10 shows the results both from the BG method and the spread from the distribution of solutions to the uncorrelated fit at T = 250 MeV. As discussed in subsection 3.1, the two BG estimators resulting from implementing or not implementing the constraint are not compatible with each other at lower momenta and the two results are plotted separately (green vs. purple dots). At higher momenta they become compatible and only one distribution is considered to estimate a systematic and statistical error. The black bars indicate the bounds we find from quoting the minimal and maximal value of the distribution of the effective diffusion constant that have χ 2 < 1 and the median of this distribution. Furthermore, the strong-coupling result from N = 4 supersymmetric Yang-Mills theory (SYM) and a weak-coupling result from leading-order (LO) perturbative QCD with α s = 0.25 are plotted as comparison.
For photon momenta kβ 3.5 the median of the uncorrelated fit distribution and the BG estimator are remarkably close and agree within the systematic and statistical error of the BG estimator. At the lower momenta the median of the distribution lies between the two BG estimators with and without the constraint. For the whole momentum range, the bounds of the effective diffusion constant from examining the minimum and maximum values of the maximum likelihood estimator consistent with the data, cover a big interval and it is not possible to discriminate between the weak-coupling and the strong-coupling scenarios. At lower momenta the spread from the model is compatible with the separation of the two BG estimators. Figure 11 shows the estimate of the effective diffusion constant deep in the deconfined phase at T ≈ 500 MeV. For this temperature, we only have one ensemble at a single lattice spacing so there is no continuum extrapolation available, and the Padé ansatz has not yet been analyzed for this data. By comparison with figure 10, however, we do not observe any strong temperature dependence of this observable.

Summary and Outlook
We presented an estimate of the photon rate from dynamical QCD based on continuum-extrapolated correlators. In order to be more sensitive to the physics on the light cone, we propose an alternative linear combination of the vector-vector correlator that eliminates the UV contamination.
We avail of two qualitatively different approaches to estimate the photon rate: the Backus-Gilbert method is a linear mapping that reconstructs a smeared estimator of the true spectral function; and a Padé fit ansatz serves as a model inspired by relativistic hydrodynamics, AdS/CFT and plausibility constraints. We can exploit the UV behaviour of the spectral function expected from an OPE and derive a superconvergent sum rule that is implemented into the model.
As the uncorrelated χ 2 -landscape is rather degenerate, we quote the median of the distribution of acceptable solutions with χ 2 < 1 and the min and max values of this distribution. The median coincides with the BG estimator at larger momenta and is compatible with the separation of the two BG estimators at lower momenta. In order to gain more control in the small k region, it can be useful to define another effective diffusion coefficient D eff (ξ, k) = ξρ λ=−2 (ξk,k) 4χ s k → D which tends to the true diffusion constant D for k → 0 at fixed ξ ∈ [0, 1).
A very similar strategy can be applied to energy-momentum tensor correlators. Let ρ µν,λσ (ω, k) be the spectral function associated with the T µν T λσ correlator. With the goal in mind to compute the shear viscosity η, consider, for k = (0, 0, k), the linear combination [ρ 12,12 − ρ 13,13 + ρ 10,10 ]. It vanishes in the vacuum, is positive for ω 2 − k 2 < 0 and falls off as k 2 /ω 2 at asymptotic frequencies. It is therefore a promising starting point for a calculation of the shear viscosity η.
Apart from increasing statistics it might be beneficial to further examine our systematics including the continuum limit. One should also take into account the correlations between data points as the exclusionary power of the correlated χ 2 is stronger. Further cross-checks on the BG estimator can be investigated, such as verifying how well the sum rule is satisfied. Finally, extending the study to a higher temperature will illuminate the temperature dependence of the photon rate.