Interpolation of hard and soft dilepton rates

Strict next-to-leading order (NLO) results for the dilepton production rate from a QCD plasma at temperatures above a few hundred MeV suffer from a breakdown of the loop expansion in the regime of soft invariant masses M2 ≪ (πT)2. In this regime an LPM resummation is needed for obtaining the correct leading-order result. We show how to construct an interpolation between the hard NLO and the leading-order LPM expression, which is theoretically consistent in both regimes and free from double counting. The final numerical results are presented in a tabulated form, suitable for insertion into hydrodynamical codes.


Introduction
Consider µ − µ + or e − e + pairs produced thermally from a quark-gluon plasma at a temperature T > ∼ 150 MeV, with the pair having a non-zero total momentum k ≡ |k| ∼ GeV with respect to the plasma rest frame, and an invariant mass M ∼ GeV. If no zero-temperature resonance lies near the M considered, non-thermal backgrounds for the production of such dileptons are expected to be smaller than for on-shell photons, and dileptons may therefore constitute a good probe of QCD interactions at finite temperature. As a particular reflection of deconfinement and chiral symmetry restoration, a relatively smooth shape is anticipated for the thermal dilepton production rate, with a characteristic overall magnitude (to be determined by theoretical computations) and an exponentially damped spectral shape.
Many different approximation schemes and kinematic regimes have been considered for thermal dilepton production in the past. First next-to-leading order (NLO) analyses were carried out long ago for k = 0 [1][2][3], finding that for M ∼ πT radiative corrections are in general small. However, moving to a "soft" invariant mass M ∼ gT with still k = 0 (here g 2 ≡ 4πα s ), a major enhancement of the rate was found after carrying out Hard Thermal Loop (HTL) resummation [4]. Most of the past work has concentrated on M ∼ gT JHEP10(2014)083 but large spatial momentum (k > ∼ πT ). In this regime the NLO rate has a logarithmic singularity, which is regulated by Landau damping of the quarks mediating t-channel exchange [5,6]. In addition, there are finite terms which all contribute at the same order because of collinear enhancement, and need to be handled through Landau-Pomeranchuk-Migdal (LPM) resummation [7] (LPM resummation incorporates HTL resummation in an approximation valid for k ≫ gT ). In order to avoid double counting, LPM resummation needs to be carefully combined with other processes [8]. A resummation beyond HTL (based on effective kinetic theory) is also needed at k = 0 for M ≪ gT [9]. In contrast, for M ≫ πT no resummation is needed at NLO, and the analysis can be greatly simplified by making use of Operator Product Expansion (OPE) techniques, with the results available in analytic form [10]. Unfortunately the OPE expansion shows convergence only quite deep in the hard regime (M ≫ 8T ) [11]. Finally, lattice simulations are being carried out at k = 0 [12][13][14] and at k = 0 [15], even though the usual issues with analytic continuation imply that the results may suffer from uncontrolled systematic uncertainties [16].
As is clear from the previous paragraph, many computations have concentrated on special regimes in which one or the other kinematic simplification can be made. The current study is the continuation of a recent project [17,18] which led to the determination of the NLO dilepton rate for generic momenta and invariant masses k, M ∼ πT [11]. The goal of the present study is to present a smooth interpolation between these hard NLO expressions, and leading-order LPM resummation in the soft regime M ∼ gT , k ≫ M . The interpolated results turn out to have a qualitatively correct behaviour even when extrapolated down to M ∼ gT , k ∼ 0. Therefore, for practical purposes, we hope that they yield a fair estimate of the thermal dilepton rate from a deconfined QCD plasma for the invariant masses and spatial momenta of interest to the current heavy ion collision program.
The plan of this paper is the following. After defining the observables to be considered in section 2, we briefly review the status of hard NLO computations in section 3 and of soft LPM resummation in section 4 (we also introduce an efficient method for the numerical solution of the LPM equations). A way to consistently combine these approaches is explained in section 5. Numerical results, meant for phenomenological use, are displayed in section 6, and we finish with a brief conclusion and outlook in section 7.

