Photon production spectrum above $T_c$ with a lattice quark propagator

The photon production rate from the deconfined medium is analyzed with the photon self-energy constructed from the quark propagator obtained by the numerical simulation on the quenched lattice for two values of temperature, $T=1.5T_{\rm c}$ and $3T_{\rm c}$, above the critical temperature $T_{\rm c}$. The photon self-energy is calculated by the Schwinger-Dyson equation with the lattice quark propagator and a vertex function determined so as to satisfy the Ward-Takahashi identity. The obtained photon production rate exhibits a similar behavior as the perturbative results at the energy of photons larger than $0.5$~GeV.


I. INTRODUCTION
Photons are important observables in relativistic heavy-ion collisions [1,2]. Because of its colorless nature, the mean free path of photons is longer than the size of the hot matter produced in heavy-ion collisions. Once photons are produced, therefore, they penetrate the hot medium without scatterings via the strong interaction. The experimental measurements of photon yield and anisotropic flow have been carried out actively [3,4].
The experimental results show that the photon yield from heavy-ion collisions has an excess at lower transverse momentum regions [3,4] compared with the reference yield obtained by the superposition of proton-proton collisions. This result suggests the importance of the photon radiation from the hot medium. The other important observable is the photon anisotropic flow. The magnitude of the observed photon flows is almost the same as that of the medium [5]. The large anisotropic flow implies that photons are dominantly emitted from late stages in which the flow of the medium has been well developed. Phenomenological analyses of photon yield based on dynamical models have not succeeded in reproducing these experimental results simultaneously so far [6,7]. Photons are produced via several mechanisms during the time evolution of heavy-ion collisions. The discrepancy between the phenomenological analysis and experimental results suggests the importance of careful surveys of the production rates through these different production mechanisms. In the present study, among these mechanisms we focus on the photon radiation from the hot deconfined medium.
In the previous phenomenological analyses on the photon yield and flows [6,7], the photon production rate calculated in perturbation theory [8][9][10][11][12] has been used for those from deconfined medium. On the other hand, the success of hydrodynamic models in describing the time evolution of the hot medium suggests that the deconfined medium is a strongly coupled system [13,14], in which a perturbative analysis would lose its validity. A non-perturbative analysis thus is desirable to evaluate the production rate from the deconfined medium. However, the direct analysis of the emission rate in lattice QCD numerical simulations is quite difficult, because this analysis requires an analysis of on-shell spectral functions.
In this paper, we analyze the photon production rate from the deconfined medium by partially using the lattice data in order to incorporate the non-perturbative effects in the analysis. The photon production rate is related to the photon self-energy in medium [15][16][17]. From the Schwinger-Dyson equation (SDE), the full photon self-energy is constructed from the full quark propagator and the full photon-quark vertex. In Refs. [18][19][20], the quark propagator in the deconfined medium is analyzed on the quenched lattice in Landau gauge with a pole ansatz. We employ this lattice propagator as the full quark propagator in the SDE. A photon-quark vertex is constructed from the lattice quark propagator so as to satisfy the Ward-Takahashi identity. Our analysis thus respects gauge invariance. The same analysis has been applied to the calculation of the dilepton production rate at zero three-momentum in our previous study [21]. The present study extends this analysis to nonzero momentum to calculate the photon production rate.
We analyze the photon production rate for two temperatures, T = 1.5T c and 3T c , where T c is the pseudo critical temperature. We show that the photon production is dominated by the pair annihilation and Landau damping of two quasi-particle, the normal and plasmino, modes. Although the pair annihilation cannot produce real photons in the perturbative analysis, in our formalism peculiar dispersion relations of quasi-particle modes in the lattice quark propagator allow the process. We find that the obtained photon production rate exhibits a similar behavior as the perturbative results both in the slope and the magnitude at the photon energy larger than 0.5 GeV, although the photon production mechanism is completely different.
The outline of this paper is as follows. In the next section, we introduce the SDE for the photon self-energy, and discuss the construction of its components, the lattice quark propagator and a photon-quark vertex satisfying the Ward-Takahashi identity. In Sec. III, we calculate the photon self-energy with this setup with and without vertex correction. In Sec. IV, we show the numerical results on the photon production spectrum. The final section is devoted to a summary.

