Improved opacity expansion at NNLO for medium induced gluon radiation

When an energetic parton propagates in a hot and dense QCD medium it loses energy by elastic scatterings or by medium-induced gluon radiation. The gluon radiation spectrum is suppressed at high frequency due to the LPM effect and encompasses two regimes that are known analytically: at high frequencies $\omega>\omega_c = \hat q L^2$, where $\hat q $ is the jet quenching transport coefficient and $L$ the length of the medium, the spectrum is dominated by a single hard scattering, whereas the regime $\omega<\omega_c$ is dominated by multiple low momentum transfers. In this paper, we extend a recent approach (dubbed the Improved Opacity Expansion (IOE)), which allows an analytic (and systematic) treatment beyond the multiple soft scattering approximation, matching this result with the single hard emission spectrum. We calculate in particular the NNLO correction analytically and numerically and show that it is strongly suppressed compared to the NLO indicating a fast convergence of the IOE scheme and thus, we conclude that it is sufficient to truncate the series at NLO. We also propose a prescription to compare the GW and the HTL potentials and relate their parameters for future phenomenological works.


Introduction
The strong modification of jet observables in Heavy Ion collisions (measured both at RHIC [1,2] and the LHC [3][4][5]) when compared to proton-proton events, provides one of the key observations of the formation of the Quark Gluon Plasma (QGP) in such events. The continuous interactions between these hard probes and the dense QCD plasma induce a cascade of gluons, which inevitably modify the jet's properties (see [6,7] for recent reviews on the topic). As a consequence, for extracting the QGP properties from the experimental study of jets, an accurate and complete understanding of the medium induced radiation spectrum is critical.
One of the first (and most crucial) theoretical steps towards this goal consisted on the study of the emission spectrum of a single parton embedded in a QCD medium. In the regime where the medium is sufficiently large such that the parton may interact with several scattering centers in the plasma, the medium induced spectrum admits a full analytic treatment, captured by the Baier-Dokshitzer-Mueller-Peigné-Schiff-Zakharov (BDMPS-Z 1 ) formalism [10][11][12][13][14]. The region of validity for the BDMPS-Z formalism is bounded from below by the single (low energy) scattering limit (Bethe-Heitler limit), where the quantum mechanical formation time of the radiated gluon is of the order of the in-medium mean free path, t f ≡ ω/k 2 ⊥ ∼ mfp which is assumed to be much smaller than the medium length L, and thus the gluon is emitted incoherently by individual scattering centers. On the opposite end, the gluon formation time is bounded from above by the medium length L. In this regime, multiple soft scattering may act coherently as a single scattering center during t f . Hence, the effective number of scattering centers is much smaller than the actual number of scattering centers, i.e., N eff ∼ L/t f L/ mfp ≡ N scatt . The transverse momentum accumulated during t f via diffusion, k 2 ⊥ ∼qt f , whereq is the corresponding transport coefficient. This allows us to solve for the formation time (1.1) Hence, the radiative spectrum is suppressed as N eff ∼ ω −1/2 for ω BH =q 2 mfp ω ω c = qL 2 . This is the QCD analog of the Laudau-Pomerantchuk-Migdal (LPM) effect [15,16]. For formation times larger than L a maximum LPM suppression is achieved.
The above parametric analysis is valid so long as one can neglect large momentum transfers from the medium to justify the application of the diffusion approximation. However, due to the large Coulomb tail in the elastic cross section the medium transport parameterq will depend logarithmically on the transverse size of the radiated gluon. In the BDMPS-Z approximation, one assumes thatq is roughly constant invoking the slow variation of the Coulomb Logarithm. This makes the problem analytically tractable, but fails to capture the correct scaling in the region of phase space where the dominant contribution comes from single hard scattering (which is correctly captured by the GLV approach).
Until recently, an analytic approach which was able to connect the BDMPS-Z and GLV regimes into a single framework was not known, although several numerical based approaches were able to solve the problem exactly [17][18][19]. In previous papers, one of us introduced a systematical way of taking into account the hard p T tail encapsulated in the medium scattering potential [20]. This approach was latter extended to also take into account the full scattering potential [21], and was shown to correctly capture the features of both regimes 2 . In this paper, we will refer to this approach as the Improved Opacity Expansion (IOE).
We extend the work presented in [20] by computing the next order contribution to the integrated medium induced emission spectrum in the IOE approach. We study the NNLO term (i.e. we allow for the possibility of two hard scattering centers) in the IOE, showing that this term gives a small contribution to the full spectrum (when compared with the LO (BDMPS-Z) and NLO terms), ensuring that the series expansion is under control.
We show, in particular that in contrast to the plain opacity expansion where high orders are suppressed by inverse powers of ω, which is indeed the case for ω > ω c , in the regime ω < ω c higher orders are only suppressed logarithmically and the leading order power scaling ω −1/2 extends to all orders. As a result, the full spectrum in this regime can be expressed in the leading order form with an effective transport coefficient that can be calculated order by order in the IOE scheme, that is, where the effective transport coefficient is calculated to NNLO in the IOÊ where the IR cut-off's, that are fully fixed at leading logarithm accuracy, read andq 0 is given by Eq. (2.5) and Eq. (2.9) for the GW and the HTL models respectively. The present manuscript is divided as follows. Section 2 and subsections therein review the work presented in [20] and introduce a general form for the IOE expansion. In section 3 we study the NNLO order term in the IOE. Finally we discuss and summarize our findings in section 4. Complementary numerical work is shown throughout the paper.

