Low Mass Thermal Dilepton Production at NLO in a Weakly Coupled Quark-Gluon Plasma

We present a computation, within weakly-coupled thermal QCD, of the production rate of low invariant mass ($M^2 \sim g^2 T^2$) dileptons, at next-to-leading order (NLO) in the coupling (which is $O(g^3 e^2 T^2)$). This involves extending the NLO calculation of the photon rate which we recently presented to the case of small nonzero photon invariant mass. Numerical results are discussed and tabulated forms and code are provided for inclusion in hydrodynamical models. We find that NLO corrections can increase the dilepton rate by up to 30-40% relative to leading order. We find that the electromagnetic response of the plasma for real photons and for small invariant mass but high energy dilepton pairs (e.g., $M^2<(300\:\mathrm{MeV})^2$ but $p_T>1 \: \mathrm{GeV}$) are close enough that dilepton pair measurements really can serve as Ersatz photon measurements. We also present a matching a la Ghisoiu and Laine between our results and results at larger invariant masses.


Introduction
Most of the particles originating from an ultra-relativistic heavy ion collision are believed to undergo significant "rescattering" (parton energy loss, scattering, fragmentation, hadronization, hadron-hadron scattering) before they emerge and fly to the detector. This degrades the information they carry about the early stages of the collision, which must be reconstructed by trying to model the scattering processes, for instance with hydrodynamics. Hard probes, on the other hand, are anything sufficiently high-energy or weakly-coupled that they can penetrate the heavy ion environment with little re-interaction. They therefore carry more direct and unprocessed information about early conditions. Electromagnetic probes are a good example, because the electromagnetic coupling strength α = 1/137 is small enough that re-interactions are rare.
Photons are one hard probe candidate. With this in mind we recently performed an improved determination of the rate of photon production from a (weakly-coupled) quarkgluon plasma [1]. But photons come with some experimental challenges. The highest energy photons are expected to come primarily from collisions between the initial state partons [2], and therefore carry little information about the QGP evolution. At intermediate energies where the Quark-Gluon Plasma is expected to contribute (say, 1 to 4 GeV), a photon excess has been observed [2]. These measurements are challenging because of the large background of decay photons which must be subtracted. For this reason, experimentalists have also focused on small-mass dileptons, which we can think of as massive off-shell photons. Provided the mass-squared of the pair is above the pion mass, the (Dalitz) pion decay background is absent and the foreground rates are under much better control. For this reason, e + e − pairs with invariant masses somewhat above m 2 π have been measured, to serve as an ersatz photon rate measurement [3].
The dilepton production rate is indeed related to the real photon production rate, but they are not quite the same, and we need some theory input to understand their relation. At lowest order in the electromagnetic coupling α and in equilibrium, 1 both the rate per unit 4-volume to produce a photon and to produce a dilepton are determined by the current-current correlation function with J µ = eQψγ µ ψ the electromagnetic current operator and with .. an average over the quark-gluon plasma ensemble. In terms of Π < , the photon production rate per unit phase space is where the photon 4-momentum K is taken on-shell as indicated. The dilepton rate exists for any timelike positive-energy K µ such as that −K 2 > 4m 2 l (with m l the lepton mass), and away from the threshold it is given by (see for instance Ref. [4]) (1. 3) The difference, besides the factor of α/3πK 2 , is that the current-current correlator Π < (K) is evaluated at a timelike (massive) value for the dilepton rate, and at a null (massless) value for physical photons. So how much does the shift from null to timelike 4-momentum change things? Fig. 1 shows Π < (K) for a fixed |k| = 6T and a range of k 0 values. The solid (black) curve is a free-theory level calculation in QCD; the dotted (red) curve is the value in strongly interacting N =4 SYM theory, from [5] (normalized to have the same large-k 0 behavior). In free (zero-coupling) QCD there is a cusp at the real-photon point. At strong coupling, Π < is smooth. If physical QCD behaves more like the red curve, then the photon rate 1 In general for a nonequilibrium system, one finds the total (not per-unit-volume) rate by replacing Π < (K) with its space-integrated version, W < (K) ≡ d 4 Xd 4 Y e iK·(Y −X) J µ (Y )Jµ(X) . The relationship between photon and dilepton production involves whether k 0 = |k| or k 0 > |k| in the same way. We will not consider the nonequilibrium case because we have so few tools for the calculation of fully nonequilibrium correlation functions. Figure 1. The current-current correlator Π < (K) for the case k = 6T , as a function of k 0 (normalized by α and its dominant k 0 dependence). The black curve is for free (g = 0) QCD, illustrating the cusp at k 0 = k. The red dotted curve shows the behavior in N =4 SYM theory at infinite coupling, the most strongly coupled QCD-like theory known, which shows no cusp-like behavior. and the small-mass-squared dilepton rate are almost interchangeable. If it behaves more like the black curve, then the dilepton rate will show a sharp dependence on the invariant mass of the dilepton, and photon production will be suppressed relative to expectations based on moderate invariant-mass dileptons (if those expectations are based on Eq. (1.2) and Eq. (1.3) and the assumption of smooth behavior in Π < ).
The goal of this paper is to provide the most complete perturbative calculation of Π < for K µ close to lightlike which is currently possible. Previously, Ref. [6] have shown how to compute the dilepton rate for K 2 parametrically in the range K 2 ∼ g 2 T 2 at leading order in the coupling. We improve this determination to the next order in the strong coupling g. We also extend the result to larger virtuality, K 2 ∼ gT 2 , and discuss the matching onto the recently completed next-to-leading order calculation at large invariant mass squared [7,8]. Our main motivation is to improve Fig. 1, showing how the finite-coupling, perturbative rate behaves near the real-photon point K 2 = 0.
Besides the phenomenological justification we have presented, there is an additional theoretical reason to be interested in doing this. It is possible to determine the Euclideantime-domain behavior of Π nonperturbatively on the lattice [9,10]. At least in principle, this can be analytically continued to determine the real-frequency behavior which is physically interesting, for instance, by applying an Ansatz [10] or using the Maximum Entropy Method [11]. Unfortunately, in practice this method is very bad at reconstructing frequency-domain functions which possess sharp features, such as that displayed by the black curve in Fig. 1. This is particularly so if the feature is not expected and is not built into the model function (priors) used in the reconstruction. Therefore, determining whether we expect such a feature would be very useful in characterizing and improving the reliability of lattice reconstructions for the photon and dilepton rates.
The outline of the remainder of the paper is as follows. We begin in Section 2 by reviewing the leading-order calculation of Π < (K) for real photons (K 2 = 0). We then show the two extensions in the existing literature: the extension to small virtuality, and the extension to next-to-leading order (NLO). Then, Section 3 shows how to apply both extensions at once, giving the NLO dilepton rate at small virtuality. Next, Section 4, presents the results in two ways. We show the results in the strict small-g limit. And we make phenomenological plots at finite values of g. To do so, we discuss how our calculation connects onto the recent calculation of Laine and Ghisoiu [7,8], which is valid to O(g 2 ) corrections for −K 2 ∼ T 2 . We end with a discussion, two technical appendices, and a third appendix collecting tabulated numerical results.
Very briefly, the two most important take-home results of our work are, that the dilepton rate at very small invariant mass is increased by up to ∼ 30% relative to the leading order one, and by a smaller amount at larger invariant mass. And we find that high-energy but low-mass dilepton pairs, such as those recently studied by the PHENIX collaboration [3], are probing essentially the same electromagnetic response as real photons are; so their use as ersatz photon measurements is valid.