Basic definitions
To leading order in α e ≡ e 2 /(4π) [19][20][21] and omitting power-suppressed corrections from Z-boson exchange, the production rate of µ − µ + pairs from a hot QCD medium, with a total four-momentum K ≡ K µ − + K µ + ≡ (k 0 , k), can be expressed as Here n B is the Bose distribution, and ρ NS and ρ SI denote spectral functions in the "nonsinglet" and "singlet" channels, respectively, with the quark flavours assumed degenerate JHEP10(2014)083 for simplicity. For N f = 3 the singlet channel drops out, and in the following we concentrate on where c denotes a connected quark contraction; η µν ≡ diag(+−−−); and X is an integral over the spacetime volume. According to eq. (2.1), ρ NS must be negative, so we mostly discuss in the following, where Π R refers to the retarded correlator. Let us inspect separately the "transverse" and "longitudinal" parts of − Im Π R . Choosing for convenience the z-axis to point in the direction of k, The remaining longitudinal part can be expressed as where we made use of a Ward identity relating Im Π R,00 and Im Π R, 33 . The physically relevant combination is

NLO dilepton rate for general momenta
The observable of eq. (2.7) (although not separately its two parts Im Π R,T , Im Π R,L ) is currently known up to NLO in a strict loop expansion [11]. However only the leading-order (LO) result can be given in analytic form: Here light-cone momenta and a photon invariant mass were defined as For future reference, it is helpful to present eq. (3.1) also in a form before a final integration. We do this separately for the parts in eqs. (2.5) and (2.6): It is seen that a substantial cancellation takes place when adding up eqs. (3.3), (3.4). At NLO, it is more cumbersome to work out analytic expressions. However a convergent 2-dimensional integral representation can be given [11]. An analytic result is obtained on one hand for the dominant logarithmic divergence at M ≪ πT [22], and on the other for M ≫ πT [10]. Let us define an "asymptotic" thermal quark mass by where n F is the Fermi distribution, p ≡ |p|, and p ≡ d 3 p/(2π) 3 . Then the soft divergence reads whereas the asymptotic expansion in the hard limit is given by

Basic equations
Leading-order Landau-Pomeranchuk-Migdal (LPM) resummation for the dilepton production rate was worked out in ref. [7]. The dilepton case is a generalization of the on-shell photon production rate that had been considered previously [23,24]. A field-theoretic derivation of the basic equations can be found in ref. [25], and yet another approach yielding the same dynamics, operating within the imaginary-time formalism, in ref. [26]. In its usual formulation LPM resummation assumes the kinematics k 0 ≫ gT and k 0 − k ≪ k 0 . Then only the leading terms in a Taylor expansion around the light cone k 0 = k are relevant. Parametrizing the kinematics through k 0 and M 2 , this means that the spatial momentum k can be expressed as and the validity of this expansion is assumed in all formal manipulations of the present section. (Numerically, however, we may at times exit the regime in which eq. (4.1) is literally accurate; the procedure to be followed in these cases is discussed below.) Because of the assumption k 0 , k ≫ gT , Hard Thermal Loop (HTL) self-energies and vertices can be simplified through a "hard-particle" approximation (cf. e.g. ref. [27]), resulting in an effective kinetic description [28] with particles carrying "asymptotic" thermal JHEP10(2014)083 masses [29]. With a minor change of conventions with respect to ref. [7] (reshuffling of imaginary units; inversion of the sign of one of the frequency variables appearing; rescaling of wave functions; and use of n F (−ω) = 1 − n F (ω)), we are then led to define a "2-particle Hamiltonian", where ∇ ⊥ operates in the two "transverse" directions. 1 The light-cone potential is [10,30] where y ≡ |y| denotes a 2-dimensional transverse separation; is an electric mass parameter in EQCD; and K 0 is a modified Bessel function. The superscript in V + is a reminder of the fact that the potential is positive for all y > 0 (it vanishes for y = 0).
We need to solve Schrödinger equations in the S and P -wave channels: Then the functions we are interested in are The Dirac-δ constraints here correspond to energy conservation, whereas the Schrödinger equations in eqs. (4.4), (4.5) can be viewed as reflecting momentum conservation.
Of the variables characterizing external kinematics (k 0 , k, M ), only two are independent (k 2 0 − k 2 = M 2 ). As mentioned above, the derivation of the LPM equations can be justified as a leading term in a Taylor expansion in M 2 /k 2 for k ≫ gT . This implies that k 0 , k, M should be related through eq. (4.1). Sometimes, it may however be convenient to also apply the LPM equations beyond their parametric validity range. There is no unique way of doing this, however one possible criterion is that there be a specific cancellation between the transverse and longitudinal contributions, namely that the terms ω(k 0 −ω) in eqs. in a specific way. This means that we have to give up either the strict M 2 /k 2 multiplying Im Π R,00 in eq. (2.6) or the strict k − k 0 that is represented by −M 2 /(2k 0 ) in eq. (4.2). We have adopted a procedure, corresponding to ref. [7], where a compromise has been made at both points; however we have verified that the numerical effect of other choices, if made consistently, is small. Following ref. [7], it is helpful for the following to define a parameter M 2 eff originating from a combination identifiable in eq. (4.2):