Medium-induced gluon spectrum
The general form for the integrated medium-induced gluon spectrum off a high energy parton (in color representation R) is given by [14,20] where ω is the gluon frequency (assumed to be much softer than the emitting parton E ω) and the second term inside the brackets subtracts the vacuum like contributions. The Green's functions K and K 0 are solutions to a 2-dimensional Schrödinger equation in the transverse plane, and obey where v(x) is the potential defined by the in-medium (elastic) scattering cross section. K 0 obeys a similar equation with v = 0.

The HTL and Gyulassy-Wang potentials
Typically, the in-medium scattering cross section can either be obtained from Hard-Thermal-Loop (HTL) theory or from the Gyulassy-Wang model (GW model). Since in this paper we only focus on the large k ⊥ tail corrections, to leading logarithmic accuracy the model choice is irrelevant. However, even at leading order both models differ by finite terms. As such, the choice for the medium parameters, in each model, has to take this into account by providing a map between the different model parameters and the set of physical parameters. Therefore, before presenting the Improved Opacity Expansion in order to solve (2.1), we compute the potential v entering (2.2) to leading order accuracy in both the GW and HTL models. As we will show, this will allow not only to have the full leading term of the potential in both models, but also provides a map between each model. The elastic cross section in the GW model corresponds to an Yukawa interaction and reads where n is the density of scattering centers in the medium, µ is an infrared cut-off related to the Debye mass in the plasma m D and g is the QCD coupling constant. In general, n is a function of time, which we will overlook for the moment. Then the potential v appearing in where we have explicitly introduced the color charge C A directly into the potential and K 1 is the modified Bessel function of the second kind of order 1. Following [20] we have introduced the transport coefficient stripped of any logarithmq 0 where we have introduced the physical scale µ related to the GW screening mass µ ≡ µ GW by µ 2 ≈ µ 2 GW 0.29. The subdominant term is suppressed by a power of x 4 µ 2 . The HTL formalism [23] predicts an elastic cross section of the form where m 2 D = (1 + n f 6 )g 2 T 2 is the QCD Debye mass (squared), T is the QCD plasma temperature and n f is the number of light quark degrees of freedom. For an equilibrated system it is well known that n ∼ T 3 and the HTL and GW give the same result when the Debye mass is only taken into account as the infrared cut-off.
Solving for v HTL analytically is similar to v GW (see Appendix A). We find v(x, t) HTL = where nowq 0 is given byq (2.10) This provides a complete and consistent (leading order) map between the set of parameters used in the GW model and the set of parameters coming from the HTL framework. In particular if we want to match the HTL and GW models to logarithmic accuracy we should have m 2 D = e µ 2 GW . (2.11) In conclusion, we have shown that the leading logarithmic behavior for the potential entering equation (2.2) is fully captured in the GW model and HTL approach by defining a physical screening mass µ and a map between the medium model parameters and this scale. This is clearly seen in figure 1, where we computed the full HTL and GW potentials with the prescription given by (2.11). It is clear from this numerical exercise, that up to around |x| ∼ 2 m D , the two potentials match almost exactly. As one goes to larger dipoles sizes, the leading logarithmic approximation is not enough and a new map taking into account higher order terms would have to be constructed 3 . Therefore, for the rest of this paper we will work with the the leading order accuracy cross section given by The propagator potential (at leading logarithmic accuracy) then reads whereq is the medium transport coefficient and µ is introduced as an infrared cut-off. Note that in generalq has a trivial dependence on time via n, which in this paper we take to be n(t) = nΘ(L − t), with L the medium length (plasma brick model).