Review of Previous Results
The real photon rate was first calculated at leading order by Arnold, Moore, and Yaffe (AMY) in Ref. [12,13]. Aurenche et al showed how to extend this treatment to virtual photons, that is, dileptons, with small invariant mass, again at leading order in the coupling [6]. Recently, Ghiglieri et al showed how to extend the AMY treatment of real photons to the next-to-leading order (NLO) [1]. In this section we will review the calculations in each of these papers. The NLO dilepton calculation will then consist of merging the innovations in Refs. [1,6].

Notation
Here and throughout the paper capital letters stand for four-vectors, lowercase italic letters for the modulus of the spatial three-vectors and the metric signature is (−+++), so that is the momentum of the photon, which we choose to be oriented along the z axis. We will work perturbatively in the strong coupling g, meaning that we treat the scale gT (the soft scale) as parametrically smaller than the scale T (the hard scale).
Throughout the paper we will often use light-cone coordinates, which we define as p − ≡ p 0 − p z and p + ≡ p 0 +p z 2 . This normalization is nonstandard, but we find it convenient because dp 0 dp z = dp + dp − , and because we will frequently encounter cases in which p − = 0, in which case p z = p 0 = p + with our conventions. The transverse coordinates are written as p ⊥ , with modulus p ⊥ . With this choice, we will treat the plus component of the photon momentum, k + , to be of the order of the temperature or larger, k + > ∼ T , while the other component is assumed to be of order g 2 T , k − ∼ g 2 T , so that −K 2 = 2k + k − ∼ g 2 T 2 .