Method for numerical solution
In order to solve eqs. (4.4), (4.5) numerically, we adapt to two dimensions a method employed in appendix A of ref. [31] for solving vector and scalar channel quarkonium spectral functions in three dimensions. The basic approach was introduced in ref. [32] for the vector channel (S-wave) case at zero temperature. Its idea is to reduce the solution of an inhomogeneous equation to determining that solution of the homogeneous equation which is regular at origin. By rescaling the transverse variable as ρ ≡ ym E ; introducing a coordinate ρ ′ as a handle on the behaviour of the solution under rotations; rescaling the wave functions into a dimensionless form; and making use of the parameter M 2 eff introduced in eq. (4.9), the inhomogeneous Schrödinger equations in eqs. (4.4), (4.5) can be re-expressed as specific limits of In polar coordinates, ρ = (ρ, φ), the solution of the corresponding homogeneous equation, can be written as (4.14) Among the two solutions for each ℓ, the one regular at origin (denoted by u r ℓ ≡ u < ℓ ) is of the form u r ℓ (ρ) = ρ 1/2+|ℓ| 1 + O(ρ 2 ) + i ζ ρ 9/2+|ℓ| ln ρ/ρ 0 + . . . , where ζ and ρ 0 are constants. The coefficient of the small-ρ asymptotics of the real part has been fixed in a particular way. With this normalization, and choosing the solution regular at infinity as u > ℓ (ρ) ≡ u r ℓ (ρ) ∞ ρ dρ ′ /[u r ℓ (ρ ′ )] 2 , the solution of eq. (4.10) can be written as Subsequently we obtain results analogous to eqs. (4.25) and (A.33) of ref. [31]: (4.18) Making use of the symmetry ω 1 ↔ ω 2 and carrying out one of the integrations, the final expression reads where the radial wave functions are to be solved from with the asymptotics at ρ ≪ 1 chosen according to eq. (4.15). 2 As a crosscheck we show in appendix A that eqs. (4.19), (4.20) reduce to the correct free results in the appropriate limit.

Re-expansion of the LPM result and matching with NLO
In order to combine LPM resummation with the NLO result, we need to identify those terms in the NLO result which are also part of the LPM resummation. Care must be taken in order not to count such terms twice. The identification can best be carried out by re-expanding the LPM result as a "naive" power series in g 2 , with the kinematic variables k, M treated formally as of O(πT ), because this is also the structure inherent to the hard NLO result.
The gauge coupling appears at two points in section 4.1: in the parameter m 2 ∞ of eq. (4.2), as well as in the potential V + of eq. (4.3). If we expand to zeroth order in g 2 , 2 The determination of u r ℓ and the integration over ρ in eq. (4.19) can be implemented as a simultaneous solution of nine real first-order differential equations: for Re u r ℓ , Im u r ℓ , Re ∂ρu r ℓ , Im ∂ρu r ℓ , ℓ = 0, 1, and the ρ-integral in eq. (4.19). A relative accuracy ∼ 10 −6 can be reached with modest expense for all k, M considered.

JHEP10(2014)083
eqs. (4.4), (4.5) can be solved in a Fourier representation. It is straightforward to check that eqs. (4.6), (4.7) then yield A cancellation of ω(k 0 − ω) as discussed in the paragraph following eq. (4.7) is readily verified. Summing together and carrying out the remaining integral, we get a limit of eq. (3.1): We also need the term of O(g 2 ) from the re-expansion of the LPM result. Note first that the contribution from the potential V + through eq. (4.2) is of O(g 4 ) in this counting. One way to see this is that before carrying out the final integral, the form of the potential is We see that if the propagator is expanded to O(m 2 E ), the term of O(g 2 ) drops out. In contrast, there are two contributions of O(g 2 ) from the mass term m 2 ∞ . Solving eqs. (4.4) and (4.5) in a Fourier representation and taking the cut needed in eqs. (4.6) and (4.7), m 2 ∞ changes the integration range for the Fourier momentum. In addition, it appears explicitly in the integrand, if the Fourier momentum originating from −∇ ⊥ · f is substituted by other variables as dictated by the Dirac-δ constraint from the cut. The latter contribution leads to a logarithmic divergence. Determining the logarithmic term explicitly, and making use of symmetries in order to simplify the finite terms, we get The logarithmic divergence on the first row of eq. (5.6) exactly matches that in eq. (3.7). For future reference, summing together eqs. (5.4) and (5.6), we define JHEP10(2014)083