II. SCHWINGER-DYSON EQUATION FOR PHOTON SELF-ENERGY
A. Schwinger-Dyson equation The photon production rate from a static medium per unit time per unit volume is related to the retarded photon self-energy Π R µν (ω, q) as [15][16][17] with the inverse temperature β = 1/T and the fourmomentum of photon (ω, q). This formula is correct at the leading order in electromagnetic interaction and all orders in strong interaction. The retarded photon self-energy is related to the imaginary-time self-energy Π µν (iω m , q) by the analytic continuation, with the positive infinitesimal number η and the Matsubara frequency for bosons ω m = 2πT m with integer m. The imaginary-time full photon self-energy is given by the SDE in Matsubara formalism with the full quark propagator S(P ) and the full photon-quark vertex Γ µ (P + Q, P ) as where P µ = (iν n , p) and Q µ = (iω m , q) are the fourmomenta of quarks and photons, respectively, and ν n = (2n + 1)πT with integer n is the Matsubara frequency for fermions. e f is the electric charge of quarks with the index f representing the quark flavor. The color, flavor, and Dirac indices of S(P ) are suppressed for notational Diagrammatic representation of the Schwinger-Dyson equation for the full photon self-energy, Eq. (3). The shaded circles mean the full quark propagator and the full photon-quark vertex.
simplicity, and Tr C and Tr D represent the trace over the color and Dirac indices, respectively. In covariant gauges, the off-diagonal elements of the quark propagator in color space vanish. As a result, the trace over the color indices gives a factor three in Eq. (3). In Fig. 1, we show the diagrammatic representation of Eq. (3). The shaded circles mean the full propagator and vertex function. In the following, we consider the two-flavor QCD with degenerate u and d quarks, in which f e 2 f = 5e 2 /9.