Leading order Photons
In a plasma without net baryon number 2 , the photon self-energy Π < (K) is determined by diagrams with two current insertions (which we represent as external photon lines) attached to a quark loop, with an arbitrary number of gluonic corrections. AMY have shown [13] that only some of these diagrams contribute, and only in some kinematical regions. In particular, the lowest-order diagram, shown in Fig. 2, corresponds to a kinematically disallowed process. Instead the lowest order diagrams involve at least one gluon line, as depicted in Fig. 3. The figures also establish the momentum notation we will use throughout: the quark lines have momentum P and P ± K, while gluon lines have momentum Q. They also show that the rate will depend parametrically on couplings as e 2 g 2 . Figure 2. The Born diagram on the left. Its cut is associated with the square of the thermal Drell-Yan process on the right. This process is kinematically forbidden for a real photon and on-shell quarks. In evaluating the cut diagrams in Fig. 3, there are contributions where the gluon line is an on-shell external gluon state. The resulting processes, shown in Fig. 4, are leading-order for quark momentum in the p ∈ [gT, T ] momentum range, and were traditionally assumed to give the leading-order photon production rate [16,17]. We will call these 2 ↔ 2 processes or elastic processes. But another possibility is that the gluon line can be soft, q ∼ gT , and spacelike off-shell. In this case it is essential to include the Hard Thermal Loop (HTL) resummations to the gluon propagator [18,19], under which the gluon has significant =⇒ 2 + Crossings   . Two-loop diagram cut through a self-energy correction on the gluon, which corresponds to scattering-induced photon radiation (crossings not shown) spectral weight in this region. This leads to a distinct contributing kinematical region which corresponds to scattering-induced emission, as shown in Fig. 5. We will call these collinear processes or collinear splitting processes. Aurenche et al [20] first showed that these processes are also leading order and can even be numerically dominant. The reason is that the process includes a kinematical region in which the intermediate quark line in Fig. 5 is nearly on the mass shell. But this near-singularity requires the inclusion of self-energy corrections, which bring in additional diagrams by gauge invariance and the necessity to correctly represent charge conservation. Therefore, in the kinematic region where gluons are soft and spacelike (representing scattering processes), one must sum over multiple gluon exchanges, such as the diagram of Fig. 6. The interference effect this generates and the associated suppression are called the Landau-Pomeranchuk-Migdal (LPM) effect.
In [13], AMY showed that these two kinds of processes (elastic scattering when one gluon is on-shell, scattering induced emission with any number of soft spacelike gluons) are both needed in the calculation, but arise from kinematically distinct momentum regions. Therefore the computation can be separated into a contribution from each process. The easiest way to see that this is true is to consider the components of the off-shell fermion's momentum P , particularly the transverse component p ⊥ and the longitudinal component p + . As illustrated in Fig. 7, the relevant momentum regions are quite distinct when viewed with this variable. Figure 6. Example of how a cut "multi-rung" diagram corresponds to interference between processes where a photon is emitted before or after a series of scatterings (gluons connecting to X's). Such diagrams must be resummed to determine the leading-order photon production rate in the collinear kinematic range.  Figure 7. Distinct momentum regions giving rise to leading-order photon production, in terms of the off-shell fermionic momentum components p ⊥ and p + . The 2 ↔ 2 processes occur when these are comparable, the collinear processes occur when p + ≫ p ⊥ .
The required diagrams in the collinear splitting region can be resummed into an integral equation for the production rate. The analysis is quite detailed, and rather than reproduce it here, we will just quote the key results.
The contribution to Π < (K), for K 2 = 0, from elastic scattering and at lowest order in the strong coupling, arises from the scattering processes q k ′ g l ′ → γ k q l and q k ′q l ′ → γ k g l (and processes with q ↔q), which contribute (writing k 0 = |k| = k) where s, t, u are the usual Mandelstam variables and for QCD with uds quarks (2.2) is the leading-log coefficient, which we will use to normalize all contributions to Π < (K) in what follows. This expression is based on massless dispersion for the external states and vacuum matrix elements.
To evaluate this expression it is convenient to work in terms of the transfer 4-momentum P µ ≡ K ′ µ − K µ , which is the momentum carried by the intermediate off-shell quark line in diagrams giving rise to the 1/t matrix element. For small p ∼ gT the on-shell condition (K ′ ) 2 = 0 enforces p − = O(g 2 T ), and we may simplify −s/t, u/t, and the statistical functions enough to perform all integrations but the p + , p ⊥ integrals: which is log divergent. The divergence is removed because for small p ∼ gT there are HTL [18,19] medium corrections to the fermionic self-energy which change the matrix elements. We showed in [1] that, by working in the p ⊥ , p + basis, these modifications can be handled analytically and the p + integral can still be performed; the infrared limit of the p ⊥ integral is modified to where m 2 ∞ is the thermal asymptotic 3 quark mass m 2 ∞ = g 2 T 2 C R /4. Matching this behavior with the large-p behavior which is found numerically, one finds The result for C hard is a numerical fit. Next it is necessary to consider interactions with one or more spacelike off-shell gluons (that is, gluon-exchange scatterings). An efficient way to think of these processes is that the scattering provides a little off-shellness such that the process of Fig. 2 becomes kinematically allowed. This is only efficient for a narrow range of transverse momenta p 2 ⊥ ∼ g 2 T 2 , leading to a g 2 phase space suppression. The relevant diagrams are shown in Fig. 8, which emphasizes the small angular spread between the photon and the emitting particle. It is because this angle is narrow that we refer to these as collinear processes. 3 m 2 ∞ = 2m 2 q the "mass" of a quark at rest, which is used in much of the HTL literature [18].
The deviation from the mass shell is δp 0 ∼ g 2 T , which allows for large off-shell propagation distances d ∼ 1/g 2 T , which is the same order as the mean inter-scattering distance. Therefore one must resum multiple scattering, as illustrated in Fig. 6. The conceptual issues involved in such a resummation are treated in [21] and the technical issues are in [13]. Skipping the details, the current correlator from these processes, for on-shell 4-momentum K and at leading order, is given by . (2.8) Here δE = E p+k − E p − k 0 is the eikonalized energy difference between having a quark of energy corresponding to a momentum of k + p and having a quark of energy corresponding to momentum p with a photon of energy and momentum k. C(q ⊥ ) is the differential soft scattering rate, which at leading order reads [13,22] where v µ k ≡ (1, 0, 0, 1) and G rr µν is the cut (symmetric) HTL gluon propagator. Here m D is the Debye mass, (1 + N f /6) at leading order. Eq. (2.6) can be understood in the following way. Real photon production involves a current operator insertion in the amplitude, followed by time evolution, and then a current insertion in the conjugate amplitude. At times between the current insertions, the density matrix contains an off-diagonal term with a quark with momentum p + , p ⊥ and a photon of momentum K in the amplitude, but a quark of momentum (p + +k), p ⊥ in the conjugate amplitude. The size for this entry in the density matrix is f (p ⊥ ) and Eq. (2.7) is the timeintegrated evolution equation for this amplitude; δE represents the phase accumulation because the state in the amplitude and conjugate amplitude do not have the same energy, while the C(q ⊥ ) term describes the effect of scattering processes on the evolution of the density matrix. The second line in Eq. (2.6) describes the overlap of the current operator on this density matrix element, times the DGLAP splitting kernel.
Eq. (2.7) and Eq. (2.6) were solved explicitly in Ref. [12]. We will not present the details here, since we will have to solve a modified version of this relation.

Leading order Dileptons
Now let us look at the modifications necessary to consider slightly timelike K µ , that is, k + = k 0 +k 2 ∼ T as before, but k − = (k 0 − k) ∼ g 2 T rather than strictly k − = 0. This frequency difference is of order the characteristic soft plasma scattering rate, meaning that no particle in the plasma is on-shell to better precision than O(g 2 T ). Therefore such a small virtuality does not change any of the parametric analyses, and the set of diagrams is unchanged. But we must re-consider their evaluations.
First consider the 2 ↔ 2 processes. The processes shown in Fig. 4 always involve an off-shell intermediate line. Because of plasma screening, the leading-order contribution only arises when this line's momentum is > ∼ gT . Therefore the leading-order 2 ↔ 2 results are not changed by k − ∼ O(g 2 T ).
The same is not true of collinear processes, because the virtuality involved is δE ∼ g 2 T . Therefore the introduction of k − ∼ g 2 T makes an important difference, even making the tree-level (Born) process kinematically allowed if k − is large enough. Aurenche et al present a complete leading-order analysis [6] (see also [23]), showing that Eq. (2.6) through Eq. (2.8) are changed to The changes are that, first, δE receives an extra contribution from k − , and, second, there is a new contribution g(p ⊥ ) accounting for the contribution of longitudinal photons, which are absent for k − = 0. We should remark that these expressions are also valid for k − < 0, that is, spacelike K 2 , which corresponds physically to probing the plasma through deep inelastic scattering. In this regime the longitudinal term g(p ⊥ ) corresponds to a timelike photon exchange and gives a negative contribution. While the DIS regime is not experimentally accessible, we still study it, because the continuity or smoothness of Π < is relevant to the lattice reconstruction of the spectral function.