Numerical evaluation
The goal now is to combine the NLO result from section 3 with the LPM result from section 4. We first need to subtract from the NLO result those terms that are resummed into the LPM expression, cf. eq. (5.7). After the subtraction, the "full" LPM result of eq. (4.19) can be added. Thereby the final result reads In a regime where the resummation has no effect, i.e. ∆ Im Π R | LPM = 0, we recover simply the consistent NLO result. On the other hand, in the soft regime where LPM resummation is important, the difference Im Π R | NLO − Im Π R | expanded LPM represents a hard "matching" contribution (involving 2 ↔ 2 scatterings and void of the logarithmic divergence visible in eq. (3.7)) which needs to be added to the soft LPM result.
Numerical results for ∆ Im Π R | LPM are shown in figure 1. It is seen that LPM resummation has no effect in the hard regime M ≫ gT . It does have a substantial influence in the regime M < ∼ gT , M ≪ k. The behaviour changes qualitatively at small k when the inequality M ≪ k is no longer satisfied. 3 However, for any fixed k > 0, the curves do reach a regime with M ≪ k if extrapolated far to the left. Therefore it is perhaps not completely surprising that they turn out to be qualitatively correct even for k < ∼ M (cf. figure 2(right)).
The results obtained after adding ∆ Im Π R | LPM to the NLO expression (figure 2(left)) are shown in figure 2(right). LPM resummation is seen to remove the logarithmic divergence of the NLO result at small M ≪ πT and leave over a smooth behaviour. For very small k the results show an increase which is in surprisingly good agreement with an effective kinetic theory computation relevant for this regime [9].

Tabulated spectra
Given values for − Im Π R , physical dilepton rates are given by eqs. (2.1), (2.3). In the following we refer to spectra for the production of µ − µ + pairs, but the corresponding results for e − e + can be obtained by a trivial change of the prefactor in eq. (2.1). Going over to physical units, viz.

Conclusions and outlook
The purpose of this paper has been to collect together ingredients from two existing computations, namely a consistent LPM-resummed computation of the thermal dilepton rate in a regime of "soft" invariant masses M < ∼ gT (ref. [7], with minor modifications [8]), as well as a full NLO computation in a regime of "hard" invariant masses M ≫ gT [11]. We have shown that the two different regimes can be "interpolated" into a result which is theoretically consistent in both regimes and should represent a fair approximation (with uncertainties of less than ∼50%) for all spatial momenta and positive invariant masses. The uncertainty estimate is based on a recent analysis [26] in which the equations of LPM resummation, analytically continued to imaginary time, permitted for a determination of vector channel screening masses and correlation functions which could be directly compared with lattice Monte Carlo data at T ≈ 250 MeV. Our results have been tabulated (cf. footnote 4) in a form which hopefully allows for their easy insertion into hydrodynamical codes such as ref. [35]. are not fully consistent even at LO, but they nevertheless display a qualitatively correct behaviour, as a numerical comparison with an effective kinetic theory analysis [9] shows (cf. figure 2(right)).
One way to improve upon our results would be to include NLO corrections of O( √ α s ) in the soft regime M < ∼ gT, k ≫ M , similarly to what has been done for the photon production rate in ref. [36]. A systematic study of the very soft regime M < ∼ gT , k < ∼ M could also be envisaged, by making use of effective kinetic theory and incorporating full HTL structures where necessary. Furthermore the results could in principle be extended into the spacelike domain M 2 < 0, which would allow for another direct comparison with lattice measurements, as has been outlined in ref. [11]. Finally, it would be interesting to consider non-equilibrium backgrounds, probably relevant for practical heavy ion collision experiments. We hope to return to some of these challenges in the future.

Acknowledgments
M.L acknowledges useful discussions with C. Gale, J. Ghiglieri and G. Vujanovic. This work was partly supported by the Swiss National Science Foundation (SNF) grant 200021-140234.

A Free limit of LPM resummation
As a crosscheck, we discuss here what happens if the potential iV + is replaced by i0 + in eq. (4.20). The correctly normalized regular solutions become (for 0 < ω < k 0 ) Im[g(y)] π = ω 1 ω 2 2πk 0 Im