B. Lattice quark propagator and spectral function
In the present study, we employ the quark propagator obtained in the lattice numerical simulation in Ref. [18][19][20] for the full quark propagator in Eq. (3). This approach was used in our previous study for the analysis of the dilepton production rate at zero threemomentum [21]. In this subsection, we give a brief review on the general property of the quark propagator and discuss the implementation of the results in Ref. [20] in our analysis.
The imaginary-time quark propagator for momentum p is defined by where ψ µ (τ, x) is the quark field with the Dirac index µ, and τ is the imaginary time restricted to the interval 0 ≤ τ < β. We show the Dirac indices of quark fields explicitly for a moment. The Fourier transform of Eq. (4) is given by Equation (5) is rewritten with the quark spectral function ρ µν (ν, p) as From Eqs. (5) and (6), the spectral function is related to Eq. (4) as e (1/2−τ /β)βν e βν/2 + e −βν/2 ρ µν (ν, p) .
In Refs. [18][19][20], the imaginary-time quark propagator Eq. (4) is evaluated on the lattice in Landau gauge with the quenched approximation. On the lattice the quark propagator is measured for discrete imaginary times. To obtain the real time information, one has to deduce the spectral function from the discrete information. In Refs. [18][19][20], the spectral function ρ ± (ν, p) is analyzed with the two-pole ansatz, where Z ± (p) and ν ± (p) are the residues and positions of the poles, respectively, which are determined by fitting to the lattice propagator for each p.
Comments on the two-pole ansatz Eq. (12) are in order. First, two particle-like states ν + (p) and ν − (p) in Eq. (12) correspond to the normal and plasmino modes, respectively, which appear in the quark propagator in the hard thermal loop approximation [25]. In fact, the analyses in Refs. [18][19][20] suggest that these modes behave like the normal and plasmino modes as functions of momentum p and bare quark mass m 0 . For example, for large p and m 0 the residue of the plasmino mode Z − (p) becomes small and the propagator approaches the one of the free quark. Second, the restoration of chiral symmetry for massless quarks above the pseudo-critical temperature T c is confirmed explicitly on the lattice by measuring the scalar term of the quark propagator [19,20]. Third, in Refs. [19,20] various fitting ansätze have been tested besides the two-pole ansatz Eq. (12), such as an ansatz which allows non-zero widths of the poles in Eq. (12). However, it was turned out that fitting with widths does not improve χ 2 , as the minimum of the χ 2 is at vanishing widths. From these results, we suppose that the existence of sharp quasi-particle states in the quark spectral function is supported from the lattice analysis. Fourth, the quark propagator depends on the gauge in general, Open symbols show the momentum dependence of the parameters ν+(p), ν−(p), and Z−(p)/(Z+(p) + Z−(p)) obtained on the lattice in Ref. [20]. The solid lines represent their interpolation obtained by the cubic spline method. The dashed line in the upper panel represents the light cone.
while the analyses on the lattice have been performed in Landau gauge. We, however, note that the positions of poles in the quark propagator are independent of the gauge [22].
In Fig. 2, we show the p dependences of the fitting parameters in Eq. (12) obtained in Ref. [20] for T = 1.5T c and 3T c . The upper panel shows the dispersions relation ν + (p) and ν − (p), while the lower panel shows the relative strength of the residue of the plasmino mode Z − (p)/(Z + (p) + Z − (p)). These analyses are performed on lattices with volume 128 3 × 16, where both the lattice spacing and finite volume effects are found to be small [20]. The two axes in the upper panel are normalized by the thermal mass m T defined by the energy of the quasi-particle excitations at p = 0 as m T = (ν + (0) + ν − (0))/2. After the extrapolation to infinite volume, m T /T is estimated as 0.768(11) and 0.725 (14) at T = 1.5T c and 3T c , respectively [20].
Although the lattice data are available only for discrete values of p because of the periodic boundary condition, one needs the quark propagator as a continuous function of p to solve Eq. (3). To obtain the fitting parameters as continuous functions of p, we perform the interpolation of the lattice data by the cubic spline method. From Eq. (11), one can show that dν + (p)/dp = −dν − (p)/dp, d 2 ν + (p)/dp 2 = d 2 ν − (p)/dp 2 , and Z + (p) = Z − (p) at p = 0 [19,20]. These properties are taken into ac-count in the interpolation. The lattice data are available only in the momentum range p/T 4.7, and we need the extrapolation to higher momentum. We extrapolate the parameters using an exponentially damping form for and for ν ± (p), which approaches the light cone exponentially for large p. The parameters R, α, and β ± i are determined in the cubic spline analysis. The p dependence of each parameter determined in this way is shown by the solid lines in Fig. 2. We tested another extrapolation form, ν ± (p) = p + β ± 1 /p + β ± 2 /p 2 + · · · , but found that the dispersion relations hardly change. Finally, we set to satisfy the sum rule Eq. (10). We note that the slope of the plasmino dispersion relation exceeds unity for p 3m T and is acausal. This unphysical behavior may come from an artifact of the lattice simulation and/or the ansatz for the spectral function. This acausal behavior gives rise to the Cherenkov radiation, as we will see in later sections. However, as we show explicitly in Sec. IV, the residue of the plasmino mode is small, Z − (p) < 0.05, in this momentum range and the contribution of the Cherenkov radiation to our final result is negligibly small except for the low energy region ω 0.1 GeV.
Using Eqs. (6), (8), and (12), the quark propagator is given by with where we define S ±1 (P ) = S ± (P ) and Λ ±1 (p) = Λ ± (p) for notational simplicity. Correspondingly, the inverse quark propagator is given by with and . Note that the inverse propagator has poles at iν n = ±V (p). These poles inevitably appear in the multipole ansatz, because the propagator S s (ω, p) after the analytic continuation to real time has one zero point in the range of ω surrounded by two poles. We will see that the poles at iν n = ±V (p) give rise to additional photon production mechanism that we call "anomalous" productions.

C. Photon-Quark Vertex
To solve the SDE, Eq. (3), one needs the full photonquark vertex Γ µ (P + Q, P ) besides the full quark propagator. The vertex Γ µ (P + Q, P ) must satisfy the Ward-Takahashi identity (WTI) One, however, cannot determine the four components of the vertex uniquely from this constraint. In the present study, as the full vertex satisfying the WTI, we choose the following simple representation: One can easily check that Eqs. (21) and (22) satisfy the WTI Eq. (20). In this choice, we construct the vertex so that Γ µ becomes the bare vertex γ µ when the free quark propagator S −1 (P ) = p is substituted. Although our vertex Eqs. (21) and (22) is not the unique choice [24], we emphasize that our analysis with Eqs. (21) and (22) satisfies the gauge invariance required by the WTI.