Real photons at NLO
Now we return to real photons, k − = 0, but we consider corrections at the next order in the strong coupling g. This case has been considered in detail in Ref. [1], so again we will summarize the results without presenting a detailed derivation.
In vacuum field theory, NLO corrections generally arise from a loop or an extra external leg and are suppressed by O(g 2 ). In thermal field theory, NLO corrections can be as large as O(g) whenever soft p ∼ gT momenta are involved, since bosonic lines with O(gT ) momenta can introduce large Bose factors n B (p 0 ) ∼ T /p 0 ∼ 1/g when p 0 ∼ gT . Therefore we must examine our calculation for such soft sensitivity.  Figure 9. Phase space regions in the p ⊥ , p + plane which are relevant at NLO. The collinear region expands towards small p + and larger p ⊥ , while the 2 ↔ 2 region expands to smaller p ⊥ . They overlap in the semi-collinear region.
The 2 ↔ 2 calculation receives O(g) corrections because an O(g) fraction of the contribution in Eq. (2.1) arises when the argument of n B (l) or n B (l ′ ) is O(gT ). In this region our treatment of the external state dispersion and spectral weight break down, so we must reconsider this O(g) contribution. In this region s, t ∼ gT 2 , so the other momentum combinations are at a relatively small opening angle θ ∼ √ g. Therefore we call this a semi-collinear contribution. There can also be O(g) corrections more generally when the exchanged fermion becomes soft, t ∼ g 2 T 2 . We illustrate this expansion of the relevant P -momentum region in Fig. 9: while at leading order we need the blue shaded region, at NLO we need the region enclosed in the larger blue contour.
The collinear calculation has O(g) corrections from several places. First, the ingredients in Eq. (2.7) receive O(g) corrections, where [25] and the explicit expression for δC can be found in Ref. [24]; both effects are summarized in [1] Eq. (B.27) and Eq. (C.1). In addition, an O(g) portion of the contribution arises from p + ∼ gT and an O(g) portion arises from p 2 ⊥ ∼ gT 2 . In each region at least one collinear approximation used to derive Eq. (2.6), Eq. (2.7), or Eq. (2.8) breaks down -though in the p 2 ⊥ ∼ g 2 T 2 region the collinear approximations only suffer O(g 2 ) correction. The expanded phase space region relevant at NLO is again illustrated in Fig. 9 as the larger red contour. Note that this contour overlaps with the blue (2 ↔ 2) contour in the soft region and in the semi-collinear region with p + ∼ T and p ⊥ ∼ √ g T . These regions will need special attention.
We can handle the corrections δm 2 ∞ and δC as follows. We will replace f (p ⊥ ) → f (p ⊥ ) + δf (p ⊥ ), with δf the O(g) correction arising at linear order from δm 2 ∞ and δC.
Specifically, Eq. (2.7) becomes . Expanding at zero order in δm 2 ∞ and δC we recover Eq. (2.7), while at first order we find Once f (p ⊥ ) has been determined, it acts as a source in an inhomogeneous equation for δf (p ⊥ ). We present a particularly clean methodology to treat these integral equations in Appendix A, but we do not present the results here because they will be superseded in the next section. The corrections remain subdominant in the large-p ⊥ and small-p + regions, so there is no overlap between these corrections and those from the kinematic edges of the collinear region.
Treating the kinematic edges, that is, the soft and semi-collinear regions, turns out to be surprisingly easy. The solution is to perform the leading-order calculations as if nothing has changed, to evaluate the NLO collinear corrections from δm 2 ∞ and δC, and to add the following next-to-leading correction from the overlap regions: where C soft+sc (k) is given by Eq.(6.3) of Ref. [1].

Dileptons at NLO
In the previous section we saw how to extend the leading-order photon production calculation of Arnold, Moore, and Yaffe in two ways. We extended it to nonzero but small photon virtuality; and we extended it to next-to-leading order. Now we want to manage both extensions simultaneously. To do this, we need to figure out how the NLO corrections detailed in Subsection 2.4 get modified by the virtuality discussed in Subsection 2.3.

NLO effect of virtuality
A key step in the leading-order dilepton treatment was showing that the 2 ↔ 2 processes are not disturbed at leading order by the introduction of a small photon virtuality k − ∼ g 2 T .
We now show that this persists to NLO (defined as O(g)).
To do so, we will have to look momentum region by momentum region. The hard scattering region (Eq. (2.1) with l, l ′ , k ′ ∼ T ) is the simplest to analyze. An O(g 2 T ) shift to k 0 changes the statistical functions, the matrix element, and the phase space only by O(g 2 ) amounts, and can therefore be neglected. The semi-collinear region is also straightforward. If we consider the u/t term and l ∼ gT , we find that the shift to k 0 modifies the value of l by an O(g 2 T ) amount. This will change the phase space, n B (k ′ ), and u/t by O(g) amounts. But since this region is already only an O(g) correction, the impact of these shifts is O(g 2 ). We can also see this by considering the semi-collinear region from the collinear side. If we consider Eq. (2.13) in the region p 2 ⊥ ∼ gT 2 , an O(g 2 T ) value for k − is an O(g) correction to the p 2 ⊥ term. And the contribution in Eq. (2.10) from g(p ⊥ ) is O(g) relative to the contribution from f (p ⊥ ), due to the p ⊥ factors associated with the latter. Therefore, in this kinematic region we again find that k − ∼ g 2 T leads to an O(g) correction to an already O(g) suppressed region, and can be neglected.
The analysis of the soft region is more subtle. We must reconsider the approximations leading to Eq. (2.3), allowing k − ∼ g 2 T . The shift has only an O(g 2 ) bearing on the statistical functions and l, l ′ integrations; but it will shift the value of P µ . Specifically, the on-shell condition ( 3), Eq. (2.4), we assumed that p − = 0; while in fact it is smaller than p ⊥ by a relative factor of g. But we showed in Ref. [1] that finite-p − corrections to Eq. (2.4) only arise at even orders in p − . Therefore a relative-order g suppressed correction to p − leads to an O(g 2 ) correction, which can be neglected. Hence the soft and semi-collinear corrections, summarized in Eq. (2.17), are not modified at O(g) order.
So the only place where k − ∼ g 2 T can affect Π < at O(g) is in the collinear phase space region. The collinear expansion which establishes Eq. (2.10), Eq. (2.11), Eq. (2.12), and Eq. (2.13) is an expansion in k − ≪ T , p 2 ⊥ ≪ T 2 , so a small O(g 2 T ) correction will only modify the structure of these equations at the O(g 2 ) level. So the only modifications we need at the NLO level are the same ones we needed for real photons; the NLO corrections δC and δm 2 ∞ . Hence, an NLO calculation of the dilepton rate requires the NLO treatment of real photons from 2 ↔ 2 and semi-collinear regions, and a treatment of Eq. (2.10) to Eq. (2.13) but including NLO corrections using Eq. (2.16) and the analogous expression for δg(p ⊥ ), We specify our numerical procedure in Appendix A.