The harmonic oscillator approximation
In general the solution to equation (2.2) is not known in a closed form even for the leading logarithmic potential Eq. (2.13), but the case of vacuum propagation and as we shall see shortly, the harmonic oscillator, admit a complete full analytic treatment. The free propagator reads (2.14) To connect Eq. (2.13) with the harmonic oscillator, following the BDMPS-Z approach, one assumes that 1/x 2 ∼ Q 2 , where Q 2 is the typical transverse scale of the process to be determined later. This makes the potential that of a harmonic oscillator (HO). Therefore a closed form solution for K ≡ K HO exists [24] K HO (x, t|y, The functions S and C are the solutions to the initial condition problems (2.16) Here we have defined the (complex) frequency of the harmonic oscillator In the case where the medium is a plasma brick we have the closed form solutions for C and S given by The time dependence on Ω has disappeared since the medium is assumed to be static. In general, the functions C and S can be obtained for any medium, by either solving the path integral K via the semi-classical approximation 4 (Pauli's formula) [25] or by following the Wronskian approach (i.e. solving (2.16); see [20,26] for details). For the present paper, the second approach is more useful.
The properties of the functions C and S can be studied in general for any medium (see [26]). In particular, one can show that they obey the following identity [20,21,26] where we also introduced the handy notation C(t, s) ≡ C t,s and equivalently for S. In the future we will extend this notation to allow the shorthand writing of C(t 2 , t 1 ) ≡ C 2,1 , with the same applying to S and only valid when the dependency in t is clear and only the sub-indices matter. Although not generally true, we will also treat C has being an even function, which is true in the plasma brick model.

The general structure of the spectrum
The Improved Opacity Expansion is realized by expanding the full medium induced spectrum around the BDMPS-Z solution, such that the leading order (LO) term in the expansion matches the known solution and the higher order (N m LO) terms take into account the hard scattering contributions. This is achieved by rewriting the potential v as [20] v(x, t) = 1 4 (2.21) 4 In fact, for quadratic Lagrangians this approximation gives the exact result [24].
Here Q 2 is the matching scale between the two regimes of the spectrum and as it is clear from (2.21), the spectrum is independent of it, when all orders in perturbation theory are taken into account. For the moment Q 2 is assumed to be arbitrary, but as we shall see the logarithmic structure of the expansion requires a specific choice that is as expected the typical transverse momentum acquired by the radiated gluon Q 2 ∼ √ ωq. This expansion can be incorporated into K by using the Dyson-like equation for the propagator where´z ≡´d 2 z. Each order in perturbation theory is then obtained by expanding the above equation in powers ofq 0 (see the detailed discussion in [20]). Doing this procedure and using equation (2.1) the full spectrum reads 5 [20,21]

The leading order (BDMPS-Z) term
The leading order term is the well known BDMPS-Z result that can be obtained by using equations (2.15) and (2.14) in the general formula for the spectrum (2.1). One then obtains the compact formula [20] ω where from this point on we always assume that we are working within the plasma brick model. It is then easy to show that the S and C functions obey the following differential .
This allows one to make the t 2 integration directly and obtain where in the first step we cancelled the divergent pieces between K HO and K 0 and we have used (2.20) in the last step. Then the spectrum finally reads Defining the characteristic frequency ω c as the typical frequency of the emitted gluon with formation time of the order of the medium length we can obtain the asymptotics of (2.27) which quantitatively shows the scalings discussed in section 1. As we shall see the HO does not correctly capture the scaling when ω ω c , i.e., ω −1 , which is dominated by a single hard scattering.

The m th order correction
The general form for the m th contribution to the full spectrum which includes the hard scattering potential is given by (2.30) Here we have ordered the times of each scattering center from t 1 to t 2 in increasing order of the sub-index, running from 1 to m. The transverse position of the i th scattering center z i is also ordered from the first scattering center (z 1 ) to the last one (z m ). Equation (2.30) is obtained by iteratively using equation (2.22) in equation (2.1). As was shown in [20], the two extreme propagators can be integrated out, after performing the derivatives and using the general relation We are then left with just the intermediate position integrals and the time integrations at each scattering center. Introducing the representation for K HO in (2.15) and using the explicit formula for δv one eventually obtains the compact formula where we use the prescriptions: s 0 = 0, σ m+1,m = 1 and z m+1 = 0. Also, the factor depending on z m and z 1 outside the product, should be understood to be integrated over (i.e. the factor enters the z 1 and z m integrals; this is denoted by the slashed integral symbol). Here the factor π m comes from the m factors of K HO present in the general formula and the factorq m 0 is due to the presence of m δv terms. The 2 3m appears as a combination of the K HO normalisation factors and the terms in δv.
We have introduced the following functions with the boundary properties C 1,0 = C ∞,1 and C m+1,m = C m,0 and the same for the S function. Also It is clear from equation (2.32) that performing the remaining integrations is non trivial when m ≥ 3, so that the NLO and NNLO are special cases where one can hope to make analytical simplifications. For completeness we give the NLO (m = 1) term, already computed in [20,21] where in the second step we have translated from the notation for general m to the case m = 1 and we used the boundary properties of the C function. This result perfectly matches the result from the previous papers 7 .

The Next-to-Leading order correction
Before computing the NNLO term in the IOE, we present the NLO contribution already computed in previous work [20]. Starting from (2.35) we use the identitŷ to get the spectrum where the angular integration was also carried out. 7 We would like to point out that in equation (2.36) there is an overall extra minus sign when compared to [20,21]. This corrects the small mistake present previously, which does not affect the results significantly.
In analogy to what was done for the LO term, we also study the limiting cases ω → 0 and ω → ∞.
In the first case, it is easy to check that k 2 (s) → −ωΩ. The NLO contribution can then be computed exactly [20] and reads which shows that at the low frequency part of the spectrum this contribution scales like the LO term but suppressed by a logarithmic contribution ∼ log −1 √q ω µ 2 . To get to this result we have assumed in the last step that Q 2 ∼ √q ω. On the other hand, the high energy limit implies that k 2 (s) → iω 2s . The NLO term becomes dominant in this region of phase space and we have from equation (2.38) which matches the asymptotic behavior of GLV [9,20]. Here we usedω c ≡ µ 2 L 2 and χ ≡q 0 L µ 2 that measures the opacity. This term is dominant compared to LO contribution (the BDMPS-Z result is power suppressed).
3 The Next-to-Next-to-Leading order correction Using (2.32) we can obtain the NNLO term directly. The angular integrations can be done in a straightforward way and in the end one is left with 4 integrations to perform.
where J 1 is the Bessel of the first kind of degree 1 and´z ≡´∞ 0 dz. From this point on the indices in σ will be dropped. We have Formally, it is still possible to further simplify (3.1) by performing one of the z integrations. This leads to the appearance of a finite sum of Bessel functions and a logarithmic contribution. Although this decreases the number of integrations by one, the result obtained is neither suitable for numerical implementation nor is it of easy analytic study. We proceed to analyze equation (3.1) in two limiting regimes. First we explore the region where ω → ∞, i.e. where the major contribution to the spectrum should come from the NLO term. Then we study the opposite limit where ω → 0.

Large frequency limit
In this regime we let ω → ∞ ≡ Ω → (1 − i) × 0. We notice that in such regime the k 1 , k 2 and σ functions can be simplified using the fact that C sa,s b → 1 and S sa, Before using this approximation, we rewrite the spectrum in such a way that the z 1 , z 2 can be integrated out. We proceed by power expanding the Bessel function This representation is advantageous since it allows to directly do the integrations in z 1 and z 2 by usinĝ where Eψ(n) = exp(ψ(n)), ψ(n) = Γ (n)/Γ(n) and Γ is the gamma function. Putting all together the spectrum then reads We see that the converge of the series is then controlled by the dimensionless quantity . We notice that in the high energy limit using (3.3), equation (3.8) reduces to which is a well posed expansion parameter since s 2 < s 1 . Therefore, in the high energy limit, the power expansion representation of the spectrum gives a convergent series, and therefore only a finite number of terms in the expansion are required to achieve a reasonable numerical convergence.
Before numerically solving equation (3.6), let us study the asymptotic behavior of the spectrum analytically. The LO term for ω ω c scales withᾱ( ωc ω ) 2 (see (2.29)). The NLO term is given by (2.40) and scales withω c ω and is therefore the leading order term. The NNLO term is better discussed in terms of the rate ω dI NNLO dLdω . Then in the same limit as before the contribution reads .
Taking the real part we have (3.10) To proceed we rescale the time integration with u = s 2 /L and only keep the leading order contribution in the logarithms ∼ log( ω Q 2 L ).
where we have neglected all terms not doubly enhanced by logarithms and the remaining (finite) numerical factor coming from the integration in u. In the last step the u dependence in the logarithms can be dropped since it is only single logarithmic enhanced. The full spectrum predicted by the IOE is then dominated by the NLO term to all orders, since all higher order terms contribute with power corrections ∼ω c ω , which are suppressed. There are also some logarithmic enhancements, but this are always small compared to the power terms. Notice that, when moving away from the strict high energy limit, the NLO (i.e. leading term) will originate corrections (coming from the k 2 expansion) which contribute at LO order ∼ ωc ω 2 , with some possible logarithmic corrections. This also applies to higher contributions, where the corrections coming from the high energy limit expansion of the k's and σ's functions come with extra power law contributions. This fact ensures that the NLO term will always be the dominant piece in the IOE.
In figure 2 (Left) we present the numerical computation of the LO, NLO (already shown in [20]) and the NNLO terms in the IOE. In addition, we present the GLV spectrum. The numerical results depict exactly what was argued before. At large ω, the NLO term becomes the dominant contribution to the spectrum. The NNLO at LO lines become almost parallel at large ω, thus showing that these two terms give the same asymptotic contribution (this is not strictly true, since they will differ by subleading logarithmic terms). Moreover, we also notice that the actual numerical values assumed by the NNLO curve are at their best only an order of magnitude smaller than the NLO contribution. This shows, that for practical purposes, in this regime, the NLO truncation already offers an excellent approximation to the full spectrum, and subleading corrections do not change the behavior of the IOE.

Small frequency limit
The small frequency regime requires a more delicate approach. This is mainly due to the fact that in this limit the BDMPS-Z solution, without any kinematic constraints [27,28], is divergent. In the case of the LO and NLO, this divergence is well under control, since the diverging pieces factorize from the remaining terms. This is no longer true at NNLO order, and thus requires a more careful treatment.
Our starting point is again equation (3.1), but we now take the limit ω → 0 ≡ Ω → (1 − i) × ∞. From the discussion present in the last subsection, it is clear that in this case using the power expansion of J 1 directly is not an optimal strategy, since at some point we would be required to resum all terms in these expansion. Therefore, we keep the integrations in z 1 and z 2 , and take the limiting forms for the C and S functions lim ω→0 Ω cos(Ωx) sin(Ωx) = iΩ , lim We apply this approximation in all the C and S terms but the ones that explicitly dependent on the time difference s 1 − s 2 . In such terms, we cannot use the above approximation 8 since Ω ∼ 1 t f , where t f is the typical formation of a BDMPS-Z gluon. Since, parametically, the support of the functions depending on the time difference s 1 − s 2 is of order t f , these type of dependencies have to be kept in full. If neglected, we would be ignoring the part of the support of the function where it is not damped or highly oscillatory. Then one gets and σ cannot be simplified. The NNLO contribution to the IOE spectrum reads (3.14) To proceed, we do the change of variables (s 1 , s 2 ) → (s 1 , τ = s 1 − s 2 ). To continue, we notice, as argued before, that the main contribution to the integral comes from the region τ ∼ t f , and therefore, the dependence of the result on the upper bound of the integral is small. Therefore, we promote the upper bound L − s 1 → ∞. The integration in s 1 is then trivial and we are left with just one time integration. This approximation is similar to approaches where the medium induced gluon emission is taken in the Markovian (classical) limit, where the all the shower is dominated by decoherent emissions [29].
In this regime, we can rescale τ → t = q 4ω τ , so that the integration is done in terms of dimensionless quantities. Finally, the functions C and S still present in k 1 , k 2 and σ have the complex argument (1 − i)t. Therefore, we Wick rotate the time integration with the transformation −iT = (1 − i)t. Then the result reads To make the integral completely dimensionless, we take the scale Q 2 ∼ √q ω as in [20]. Then the remaining integral can be computed exactly (3.16) 8 In case this was done, the result obtained would be divergent.
Thus the scaling for the NNLO term at small frequencies reads The NLO contribution is given by equation (2.39) and exhibits the same scaling when Q 2 = √ ωq. Unlike the high energy limit, where we observed that moving away from the strict ω → ∞ limit originated terms which have to be incorporated in the all orders expansion, at small frequencies (and evaluating Q 2 = √ ωq ≡ Q 2 c ) an all orders expansion can be written and reads where the spectrum at LO is understood to be taken in the small frequency regime and the coefficients c 0,0 ≡ 1, c 1,0 , c 2,0 , · · · , are pure real numbers, computable order by order. Notice that this expression is consistent at all orders, since every term exhibits the same scaling, up to logarithmic enhancements. The sub-indices of the c coefficients comprise two numbers, the first indicating the power of log Q 2 c µ 2 in the expansion. The role of the second index will shortly become evident.
The main difference between (3.2) and the high energy scaling is that in this case the subleading terms in the expansion (in which the LO contribution is the leading term) can become large, due to logarithmic enhancements. In the high energy limit, this was not possible since each sub leading contribution was power suppressed and the logarithmic enhancements were not dominant.
For a general choice of scale Q 2 equation (3.2) can be written as where the LO spectrum is taken at a general scale Q 2 . In the last line, we have introduced the function W which formally resums the all-order terms. Notice, that although order by order W exhibits a dependence on Q 2 the all order result is independent of the choice of the matching scale. We will discuss the properties of W later on.
We see that the second index in the c i,j coefficients denotes the expansion in powers of log j √ ωq Q 2 , opposed to the first index which denotes the terms proportional to powers of log −i Q 2 µ 2 . We have computed the coefficient c 1,1 = 1 2 explicitly in section 2.6. It is interesting to note the role that the matching scale plays in (3.19). First, suppose that we fix the matching scaling at some constant value Q 2 ≡q 0 L, which is a higher momentum scale at ω ∼ ω c . Then the logarithms scaling with Q 2 µ 2 are fixed and the evolution with ω is encoded in the logarithms of √ ωq Q 2 . This implies, that at small ω, there is a breakdown of the series since while the LO contribution diverges with ∼ ω c /ω the N m LO contributions diverge (the most diverging piece) ∼ log m √ ωq Q 2 , where we have neglected (for this discussion) the different power of logarithms in the denominators since they are constant. In fact, we expect that when ω ∼q 0 L 2 / log 2 Q 2 µ 2 the expansion breaks down. Note that this scale is parametrically much larger than ω BH ∼q 2 mfp , and thus while the LO term gives a constant contribution all other orders strongly diverge. This clearly shows that the matching scale has to be chosen such that there is a correct interpolation between the GLV and BDMPS-Z limit, which implies that Q 2 ≡ Q 2 c (ω) ∼ √ ωq. This choice will allow for mutual cancellations between the different orders in the IOE so that the spectrum does not depend on the matching scale when all orders are resumed and the correct spectrum is recovered (while still away from the Bethe-Heitler limit).
From (3.19), it seems that the N m+1 LO contributions can impact the terms at order N m LO. However, let us suppose that we choose a scale Q 2 ≡ a 2 Q 2 c , where a is dimensionless factor that rescales Q 2 c ∼ √ ωq. Then, to leading logarithmic accuracy, (3.19) becomes where we neglected the dependency in a in the logarithmic correction in the NLO term, since it is easily seen that it only contributes at higher orders. Thus, different choices for the multiplicative factor of Q 2 , at NLO accuracy, only give rise to higher order logarithmic corrections. This observation has to hold to all orders in perturbation theory, since when all terms in the series are resumed, the spectrum is independent of the choice made for the matching scale. Therefore, the expansions differing in the choice of the matching scale and truncated at some order, can only differ by higher order corrections. This fact also allows to reduce the number of independent coefficients c i,j to be computed.
In summary, we have shown that not only one has to allow for a dependence on ω in the matching scale for the perturbative expansion to be meaningful, the natural choice for this scale is Q 2 ∼ Q 2 c ≡ √ ωq, and other choices for the matching scale only differ by subleading factors (assuming one uses Q 2 ∼ √ ω, which is the only physically reasonable scaling law for this problem).
In Figure 3 we compute the spectrum at NLO accuracy while fixing the scale Q 2 ≡ Q 2 c (see figure for details) and then varying it by factors of 2. We clearly see that the variation in the matching scale, in the low energy regime, lead to minimal modifications of the spectrum. In fact, we can see that this happens because while Q 2 increases the LO contribution increases but the NLO term becomes smaller, such that the contributions balance each other out, as shown in the analytical study. We would like to point out that this study is distinct from the one perform in [20] where one varied µ 2 . It is easy to see from the above expressions, that this does not lead to the same evolution as the one presented here. Figure 3 also shows that unlike the case where Q 2 is a fixed scale here, the spectrum does not diverge around ω ∼q 0 L 2 log 2 Q 2 µ 2 and the LO and subleading terms balance each other out. This clearly shows that the interpolation problem between the GLV and BDMPS-Z regimes requires a non trivial fix to how one defines the matching scale between the soft and hard regimes.
Finally, it is also interesting to know how close to µ 2 the matching scale Q 2 , that decreases with ω, can get, such that the NNLO is still significantly smaller than the LO and NLO terms. This translates into the sensitivity of the IOE to the approach to ω BH ∼ µ 4 /q. It is clear that when Q 2 → µ 2 (or equivalently ω → ω BH ) every term in the expansion diverges. Starting from equation (3.2), we normalize the full spectrum in the low energy regime to the LO result and obtain lim ω→0 ω dI dω norm. = 1 + 0.508 where β = log Q 2 c /µ 2 = 1 2 log(ω/ω BH ). If we want to compare the contribution of the NNLO versus LO+NLO we just need to compute where α ≡ NNLO 1+NLO gives the percentile contribution of the NNLO term compared to the LO+NLO (up to NNNLO corrections).
To proceed, we wish to discuss this scaling in terms of the Bethe-Heitler frequency ω BH ≡ µ 4 /q 0 (recall Q 2 /µ 2 = ω/ω BH ). Then we can rewrite (3.22) as where we have chosen the positive root since it is the one of physical relevance.
We then have that for α = 1%, ω ≥ 18.83 ω BH ; α = 10%, ω ≥ 1.98 ω BH and for α = 50%, ω ≥ 1.21 ω BH . The inequality symbol comes from the fact that the above equation gives the lower limit for ω below which the ratio NNLO/(1 + NLO) exceeds the value of α. We see that the evolution with α is quite fast: when one requires α ∼ 1% the limit frequency has to be one order of magnitude larger than ω BH , but when α ∼ 10% the limit frequency is of the order of ω BH . This shows that for the NNLO terms to be negligible (say giving less than 10% of the total contribution to the spectrum) compared to the LO and NLO terms, is not strongly dependent on low momentum tail and any typical energy scale would satisfy the inequalities presented above. Conversely, choosing matching scales which are essentially of the order of the Bethe-Heitler scale leads to the breakdown of the perturbative expansion, as expected (notice that when ω = ω BH , equation (3.2) becomes meaningless).

Discussion and Outlook
In this paper we have provided an analytical and numerical study of the Improved Opacity Expansion at up to NNLO accuracy. In addition, we have presented, for the first time, a map between the GW and HTL models for the elastic in-medium cross section and the set of physical parameters at leading logarithmic accuracy. This results are best summarized in figure 1, where it is clear that using our map with the GW full potential gives back an extremely good approximation of the HTL potential, up to dipole sizes |x| ∼ 2/m D . The combination of both these results, guarantees that we have a complete and systematic control over the analytic structure of the emission spectrum (2.1). This mapping is crucial since it gives meaning to comparisons between emission spectrums using different medium models. In the particular case of the IOE, we showed that this allowed us to have a full control over the accuracy of our result.
Moving on to a more detailed discussion of the IOE, our study allowed us to show that in the large frequency domain the spectrum is strongly dominated by the NLO term, which follows the well known GLV scaling. All other orders in the expansion, are power suppressed by factor ofω c /ω. In particular we showed that the LO and NNLO terms are of the same order. However, an all order closed form formula is not possible to write down since, as argued above, as one moves away from the strict high energy limit, new contributions appear which are not power suppressed.
On the opposite end of the spectrum, we found that the IOE has an extremely rich and interesting structure. In this limit, the LO term is the dominant contribution to the expansion, but important logarithmic contributions appear, order by order. This is in opposition to the high frequency regime, where NLO term is dominant over power suppressed contributions.
In order to better understand the structure of the IOE in the small frequency domain we first noticed that for a fixed matching scale the expansion is ill defined and this lead us to conclude that there exists a natural scale Q 2 c = √q ω which guarantees that the ω dependency of the matching is such that mutual cancellation between the many orders of the IOE guarantee that the full spectrum is finite. In addition, we showed that rescalings of Q 2 c only affect higher order terms in the IOE (see (3.20)). This was numerically confirmed by the results in figure 3. Additionally, we want to point out that the exercise shown in figure 3 clearly demonstrates that the interpolation between the GLV and BDMPS-Z regimes requires a proper treatment as the one provided by the IOE, and does not allow for a simplistic interpolating procedure. In fact, we have shown that the correct contribution to the spectrum in the region ω < ω c needs both the LO and the NLO terms in order to describe the correct result Both these results are a direct consequence of the fact that the spectrum's dependence on matching scale must vanish when all orders in the IOE are taken into account. In fact, this observation means that after all terms are taken into account the spectrum must be of the form (for a general Q 2 scale; see equation (3.19)) where W is a general (unknown) function (introduced in (3.19)), which captures all the finite corrections to the spectrum. Notice that the dependency in Q 2 disappears. From equation (3.2) we can construct the W function order by order as where we have chosen the scale Q 2 = Q 2 c . It is straightforward to obtain to corresponding expansion of W , These results show that the IOE admits to be written in a simple closed form for a fixed accuracy level with an additional prescription for the matching scale. All the results are valid so long as the matching scale is chosen sufficiently higher than the Bethe-Heitler scale ω BH .
Before moving on, we wish to point out that in (4.1), although the leading logarithmic behavior between each order truncation is well under control, there are logarithmic contributions order by order which might spoil the behavior of the series. Recall from above, we first showed that to have a proper converging series one has to require the matching scale to evolve with ω and then we showed that there is a natural choice for this scale, with other choices (with the same scaling) varying only by subleading terms. However, before we ignored that when varying the scale Q 2 c subleading terms (like log log Q 2 c µ 2 ) can be subleading in the number of logs but be of the order O(1) 9 . This is a direct consequence of the fact that the matching scale is given by a recursive equation and one has to expand the recursive equation to a certain degree of accuracy. Then in the regime ω BH ω ω c , 9 For instance, from equation (3.19), when expanding the LO term we obtain the leading contribution ∼ log ω ω BH , while the subleading term reads ∼ log log ω ω BH . Therefore, normalizing to the LO term, the subleading term can contribute at NLO order (i.e. when counting the denominator logarithms) and can be an important factor since log log might be of order of the leading coefficient c1,0. This discussion follows to all orders and is a direct consequence of the fact that the matching scale is defined by a recursive equation.
the full spectrum should read  [20,21]. Then we have the remarkable result that in the small frequency regime, the full spectrum is captured by the BDMPS-Z solution with a renormalization ofq. Notice, that the above discussion where Q 2 c ∼ √ ωq 0 still holds when comparing the different orders of the IOE, but they might fail due to log log contributions coming from the definition of the matching scale. Again, this exercise explicitly shows that the definition of the matching scale between the GLV and BDMPS-Z is a non-trivial problem and it can not be simply fixed ad-hoc.
Another important point in the work presented in this paper, is the evidence that the contributions to the spectrum coming from the NNLO order correction are parametric and numerically small. This ensures that, for example, in phenomenological applications, the LO+NLO truncation is sufficient. We have thus shown that the IOE provides a complete, systematic and self-consistent interpolation procedure between the GLV and BDMPS-Z pictures and this results holds as long as ω ω BH 11 . Figure 4 explicitly shows that including the NNLO term does not give any significant correction to the full spectrum. In this plot we have extended the small frequency regime to large energies via the LO scaling at low energies: at intermediate ω the NNLO spectrum is obtained by extending the low energy result up to a matching scale, after which the high energy evolution at NNLO is used. We have tested this numerical procedure for several choices parameters L, µ andq 0 and verified that, for reasonable matching scales ω match ∼ ω c , the two ends of the spectrum nicely match each other. In addition, for several choices of parameters we have also seen that the NNLO contribution is always much smaller than the LO and NLO terms.
This work ensures that for future endeavours, it is sufficient to just keep track of the LO and NLO terms of the IOE. Therefore, taking into account all the results presented, in the future, we will be able to explore the single inclusive emission spectrum at NLO 10 Notice that this truncation includes the first double logarithmic contribution, as discussed in the previous footnote. 11 Extrapolating to near the scale ωBH has already been studied at NLO accuracy [21], although many questions are still to be answered. In addition, we also include the spectrum up to NLO (black line) and NNLO (pink line) accuracy. All the curves are the same as in the previous plot, except the NNLO solution which is obtained by matching the scaling as ω → 0 at a frequency cut off ω cut = 0.5 × ω c , after which the solution is obtained by using equation (3.6). The scaling for smaller frequencies is obtained by making use of equation (3.16) and the LO scaling law ωc ω . This procedure is indicated by the matched tag.
accuracy, while being able to have full control over the accuracy of the result. This is a key step for phenomenological implementations of the IOE.

A Useful integrals
In this appendix we shall calculate the following integral that is related to the GW and HTL models by letting b = a = µ and b = 0, a = m D , respectively. First we decompose the integrant as followŝ (1 − J 0 (ux)) .  (1 − J 0 (ux)) = 1 (a 2 − b 2 ) [K 0 (ax) − K 0 (bx) + log a − log b] . (A.5) There are two special cases that will correspond to the two models of interest. First, a = b