III. PHOTON PRODUCTION RATE
Next, we construct the photon self-energy Eq. (3) from the lattice quark propagator and the gauge invariant photon-quark vertex obtained in the previous section. In this section, we analyze the photon self-energy for two cases. First, we calculate the photon production rate with the lattice quark propagator but with the bare vertex in Sec. III A. This result will help us understand the effect of the vertex correction. Then, we perform the full analysis in Sec. III B.

A. Photon Production Rate without Vertex Correction
In this subsection, we calculate the photon self-energy without the vertex correction with Γ µ = γ µ . The photon self-energy in this case is given by By substituting the decomposition Eq. (16) into Eq. (23), one obtains with being the momenta of the quarks, and q = |q|, and we used Tr In Eq. (25), we changed the three-momentum integral as where the integral is carried out in the region R, which is surrounded by the lines p 1 + p 2 = q and p 1 − p 2 = ±q as shown in Fig. 3. The region R is defined by the p 1 and p 2 which can satisfy |p 2 − p 1 | = |q|.
Taking the analytic continuation iω m → ω + iη of Eq. (25) and the imaginary part, one obtains the imaginary part of the retarded self-energy, Owing to the δ-function, the momentum integral in Eq. (28) is reduced to the one-dimensional one. In Fig. 4, we show an example of the path determined by the δfunction by the red line for ω = |q| = 5m T and s = t = −η 1 = η 2 = −1 at T = 1.5T c ; the line in this case satisfies ν + (p 1 ) + ν − (p 2 ) = 5m T . In our numerical analysis, after determining the line in R numerically we carry out the line integral.
Equation (28) consists of 16 terms with s, t, η 1 , η 2 = ±1. These terms are interpreted as the photon production rates from the decay and scatterings of two quasiquarks as follows. First, the δ-function in each term expresses the energy conservation, and the residues Z ± represent those of the corresponding quasi-particles. Next, the thermal factor f (sη 1 ν η1 (p 1 )) − f (tη 2 ν η2 (p 2 )) is transformed as −[1−f (ν 1 )−f (ν 2 )] or f (ν 1 )−f (ν 2 ) for a given set of (s, t, η 1 , η 2 ) with ν 1 = ν ± (p 1 ) and ν 2 = ν ± (p 2 ). When the thermal factor becomes −[1 − f (ν 1 ) − f (ν 2 )], this factor can be rewritten as Because this is the difference between the products of thermal distributions and Pauli blocking effects, this term is interpreted as the difference of the photon emission from the pair annihilation of quasi-quarks shown in Fig. 5(a) and the photon absorption mediated by its inverse process. Similarly, we have The photon emission and absorption of this term is diagrammatically represented as in Fig. 5(b), which is called the Landau damping [23].
In the perturbative calculation, the pair annihilation with two quasi-particles cannot produce a real photon [8][9][10][11]. This is because the quasi-particle states are in the time like region and these processes cannot fulfill the onshell condition in the processes in Fig. 5. On the other hand, as in Fig. 2 the plasmino mode in the lattice quark propagator enters space-like region. Because of this behavior, it is possible that the pair annihilation produces a real photon directly. This is the most remarkable difference of our analysis from the perturbative one.
Note that we are interested in the photon production rate, while ImΠ R,µ µ represents the difference of the production and absorption rates as shown above. To obtain the former, as in Eq. (1) one multiplies the Bose distribution function 1/[e βω − 1] to ImΠ R,µ µ (ω, q). Indeed, we have

B. Full Photon Production Rate
Now, we consider the photon self-energy with the full vertex Eqs. (21) and (22).
The last piece i ImΠ Finally, by taking the trace of ImΠ R µν (ω, q), we obtain the full photon self-energy as [f (−tV (p 2 )) − f (tη 2 ν η2 (p 2 ))]δ(ω − tV (p 2 ) − tη 2 ν η2 (p 2 )) − sη 1 Zη 1 (p 1 )ν(p 1 ) tη 2 ν η2 (p 2 ) + sV (p 1 ) [f (sη 1 ν η1 (p 1 )) − f (−sV (p 1 ))]δ(ω + sη 1 ν η1 (p 1 ) + sV (p 1 )) The photon production rate is obtained by substituting this result into Eq. (1). Let us inspect the physical meaning of each term in Eq. (42). Similar to the discussion in Sec. III A, the terms in the last two lines in Eq. (42) are identified to be the pair annihilation and Landau damping of quasi-particles from the structures of the δ-function and thermal factor. These terms have different coefficients from those in Eq. (28) owing to the vertex correction. The terms in the first three lines in Eq. (42) seem to be photon emissions mediated by the quasi-particles having energies ν η and V . However, the lattice quark propagator does not have the quasi-particle mode with the dispersion ν = V (p). Therefore, the photon productions expressed in these terms cannot be interpreted as simple reactions between quasi-particles and a photon. These photon productions come from the anomalous pole included in the inverse propagator S −1 s in the vertex function Eqs. (21) and (22) [21]. Although these terms are peculiar at a glance, they naturally appear to satisfy the gauge invariance and thus must be included in the full photon production rate.