Emergence of the Born term
As shown in [6], Eq. (2.10) includes the Born term, i.e., the contribution from the thermal Drell-Yan process show in Fig. 2. Let us explore how this occurs, and use it to study how the solutions behave as k − becomes larger, in particular for is treated as formally much smaller than δE, Eqs. (2.11) and (2.12) can be solved by substitution order by order in the collision operator. The zeroth order, f (0) (p ⊥ ) and g (0) (p ⊥ ), corresponds to having no exchanged soft gluons. In more detail one has for the transverse function so that where we have included an iǫ prescription to account for a very small contribution from C. The real part only arises due to the iǫ prescription, and is a delta function at the point where the denominator vanishes. One then obtains (3.4) The first term is the coefficient on p 2 ⊥ in the denominator of Eq. (3.3), the second term is the value of p 2 ⊥ in the numerator, and the Heaviside function makes sure to consider only cases where the denominator vanishes for positive p 2 ⊥ . The longitudinal mode is evaluated in the same fashion, leading to a complete Born term of [6,26] . Because the quarks have (effective) mass m 2 ∞ , the Born term is only kinematically allowed for −K 2 > 4m 2 ∞ , as reflected by these integration limits. Next consider k 0 < 0 (the DIS regime). Eq. (3.4) is valid for either sign of k − , but the Heaviside function picks out a different range of p + , namely p + > p + DIS = k + /2( 1 + 4m 2 ∞ /K 2 − 1) (and p + < −k − p + DIS , which is symmetrical). Then Eq. (3.5) is replaced with (3.6) The kinematic edge occurs at k − = 0; but for |k − k + | ≪ m 2 ∞ p + must be large, leading to exponentially small statistical factors. Also note that the second term in square brackets is negative. It actually dominates at large p + , and Π < coll,Born,DIS itself is negative for small values of |k − |. This apparently unphysical result is because the negative second term in Eq. (3.6) arises from the real part of a quark HTL self-energy and vertex; the imaginary parts of the HTLs appear in the 2 ↔ 2 processes, and the sum is always positive as it should be.
To understand the behavior at larger K 2 and in particular to make contact with the calculations [7,8] which have been carried out there, it is instructive to expand Eq. (3.5) in large −K 2 /m 2 ∞ . The first step is to estimate the importance of including the collision operator in Eq. (3.3). It provides a width to the on-shell peak, and it provides a small imaginary part off-peak. Both effects are suppressed by O([m 2 ∞ /(−K 2 )] 2 ). So at NLO in this expansion we can neglect the collision operator, and start with Eq. (3.5). The k − /m 2 ∞ term in Eq. (3.5) is larger than the (p + ) 2 +(p + +k + ) 2 2k + p + (p + +k + ) term by one factor of −K 2 /m 2 ∞ , so the leading-order contribution is obtained by integrating k − /m 2 ∞ from p + = −k + to p + = 0 (the small-m 2 ∞ limits of p + min,max ), obtaining which we have labeled "collinear free" because it is precisely the Born term that would arise in the free g → 0 limit and in the collinear limit k − /k + → 0.
At the next order we must account for the narrower limits of the p + integration in this dominant term; p + max ≃ − m 2 ∞ 2k − and p + min = −k + − p + max . This cuts off each edge of the integration range, reducing the result by 2k − is the amount of the integration region which is removed, and the remaining expression is the integrand evaluated at the edge (p + = 0 or p + = −k + ). We must also account for the subdominant term in the integrand of Eq. (3.5): where the extrema arise from having exploited the symmetry around p + = −k + /2 and in having expanded p + max as before. At small p + the statistical functions become 1 2 and contributing a log to the integral. Adding and subtracting this limiting behavior, we find The integration on the r.h.s, which can now be safely pushed to p + = 0, needs to be performed numerically. We note that the same integration was dealt with in the context of the NLO real photon rate in [1] and an accurate fit of the result as a function of k + /T was provided. Making up for slight differences in notation we have 11) where the fit for C pair can be read off from Eq. (D.17) in [1]. Our result agrees with the one in [8].
Combining the pieces, coll,Born,lim (K) + Π < coll,Born,θ (K) . Similarly, if we expand Eq. (3.6) in |k − k + | ≫ m 2 ∞ , we find (at leading order) Note that Eq. (3.5) is dominated by p + ∼ k + /2, whereas Eq. (3.6) is dominated by p + ≤ T . Therefore for k + ≫ T the collinear expansion works much better in the dilepton than in the DIS regime, and Eq. (3.13) will have a narrower range of validity as we go to larger |k − |.

Results
Let us start by collecting the leading order result, which, as we argued, is given by the 2 ↔ 2 rate, unaffected by a nonzero k − , and the collinear rate. In detail 2↔2,LO (K) was introduced and C hard was given in Eq. (2.5) and Π < coll (K) (and correspondingly C coll ) is given by solving Eq. (2.10). We have made clear the dependence of collinear processes on k − and on m 2 ∞ /m 2 D = 3C F /(4(N c + N f /2)). At NLO we must combine the soft and semi-collinear contributions, given by Eq. (2.17), with the collinear ones, obtained by perturbing Π < coll (K) with δC and δm 2 ∞ , analogously to Eq. (2.16). In detail with the O(g) correction δΠ < (K) reading with C soft+sc given in Eq.(6.3) of Ref. [1]. These first two terms come from the soft+semicollinear contribution, Eq. (2.17), and are the same as for the real photon case. δC coll arises from solving Eq. (2.16) modified by using Eq. (2.13) for δE, plus an analogous contribution from δg. Finally, we remark that m 2 ∞ /m 2 D in Eq. (4.3) is to be intended at LO, as given in the previous paragraph.

Formal weak-coupling limit
First consider the formal weak coupling limit, where all results have been derived. In this case k − ∼ g 2 T can be neglected compared to k + ∼ T in evaluating all statistical functions, and we need not distinguish between k, k + , and k 0 . Then we can evaluate C coll and δC coll using Eq. (2.10) to Eq. (2.13) and with Eq. (2.16) (and an analogous expression for g).
We present results for the specific case N f = 3 in Fig. 10 and Fig. 11, which show C coll and (1/g)δC coll respectively at a few values of k + . We also tabulate the results in Appendix B, and provide a code for evaluating C coll and δC coll in the extra included files. Fig. 10 shows that the result rapidly approaches the Born value as k − /g 2 T becomes large. Figure 10. Contribution of collinear processes to dilepton production at leading order, expressed as C coll . That is, Π < is given by B times the indicated value. In the left plot we display the full value; on the right side we subtract the (free) Born contribution, Eq. (3.7) and Eq. (3.13). The cusp in the right plot is because the Born term has a cusp (see Fig. 1) which is subtracted from a smooth function. Figure 11. Subleading correction to the collinear contribution, expressed as (1/g)δC coll . In an NLO computation, g times this result should be added to C coll . Therefore we have also presented the difference between the LO result and the Born value. The subtraction introduces a cusp at k − = 0, since the Born value possesses such a cusp. As Fig. 11 shows, the NLO correction always increases the rate, but it is a complicated function of k − /g 2 T .

Treating the coupling as finite
Phenomenologically, we are interested in obtaining results for specific finite values of gespecially we want to push towards values of physical interest in heavy ion collisions, even though the coupling expansion will break down in this limit. When treating k − as formally infinitesimal, we did not need to distinguish between k, k 0 , and k + at NLO. However we must make some sensible choice in which one to use when we treat k − as a finite number times T . Our choice is the following. Consider a plot of Π < for a fixed k as a function of k 0 . We treat B = B(k 0 ) (so the main statistical function controlling the overall size of our results is in terms of k 0 as it should be). Then we treat k + as equal to k where they appear inside the integral in Eq. (2.10) and in Eq. (2.11), Eq. (2.12), and Eq. (2.13); and k − is determined as expected, k − = k 0 − k.
In addition, once we start considering k − ∼ T , the collinear approximations we have made stop being as accurate. In this regime it makes more sense to perform an expansion which treats k − and k + on equal footing. Of course, this has already been accomplished by Laine [7]. So it would be valuable to improve our result in this way. As we have shown in Eq. (3.7), the large-k − limit of the collinear part is enhanced by k + k − /m 2 ∞ over the remainder of the LO result, giving rise to the collinear free (Born) term. We can subtract this term and add the free Born term, computed taking finite angles into full account. The finite-angle Born term is well known and reads [4] which indeed reduces to Eq. (3.7) at first order in k − /k + and is exactly the leading-order result in the K 2 ∼ T 2 calculation of [7]. Similarly, in the DIS region it reads which again reduces to Eq. (3.13) at first order in k − /k + . Therefore, our first step in making our results reliable at larger k − is to subtract Eq. (3.7) and replace it with Eq. (4.4). That leads to 4 and leaves the NLO correction δΠ < unmodified, so that Such a prescription will thus ensure that the leading, large-K 2 behavior of our LO and NLO calculations agrees with the leading order in [7]. The next-to-leading order corrections to that calculation arise from the diagrams in Fig. 3, evaluated without resummations for all possible cuts, thus including both virtual corrections to the Born term and real corrections in the form of 2 ↔ 2 processes, properly taking into account the no-longer negligible virtuality. A procedure to patch that calculation with an LPM-resummed one for small k − was recently introduced by Ghisoiu and Laine in [8]. From the point of view of our calculation, it requires, on top of the previous substitution of the free term, the subtraction of the remainder of the large-K 2 behavior of the leading order, i.e., Π < coll,Born,m 2 ∞ (K) and Π < 2↔2,LO , and its replacement with the NLO correction for large K 2 obtained by Laine in [7]. In detail, where Π < N LO,Laine is the O(g 2 ) correction to Π < free (K) computed by Laine. As shown there, it diverges as ln(T 2 /−K 2 ) at small K 2 . Similarly, Π < coll,Born,m 2 ∞ (K) contains a ln(m 2 ∞ /−K 2 ), which is correct at large k − but divergent and unphysical at small k − . The subtraction leaves a ln(T 2 /m 2 ∞ ), which is exactly the logarithmic term we encounter at LO.  Figures 12 and 13 show a comparison of the two leading-order prescriptions (4.6) and (4.8) for two different values of k/T for the case of 3 light flavors and at a few gauge couplings. For practical reasons [to eliminate a trivial overall e −k 0 /T behavior], we plot the spectral function ρ J rather than Π < (we recall that, at equilibrium, Π < (K) = n B (k 0 )ρ J (K)). In both figures, the continuous lines represent the results obtained from Π < LO,large at three different values for the coupling, whereas the dashed lines correspond to Π < LO,GL , which is unfortunately only available for k 0 > k.  Figure 13. The same as Fig. 12, but for the case k = 3T .
Interestingly, the two prescriptions agree at k 0 = k but seem to deviate significantly already at very small k − . We believe that at small k − the discrepancy can be attributed in part to our arbitrary choice regarding the k 0 , k + , k ambiguity, in part to the subleading dependence of the 2 ↔ 2 processes on k − , which we omit but Laine includes 5 , and in part to a potential unreliability of the Ghisoiu-Laine procedure at −K 2 ≪ 4m 2 ∞ (see footnote 5 in [8]). We furthermore remark that in App. C we show by an explicit calculation that, for k − ∼ gT , Π < LO,large is correct at O(e 2 g 2 T 2 ), that is, our calculation is fully NLO in this regime as well as for k − ∼ g 2 T .
On the other hand, a significant deviation at larger values of k − /T is to be expected and there the reliable curve is certainly the dashed one, which for k + , k − ≫ T is dominated by the (finite-angle) Born term and its finite vacuum corrections, which are completely absent from our calculation. Conversely, our calculation asymptotes in this region to the Born term supplemented by a g 2 T 2 and a g 2 T 2 ln(kk − /T 2 ) correction. Note that, at least for very large k − ≫ T , we know based on Operator Product Expansion techniques [27] that g 2 T 2 terms, with or without accompanying logarithms, will be absent.
In order to better highlight the differences between the two prescriptions, in Fig. 14 we subtract off the Born term Π < free from each curve in Fig. 12. One can then more easily see the deviations at small k − at intermediate and large coupling, as well as the emergence of the vacuum correction to the Born term in the dashed lines that give rise to the approximate linear behavior at larger k 0 .
Of course, the most accurate estimate of the small virtuality region we can make is the one where we also include the NLO corrections we have found: 5 At finite k − one also has 3 ↔ 1 processes, which are also included in Laine's calculation. Note that none of the effects which give rise to δΠ < (K) involve diagrams which are included in Laine's calculation [7], and no NLO contribution grows as a power 6 at large k − /g 2 T , so no new subtraction is called for here. Figures 15 and 16 are then the NLO counterparts of Figs. 12 and 14, with δΠ < added to both full and dashed lines. These represent our best estimate of the spectral function relevant for dilepton production in the small virtuality region.  a few percent across the light cone, consistently with what is observed in the real photon case. At larger couplings the dip disappears, to be replaced by a maximum at positive k − and a minimum on the opposite side. The large corrections observed at negative k − are not to be trusted; as discussed below Eq. (3.13), the collinear approximations in our treatment break down for k − < −T .