IV. NUMERICAL RESULTS
Next, let us see the numerical results on the photon production spectrum. In Fig. 6, we show the photon production rate with the lattice quark propagator and the gauge invariant vertex as a function of the photon energy ω for T = 1.5T c by the solid line. The production rate with the bare vertex obtained in Sec. III A is also plotted by the dotted line. The photon production rate has been calculated in the perturbative analysis [10][11][12]. In the figure, we show fitting formula for the perturbative result in Refs. [10,11] for three values of the strong coupling constant α s = 0.1, 0.2, and 0.3 by the dashed lines; this fitting function well reproduces the complete leading log result for 0.2 ≤ ω/T ≤ 50 [11]. The next-to-leading order result [12] enhances the rate at most 20% around ω/T ∼ 2.
From Fig. 6, one finds that the production rates with full and bare vertices behave similarly for ω 0.2 GeV. This result suggests that the effect of the vertex correction is not large. The figure also shows that the production rate with the lattice quark propagator behaves quite similarly with the perturbative result with α s 0.1 − 0.2 except for low energy region. This result is unexpected, because the production mechanisms of photons are different between these analyses.
As discussed in the previous section, the photon pro- duction rate can be decomposed into contributions of different processes. In Fig. 7, we show the contributions of each term in Eq. (42) on the production rate. In the figure, the lines show the rate from the pair annihilation of a normal and a plasmino modes (PA(ν + ν − )) and two plasmino modes (PA(ν − ν − )), the Landau damping between a normal and a plasmino modes (LD(ν + ν − )) and two plasmino modes (LD(ν − ν − )), the rates of two anomalous productions (Anom.(V ν + ) and Anom.(V ν − )), and the total (full) rate. The figure shows that the photon production for ω 1.5 GeV is dominated by the Landau damping LD(ν + ν − ), while for higher energies two production mechanisms LD(ν + ν − ) and PA(ν + ν − ) have almost the same contributions. The process LD(ν − ν − ) is interpreted as the Cherenkov radiation from the acausal dispersion relation of the plasmino mode. From Fig. 7, one finds that the photon production with this process is well suppressed except for ω 0.2 GeV. This result shows that the effect of this unphysical process is negligible. The figure also shows that the contribution of the anomalous productions (Anom.(V ν + ) and Anom.(V ν − )) are small compared to other processes except for the low energy region ω 0.2 GeV. The decomposition of the photon production rate obtained with the bare vertex, Eq. (28), shows the similar tendency as the one of the full result. Figure 7 shows that the processes PA(ν + ν − ) and PA(ν − ν − ) take nonzero values for ω > ω +− and ω −− , respectively, with ω +− 0.86 GeV and ω −− 0.63 GeV. Owing to the emergence of these processes, the total production rate has non-smooth behaviors at these energies.
To understand the reason why these processes emerge at ω = ω +− and ω −− , in Fig. 8 we show the minimum of the photon energy produced by these processes with a fixed photon momentum q. In the figure, the minimum energy for the pair annihilation of two normal modes (ν + ν + ) is also plotted for comparison. In the regions above the lines, virtual photon production is allowed kinematically. From the figure, one finds that the region covers even a part of the space-like region for ν + ν − and ν − ν − , and thus these processes can produce a real photon, while the region for ν + ν + does not have an overlap with the light cone. The real photon production from two quasi-particle pair annihilation does not occur in the perturbative calculation based on the hard thermal loop resummation technique, because the normal and plasmino modes in this formalism exists only in the time-like region. The lattice propagator, on the other hand, has the plasmino mode which enters space-like region as in Fig. 2. This behavior is responsible for the emergence of the nonzero production rates by PA(ν + ν − ) and PA(ν − ν − ). The photon spectra shown in Fig. 6 has a cusp structure at ω = ω max with ω max 0.2 GeV. According to Fig. 7, this cusp comes from the Landau damping LD(ν + ν − ). The origin of this structure is understood from the kinematics between quasi-particle pair and photon. Let us denote the momenta of the normal and plasmino modes p + and p − , respectively, and the photon momentum q. Then, these momenta satisfy p + = p − + q.
For ω < ω max , the decay process is allowed for all directions of quark momenta. At ω = ω max , the process with perpendicular p + and q, i.e. the case with |p + | + |q| = |p − | starts to be forbidden kinematically, and the phase space becomes narrower as ω = |q| becomes larger. This constraint in the phase space gives rise to the cusp in the production spectrum. It is a surprising result that the photon production spectrum in our analysis has a similar behavior, not only in the slope but also the absolute value, as the perturbative results for almost all energy range. The slopes in both results are approximately given by the Boltzmann factor e −βω . The origin of this slope in our analysis can be understood as follows. From Fig. 7, the photon production is dominated by two processes, PA(ν + ν − ) and LD(ν + ν − ). These processes include the plasmino mode whose residue Z − is exponentially damping for large momentum. Because of this factor, the momentum of the plasmino mode must be small, p − 0. Then, to produce a real photon with large momentum |q| m T , the momentum of the normal mode must satisfy p + = q ± p − q, and the photon emission takes place via the decay of this mode. The rate of this process is proportional to the thermal probability of the existence of the normal mode given by e −βν+(q) e −βq = e −βω . For sufficiently large q, this factor determines the production rate, because other effects, such as the phase space determined by kinematics and Z + (|p + |) 1, are insensitive to q. In the perturbative analysis, the photon self-energy at high energy is dominated by two-loop diagrams and collinear enhanced diagrams [11]. In these processes, an initial quark with large momentum q must decay into soft degrees of freedom, and the Boltzmann factor appears for the same reason as ours, while it is known that other factors are insensitive to q [11]. To summarize, the coincidence in the slope comes from the thermal factor for the initial quark with large momentum. Because the decay processes in these analyses are completely different, however, the coincidence in magnitude would be accidental.
In Fig. 9, we show the result for T = 3T c . Its behavior is qualitatively the same as the production rate for T = 1.5T c , although the spectrum has a clearer dip structure at low energy ω 0.5 GeV. The contributions of different decay processes behave similarly, too.
The results at both temperatures have a similar form and magnitude with the perturbative results. If a closer look is taken, however, we recognize that our nonperturbative results are almost reproduced by perturbative results with α s 0.1 − 0.2, which is rather on the lower side. This has an important phenomenological consequence; the discrepancy between the observed photon yield at RHIC and the corresponding theoretical result is not resolved, but is enlarged. Thus, it is confirmed also by this investigation that there will exist still unknown mechanisms for photon production in relativistic heavy ion collisions.

V. SUMMARY
In this paper, we analyzed the production spectrum of thermal photons from the deconfined medium near the critical temperature T c . In order to include the non-perturbative effect in the analysis, we analyzed the Schwinger-Dyson equation of the photon self-energy with the quark propagator obtained on the lattice for two temperatures. The photon-quark vertex is constructed so as to satisfy the WTI. The effect of the vertex correction is discussed. We showed that the vertex function satisfying the WTI generates anomalous photon production mechanisms but its contribution is small. In comparison with the perturbative analysis, the spectra with the lattice quark propagator show almost the same behavior for ω 0.5 GeV. We discussed that the coincidence in the slope can be understood from the production mechanism. However, it is a nontrivial result that the our result is similar to the perturbative one even in the magnitude. The gap between experimental data and theoretical calculations on the yield and flows of photons still remain. Our result indicates the existence of unknown mechanisms of photon production in relativistic heavy ion collisions.
The authors thank T. Shimoda, K. Ogata, and K. Oda for encouragements. T. K. is supported by Grant-in-Aid for JSPS research fellows (No. 16J00981) and Osaka University Cross-Boundary Innovation Program. This work is supported in part by JSPS KAKENHI Grant Numbers 26400272 and 17K05442.