Conclusions
We have computed to next-to-leading order the production rate for dileptons with a small virtuality, parametrically of order g 2 T 2 . This required a careful merging of the NLO result for real photons and of the LO result for small-mass dileptons. To this end, we have reviewed these calculations in Sec. 2. We showed how the LO photon rate arises from two distinct processes, elastic scattering and collinear, inelastic splitting, corresponding to separate kinematical regions. The LO dilepton rate at small K 2 requires only a modification of the latter processes to account for the non-vanishing virtuality. When computing the photon rate to NLO, the two processes and the corresponding kinematical regions start to blur and merge, requiring careful analyses which we reviewed. In particular, an intermediate process/region arises, the semi-collinear one. In Sec. 3 we have shown that, similarly to LO, the extension of the NLO photon rate to finite virtuality requires only a modification of the purely collinear part of the NLO photon rate. The effect of the finite virtuality on the other processes is negligible to this order. The collinear rate to leading-and next-to-leading order is obtained by solving a set of differential, Schrödinger-like equations that resum the effects of the multiple collisions of the quarks that emit/annihilate into the virtual photon with the lightlike medium constituents. We have shown in detail how the equations governing the photon rate are modified and we have introduced in App. A a new technique for solving them. We also provide code and tables of the results in App. B. We have furthermore shown that the zeroth order in the resummation of collisions, i.e., no collisions, is (the collinear limit of) the Born term, that is the direct annihilation of a quark and an antiquark into the virtual photon. Due to the presence of thermal masses for the quarks, such a term only exists for −K 2 > 4m 2 ∞ , and it becomes dominant when −K 2 ≫ 4m 2 ∞ . Numerical results from the solution of these equations, together with the other terms that are taken, unmodified, from the photon calculations, are presented in Sec. 4. We remark that, while the dilepton rate is defined only for positive k − , the related electromagnetic Wightman function Π < is defined for any k − and our calculation is valid for k − ∼ g 2 T , irrespective of its sign. We are then able to analyze the behavior of this function across the light cone, which is of importance also for the ongoing efforts at reconstructing this function from lattice measurements of Euclidean correlators.
We find that at small couplings Π < resembles the free one shown in Fig. 1, with the cusp smoothed by the small (O(g 2 T 2 )) finite photon rate. At larger couplings the transition from the spacelike to the timelike region is smoother, resembling more the results obtained at strong coupling in SYM. These findings are highlighted by Figs. 12 to 16. These figures also show the two different prescriptions we have adopted to stretch the validity of our calculation to larger values of k − . The first entails the replacement of the collinear limit of the Born term with its full counterpart, whereas the second, originally due to Ghisoiu and Laine, patches Laine's calculation at large k − with the LPM-resummed one at small k − . We observe deviations between the two prescriptions. At large k − the latter one is to be trusted, whereas at smaller values we believe the first one to be more accurate. The NLO corrections we have found generally increase the dilepton rate, by up to 30-40% for α s = 0.3 and small virtuality, but they make little difference at larger virtuality. (This is in contrast to the large-virtuality NLO correction found by Laine, which is more important at large virtuality but which is subsumed into our LO calculation at small virtuality.) The tabulated form and the code we provide can be used to readily incorporate our results into phenomenological analyses taking into account the hydrodynamical evolution of the near-equilibrium plasma produced in heavy-ion collisions.
Finally, we should address two of the questions which motivated this study. First, experimentalists have sometimes interpreted the dilepton rate at small mass but large energy as a valid substitute for the real photon production rate. For instance, in Ref. [3], the PHENIX collaboration considers di-electrons with K 2 < (300MeV) 2 and k ∈ [1, 5]GeV. Assuming for the sake of argument a temperature of 250 MeV, this corresponds to k ∈ [4,20]T , and k − = K 2 /(2k + ) < 0.12 T for k + = 6T . So for an estimate of whether Π < is the same for a real photon and a dilepton pair in this virtuality range, we should compare the value of the top solid (blue) curve in Fig. 15 at k 0 /T = 6 and k 0 /T = 6.12. The solid blue curve shows a change of 4%; the dashed curve (which is probably the less reliable one in this range!) changes 25%. Even using the unrealistically small α s = 0.1 (green curves), the change is 20% or 40%. Therefore we find that the difference in electromagnetic response between photons and such small mass dileptons is negligible, within our perturbative computational framework. This supports the interpretation of such low-mass dileptons as a legitimate replacement for a real photon measurement.
The other question we wanted to address is whether to expect sharp features in the spectral function, which could upset an analytical continuation of lattice data. Fig. 15 suggests that the spectral function "turns on" over a range of 1-2 T width, and shows no true "notch" feature. This is much smoother than we might have worried by looking at the free case, Fig. 1. Whether it is smooth enough is a question we will leave to the experts in the field. But we close with the remark that it may be easier to determine this with a better understanding of the small virtuality behavior of Laine's un-resummed NLO calculation [7] and an extension of this calculation to k − < 0 (the DIS regime). the observation of Aurenche et al [6] that the equations are most easily solved by Fourier transforming from p ⊥ to b (impact parameter space). The Fourier-transformed versions of Eq. (2.11) and (2.12) read . We can make Eq. (A.1) a scalar equation and implement the LHS of Eq. (A.1), Eq. (A.2) as boundary conditions at b → 0: The quantities one has to solve for, appearing on the second line of Eq. (2.10), are A further boundary condition is given by the requirement that f and g decay to zero as b → ∞. Our strategy is to start at large b and use a standard ODE solver to evolve Eq. (A.3) inwards, multiplicatively adjusting the solution at b → 0 so as to obey Eq. (A.5). Let us define m 2 eff as As shown in [6], generic solutions of f and g at large values of b are given by a superposition of exp(± m 2 eff b), with slight perturbations caused by C(b → ∞) ∼ log(bm D ). As long as m 2 eff > 0, the exponents are real and it suffices to initialize the numerical solver at an appropriately large value of b to pick the shrinking solution. On the other hand, when the mass is negative, i.e., where a Born term is present (see Eq. (3.5)), the behavior at large b is mostly oscillatory, with a small damping introduced by the presence of C(b). This can cause stiffness in the numerical solution, making numerical precision challenging. This is especially so for −K 2 ≫ 4m 2 ∞ . It can also be numerically delicate to treat the f (b) equations near b = 0, where Eq. (A.5) shows the solution blows up, but Eq. (A.7) shows we need to extract a finite piece. We handle this problem by solving the ODE via a technique which takes into account the dominant oscillating behavior analytically and solves only for the interesting corrections due to C(b), δm 2 ∞ . We start by rescaling f (b), g(b) to bring Eq. (A.3), Eq. (A.4) into a standard form, Specifically, and similarly forg; d = 3 for f and d = 1 for g.
The independent solutions to the K(b) = 0 equation are known analytically. Call them 7 J(b) and Y (b), obeying and likewise for Y , with Wronskian with c a constant which depends on the normalization convention for J, Y . For A > 0 J, Y are expressed in terms of Bessel and Neumann functions, while for A < 0 they involve growing and shrinking modified Bessel functions I, K; but we will write J, Y in the following for notational simplicity. In either case we choose Y to be the divergent and J to be the finite solution at the origin. The most general solution to the K(b) = 0 equation is With this in mind, we parametrize the actual solution as This form is underdetermined, as we have doubled the number of freedoms. We eliminate the extra freedom by imposing the condition which together with Eq. (A.12) turns Eq. (A.10) into the relation Multiplying by either Y (b) or J(b) and using Eq. (A.15) and Eq. (A.13) we obtain coupled first-order equations for c 1,2 : The For the NLO corrections one must perturb the original equations: . Thenf (b) still obeys Eq. (A.10) and acts as a source term for δf : Writing δf = δc 1 J(b) + δC 2 Y (b) and applying an analogous procedure, we find (A.20) When solving this numerically in the transverse case, one must take into account that δK(b → 0) ∝ b, which makes δc ′ 1 nonzero at the origin. Note that both δc 1 (0) and δc 2 (0) give rise to a correction: at linearized order Im (A.21)

B Tabulated results
Here we present our strict leading-order and next-to-leading order results, for N c = 3 and N f = 3 and several values of k. Specifically, we present results for C coll and for (1/g)δC coll for k = 3T, 4T, 5T, 6T, 8T, 10T and a range of k − values, all parametrized in terms of k − /g 2 T . These results are independent of the value of the coupling and can be readily plugged in Eqs. (4.1) and (4.3).
We have also provided a C code which evaluates these quantities for input values of k and k − (in units of T and g 2 T ), and which allows the user to specify N f and N c . Users can use this code to develop a table of results for interpolation.

C Extension to k − ∼ gT
Here we will show that, if we apply our calculation outside of its range of validity by considering k + > ∼ T but k − ∼ gT , then our leading-order calculation actually produces the  Table 2. Next-to-leading correction to the collinear contributions to the dilepton rate, for N f = 3 QCD (N c = 3), for 6 values of momentum k and several virtualities.
correct leading and next-to-leading contributions if we subtract off the free collinear Born term and add back in the free, finite-angle Born term.
The leading-order contribution is given by evaluating the Born diagram, Fig. 2. Mo-mentum conservation and the requirement that both propagators be on shell forces P to be semi-collinear, i.e. P 2 ∼ gT 2 . Hence, thermal masses can be dropped at LO and the result is exactly the one given in Eq. (3.7), which is of order αK 2 ∼ αgT 2 . We however note that, in obtaining Eq. (3.7), we have stretched p + to 0 and −k + . Since we encounter no singular behavior, we are allowed to do so at leading order, but it will require a subtraction in the soft region at the next order to avoid a double-counting. The next order, O(αg 2 T 2 ), comes from the hard 2 ↔ 2, semi-collinear and soft regions. As we showed at the end of Sec. 3.2, the purely collinear region p ⊥ ∼ gT , p + ∼ T no longer contributes, as these processes are now off shell, and collision corrections are suppressed by O((m 2 ∞ / − K 2 ) 2 ). In the hard region, and thus for p 2 ⊥ ≫ gT 2 , only 2 ↔ 2 processes are kinematically allowed. They are still not affected by the virtuality in a first approximation, since P 2 ∼ T 2 ≫ K 2 . However, one needs in principle to treat the approach to softer momenta with greater care, as there are changes both to the kinematics and to the matrix elements. Eq. (2.3) turns into where on the first line the θ-function comes from accounting for k − in the kinematics in Eq. (2.1) and the extra term in round brackets likewise arises from a modification to the Mandelstam variable s for p + , p ⊥ , k − ∼ gT . Strikingly, the two corrections cancel, so that the approach to the soft sector is the same as in Eq. (2.3). This cancellation relies on our choice to leave the p ⊥ integration for last, which proves very convenient in this case. Cutting off the hard sector at a transverse regulator µ ⊥ gives Π < 2↔2,LO (K) = B ln For what concerns the O(g) corrections to the semi-collinear region, they come from two sources. The first is by considering the next order in k − /k + . The replacement of the collinear free term with the exact one, as introduced in Eq. (4.6), automatically takes into account all orders in that expansion. We can employ the same procedure here. The second source comes about by considering the virtual cuts of the diagrams in Fig. 3 for P 2 ∼ gT 2 . These diagrams can be computed in a rather straightforward way, since for a semicollinear P the effect of the Q integration is to give rise to the on-shell limit of the hard self-energy and vertex correction, which are equal to the same limit of their HTL counterparts, giving rise to terms proportional to the asymptotic mass. In detail Π < N LO,semi−coll (K) = B ∞ −∞ dp + n F (k + +p + )[1 − n F (p + )] n F (k + ) where the δ ′ terms come from the self-energy diagrams and the δ term comes from the vertex correction diagram. We are regulating the transverse integration with gT 2 ≫ µ ⊥ ≫ g 2 T 2 to avoid overlap with the soft region. Using where we have neglected the boundary term at p ⊥ = ∞, as it cannot contribute to the p + integration. We then obtain (C.5) A comparison with Eq. (3.10) shows that this amounts to having replaced m 2 ∞ with µ 2 ⊥ in the extrema of the integration. Hence, one can easily modify the results obtained for Π < coll,Born,lim and Π < coll,Born,θ in Eq. (3.12), obtaining