Medium-induced radiative kernel with the Improved Opacity Expansion

We calculate the fully differential medium-induced radiative spectrum at next-to-leading order (NLO) accuracy within the Improved Opacity Expansion (IOE) framework. This scheme allows us to gain analytical control of the radiative spectrum at low and high gluon frequencies simultaneously. The high frequency regime can be obtained in the standard opacity expansion framework in which the resulting power series diverges at the characteristic frequency $\omega_c\sim \hat q L^2$. In the IOE, all orders in opacity are resumed systematically below $\omega_c$ yielding an asymptotic series controlled by logarithmically suppressed remainders down to the thermal scale $T \ll \omega_c$, while matching the opacity expansion at high frequency. Furthermore, we demonstrate that the IOE at NLO accuracy reproduces the characteristic Coulomb tail of the single hard scattering contribution as well as the Gaussian distribution resulting from multiple soft momentum exchanges. Finally, we compare our analytic scheme with a recent numerical solution, that includes a full resummation of multiple scatterings, for LHC-inspired medium parameters. We find a very good agreement both at low and high frequencies showcasing the performance of the IOE which provides for the first time accurate analytic formulas for radiative energy loss in the relevant perturbative kinematic regimes for dense media.


Introduction
High-energy collisions of heavy nuclei provide the necessary conditions for creating an extended medium of hot and dense nuclear matter, referred to as quark-gluon plasma (QGP). The appearance of a short-lived stage of this exotic state of matter leaves a strong imprint on particle production at all momentum scales and is quantified by high-precision experimental measurements. In this work, particular attention is devoted to the high-energy particles that traverse the medium and can be used as perturbatively well controlled probes -1 -of the microscopic properties of the QGP [1,2]. In terms of experimental observables, these objects emerge in the detectors as collimated sprays of particles and energy, colloquially referred to as jets. Jet modifications, quantified with respect to a baseline obtained in proton-proton collisions (or "vacuum"), have been intensively studied at RHIC [3][4][5] and LHC [6][7][8][9][10][11][12] since more than a decade.
The suppression and modification of jets produced in heavy ion collisions, commonly known as jet quenching, is driven by two main phenomena: transverse momentum broadening and energy loss. The former refers to the acquisition of transverse momentum by the highly energetic partons that make up the jet through elastic interactions with the medium following Brownian motion. This diffusion in momentum space is characterized by the transport coefficient:q where k 2 ⊥ typ is the typical squared transverse momentum transfer. An important role is also played by induced energy loss, such as caused by drag and inelastic, or radiative, processes. The latter component is a result of bremsstrahlung radiation triggered by collisions with medium constituents. The associated mean energy loss was found to scale as L 2 , where L is the length of the plasma, and thus constitutes an important, if not the dominant, source of jet quenching for large media, even though it is "naively" suppressed by a power of the coupling constant [13].
The above physical picture applies for the regime of multiple scattering during the passage through the medium. This is the case when the opacity χ = L/ mfp , defined as the ratio between the medium length, L, and the mean free path, mfp = (ρσ el ) −1 (where ρ is the density of scattering centers and σ el the total elastic cross section), is of order one or larger, i.e. χ 1. In effect, many interactions, i.e. those that occur within the formation time, coherently participate in inducing gluon radiation, in an analogous way to the well-known Landau-Pomeranchuk-Migdal (LPM) effect in QED [14,15]. It was soon understood that neglecting these interference effects would fail in adequately describing the radiative spectrum in a substantial region of phase space in such media. Focusing on the diffusion approximation, and thereby neglecting the Coulomb tail that captures the physics of rare hard momentum transfers, the radiative spectrum could be analytically computed albeit with an ambiguity in setting the upper bound for the Coulomb logarithm inq [13,[16][17][18][19][20].
More precise calculations can be performed in the dilute regime, where at most a few scatterings contribute and the full parton-medium interaction potential can be used [21,22,[22][23][24][25][26]. However, their domain of applicability is limited to low opacity χ 1 or large momentum transfer. Due to the finite radius of convergence of the opacity expansion such an approach diverges [27] and thus cannot be applied in the case of a dense medium.
As we have mentioned, in both low and high opacity regimes, under a set of physically motivated assumptions that will be revisited in what follows, the medium-induced spectrum can be computed analytically. However, the scenario explored in current experimental facilities, such as the LHC or RHIC, is expected to be the one where the jet undergoes a handful of scatterings, O (1 − 50), with the medium [28]. As such, neither of the above limiting forms is in its exact domain of applicability, thus hindering quantitative theoryto-data comparisons and the extraction of the medium parameters through the jet physics program in current colliders.
A first step towards a more precise control over the accuracy of analytic calculations was first taken in Ref. [49], where the next-to-leading logarithmic corrections toq were computed in the multiple soft scattering regime or, equivalently, in the infinite length medium limit. More recently, substantial progress was made in unifying the low and large opacity regimes while recovering the results of Ref. [49] in the soft regime. This new scheme, dubbed the Improved Opacity Expansion (IOE), has been shown to successfully meet this goal when applied to the description of the single particle broadening probability [50] and to the medium-induced gluon energy spectrum [51][52][53], which constitute the major tools used in a multitude of well established phenomenological models of jet quenching [54][55][56][57][58][59][60]. In a nutshell, this framework is built as a series expansion of the in-medium scattering cross section where the zeroth order term encodes the multiple scattering solution and higher N n LO orders 1 in the series account for n hard scatterings with the medium.
This paper aims at computing the fully differential medium-induced spectrum at NLO accuracy in the IOE framework. For different values of the gluon energy and transverse momentum, (ω, k), 2 we provide an analytic formula for an arbitrary medium profile that requires numerical integrations of the same order in computational complexity as the ones encountered in the multiple soft scattering limit [35,61]. In the simplified scenario in which the medium is treated as a brick of constant density, we have obtained closed, analytic formulas for the asymptotic behavior of the spectrum computed in the three different setups considered in this work: the IOE at NLO, the single-hard scattering approximation (SH) and the multiple soft scattering regime (MS). In addition, we make a phenomenologically oriented comparison with the all-orders numerical spectrum presented in Ref. [33]. The numerical routines used in this publication are provided as ancillary files that can be found in Ref. [62].
The remainder of this paper is structured as follows. Section 2 revisits the Improved Opacity Expansion framework in full generality, including its application to the single particle momentum broadening and the medium-induced energy spectrum calculations. The core of this paper is Section 3, where the fully differential spectrum is calculated with a high level of detail. Those readers more interested in the final result and not so much in the technicalities can find a summary of the ready-to-use formulas in Section 3.3. Fi- 1 The nomenclature used here to denote the orders in the IOE should not be confused with the more familiar perturbative expansion in powers of the coupling constant. 2 Bold letters denote 2D transverse vectors in this paper, while their modulus is written as |k| ≡ k ⊥ .
-3 -nally, numerical results for LHC-motivated medium parameters are presented in Section 4, including a comparison to BDMPS-Z, GLV and the resummed to all orders spectrum. Further details on the analytic calculations can be found in Appendices A, B, C and D.
2 The Improved Opacity Expansion: an overview The Improved Opacity Expansion draws on the seminal 1948 work by Molière [63], where the transverse momentum broadening of charged particles in QED was described in such a way that the multiple soft scattering solution, well described by a Gaussian distribution, and the Coulomb power-law tail were reproduced in the appropriate limits. The IOE program consists in extending Molière's original approach to QCD and to more complex observables. So far, this strategy has been successfully applied to compute the mediuminduced gluon energy spectrum [51][52][53] and the transverse momentum broadening distribution of an energetic parton propagating through a dense QCD plasma [50].
As an introduction to the IOE, it will be instructive to first revisit how it applies to the two aforementioned observables: transverse momentum broadening and medium-induced gluon radiative spectrum. This will also serve us to lay the basis of the calculation of the fully differential spectrum.

Transverse momentum broadening
The elementary in-medium process that underlies the observables that we discuss in this work is the elastic collision rate γ el ≡ dσ el /d 2 q, where q ≡ (q 1 , q 2 ) corresponds to the transverse momentum transfer in the t-channel between the hard probe and the medium. At leading order in the coupling the rate reads γ el ∼ g 4 n/q 4 ⊥ , where n corresponds to the density of scattering centers in the medium and the 1/q 4 ⊥ dependence denotes that, at short distances, the interaction is Coulomb-like. On the other hand, when q ⊥ → 0, the power law divergence should be screened by the medium at, roughly speaking, the Debye mass m 2 D in the plasma. Equipped with γ el , we can readily write a rate equation for the transverse momentum broadening distribution P(k; t), which gives the probability for a parton in color representation R to acquire transverse momentum k due to in-medium propagation during a time t, where the final time corresponds to the length of the medium t = L and C R is the color factor associated to a representation R of SU(3). 3 The boundary condition at initial time t = 0 is simply P(k, 0) = (2π) 2 δ (2) (k). The first term in Eq. (2.1) accounts for the gain in transverse momentum of the initial parton while the second term reflects the loss of probability for finding said parton with the measured momentum k. Notice that due 3 In this paper we use the notation´q =´d 2 q (2π) 2 to describe transverse momentum space integrals and x =´d 2 x for integration in position space.
-4 -to rotational symmetry the broadening probability is a function of the modulus of the transverse momentum vector, i.e. k ⊥ ≡ |k|.
The integral of the collision rate yields the inverse mean-free-path between two collisions, i.e −1 mfp ≡´q γ el (q). At low opacity, χ ∼ L/ mfp 1, the distribution is dominated by at most a single hard scattering (SH) and one finds Conversely, at high opacity, multiple (soft) scatterings occur with order one probability and Eq. (2.1) can be approximated by a diffusion equation for which analytic solutions exist. This is done by expanding in gradients for q ⊥ k ⊥ . The first non-vanishing contribution involves the jet quenching parameterq, where the integral over q is divergent in the ultraviolet and thus must be regulated, giving rise to the standard Coulomb logarithm, while the infrared region is cut-off by the screening mass that we denote by µ 2 * . Assumingq to be constant in time, the solution to the diffusion equation is a Gaussian [50] and the associated broadening distribution reads Although this result describes the physics of multiple soft scattering (MS) of the probe in the medium, the diffusion approximation has two major drawbacks: (i) it misses the heavy 1/q 4 ⊥ tail associated with large momentum exchanges and (ii) the transport coefficient depends, logarithmically, on an undetermined ultraviolet cutoff scale.
The IOE overcomes these two limitations by shifting the expansion point of the opacity scheme from the vacuum to the harmonic oscillator potential, resulting in the Gaussian distribution presented in Eq. (2.4). This shift in the expansion is easily performed in position space and thus we should consider the Fourier pair of P(k, t), (2.5) In position space, Eq. (2.1) becomes local where the scattering potential v(x) combines the gain and loss terms and is thus ultraviolet finite (2.7) In this example, Eq. (2.6) can be directly integrated, but this is not generally possible as we shall see in the case of the radiative spectrum. Furthermore, one still needs to invert the Fourier transform and this is where the IOE scheme will be particularly useful as it allows us to reduce the Fourier transform to a sum of standard integrals. Let us recall the main difference between the Improved Opacity Expansion procedure and the usual Opacity Expansion (OE) strategy [21][22][23]. The latter performs an expansion directly in powers of v(x) of Eq. (2.6), yielding a series in powers of q −2 ⊥ once introduced in Eq. (2.5), with the leading contribution given by Eq. (2.2). In the IOE, one shifts the expansion point to be a solution to Eq. (2.6) with the potential v = v HO ≡qx 2 ⊥ /4, whose Fourier transform can be carried out to yield Eq. (2.4). If we denote such a solution by P LO , then the aforementioned shift of the expansion point leads to Here the scattering potential is split into two terms, i.e. v = v HO + δv, such that |δv| |v HO |, in which case δv can be regarded as a perturbation around the potential v HO . In doing so we aim to tame the divergence of the plain Opacity Expansion series at low enough q, typically when q 2 ⊥ <qL. This separation of v into v HO and δv is in general arbitrary and requires the introduction of a matching scale Q. Clearly, truncating the IOE series at a fixed order introduces a residual dependence on the separation scale that is of the order of the remainder and thus can be safely neglected. It can nevertheless be used to gauge the uncertainty associated with the fixed order calculation very much like scale dependence encountered in standard perturbation theory calculations.
To illustrate this point consider the leading logarithmic form given in Eq. (2.7). One would trivially write ⊥ Q 2 , up to overall time-dependent factors. A natural candidate for the separation scale in the case of momentum broadening is Q 2 ∼ qt, corresponding to the average momentum squared accumulated by the probe due to multiple soft momentum exchanges with the medium. In general, when considering other observables, the LO provides guidance to what scale should be chosen for Q 2 . We shall see below how to make this observation more precise, in particular regarding the ultraviolet behavior ofq.
Not only the IOE allows to fix the diverging behavior of the Opacity Expansion, but also provides a good approximation of the exact result at low transverse momentum provided the following hierarchy of scales is met This ensures that the Coulomb logarithm is large, i.e. log(Q 2 /µ 2 * ) 1, and since at low k ⊥ , the x integral is dominated by the region x 2 ⊥ ∼ 1/Q 2 , we also have that log 1 On the other hand, at large k ⊥ the rapidly oscillating Fourier phase implies that x ⊥ k −1 ⊥ Q −1 which flips the relative order of the LO and its correction, i.e. |δv| |v HO |, and thus the logarithmic function in v(x) can no longer be neglected. Since large momentum transfers are associated with steeply falling cross sections, such a case is associated with -6 -rare hard scatterings in the medium, and perturbation theory is applicable, recovering the standard Opacity Expansion.
Following this more qualitative discussion that highlights the strengths of the IOE approach, let us make the discussion more quantitative and rigorous by recalling some of the results presented in Ref. [50]. First, in jet quenching phenomenology two models for the in-medium scattering rate are typically considered. One option, referred to as Gyulassy-Wang (GW) model [64], is to describe the medium as an ensemble of static scattering centers with Yukawa like potentials, with the in-medium rate given by where µ is the GW screening mass and n the density of scattering centers in the medium. This leads, see Eq. (2.7), to a scattering potential of the form where we have introduced the bare jet quenching parameterq 0 (t) = 4πα 2 s C A n(t) and K 1 is the modified Bessel function of the second kind of order 1. Another popular choice is to describe the medium as a thermal bath, so that the scattering potential can be perturbatively computed using Hard Thermal Loop (HTL) effective theory [65], with . (2.13) Here T is temperature of the medium and m D the Debye screening mass. The corresponding scattering potential reads where nowq 0 (t) = α s C A m 2 D (t)T , γ E = 0.577216 . . . is the Euler-Mascheroni constant and K 0 is the modified Bessel function of the second kind of order 0. 4 The differences and similarities between these two models have been extensively discussed in Refs. [50,53]. To leading logarithmic order, they can be unified in an universal form, in accordance with where µ * is a universal screening mass that can be mapped to the masses of both models considered above. 5 Applying the IOE prescription to split the potential as v = v HO + δv, see Eq. (2.9), and inserting it back into Eq. (2.6), we obtain, after expanding in powers of the perturbative 4 Here we defined the jet quenching parameter for gluons, i.e. CR = CA. 5 The GW mass µ is related to the universal mass µ * by 4µ 2 * = µ 2 e −1+2γ E , and the Debye mass mD in HTL corresponds to 4µ 2 -7 -potential δv, (2.16) where we identify the next-to-leading order (NLO) term with the contribution O(δv), the next-to-next-to-leading order (NNLO) with the O(δv 2 ) term, and so on. 6 In Eq. (2.16), we have introduced the bare saturation scale where we allow the bare jet quenching parameter to vary in time. In addition, we define the effective jet quenching parameterq(t) =q 0 (t) log Q 2 , where the logarithmic dependence appears naturally from the splitting of v(x). As discussed above, the definition of the matching scale, Q, can not be cast in a closed form, since it enters the definition ofq as well as depends on it directly. In turn, it is obtained by solving the transcendental equation where, following our previous reasoning, we have identified Q 2 ≡ Q 2 b with the effective saturation scale Q 2 s . We truncate Eq. (2.16) at NLO accuracy, since already at this order both the hard and soft regimes should be well described. The resulting broadening distribution reads [50] where x = k 2 ⊥ /Q 2 s , and is the expansion parameter of the series in the regime k 2 ⊥ Q 2 s . 7 At large momentum exchanges, k 2 ⊥ Q 2 s , one obtains from the NLO term, and in accordance with Eq. (2.2), that The series is truncated at nmax ∼ Q 2 s0 /µ 2 * since formally this is a divergent asymptotic series; the divergence is physically associated to the fact that x ⊥ can not be smaller than 1/µ * -see Ref. [66] for a further discussion on this truncation. 7 The exponential integral function is defined as -8 -while the LO term is exponentially suppressed. In this high momentum limit, we recover the Coulomb tail encoded in a single scattering in the medium. On the other end, when k 2 ⊥ Q 2 s , we find The first term corresponds to the LO contribution. Thus, the NLO term, up to a small constant logarithm, is of the same functional form as the LO but power suppressed by λ 1. In fact, one can show that, in this regime, perturbative corrections in the IOE scale as the LO term, each increasing order suppressed by an extra power of λ =q 0 q . Hence, in this limit the LO term dominates and one recovers the multiple soft solution, which correctly describes the physics at play.
In Fig. 1 we numerically compare the broadening distribution P(k, L), for a medium with constantq 0 , computed up to LO and NLO in the IOE, with the full P obtained using Eqs. (2.5), (2.6) and the GW potential in Eq. (2.12). The result follows the above discussion: at large momentum transfers, k 2 ⊥ qL, the NLO term dominates and converges to the full result dominated by the single hard scattering result (k −4 ⊥ ). On the other hand, at low momentum transfers the LO and LO+NLO become comparable, reproducing the full result within an uncertainty band associated to the remaining freedom in the definition of Q 2 b . The biggest mismatch between the LO+NLO result and the full distribution happens near the peak of the distribution and could be eventually improved by adding more orders in the series. Nonetheless, it is clear that the IOE approach provides a neat interpolation between the soft and hard regime, instead of properly describing just one of these regions.

The energy spectrum
As a second illustrative example, we consider the application of the IOE to compute the medium-induced gluon energy spectrum. The in-medium emission spectrum of a soft gluon with energy ω from a hard parton with energy E ω in color representation R can be compactly cast as [67] ω dI dω = 2ᾱπ Hereᾱ = α s C R /π and K(x, t 2 ; y, t 1 ) is an effective emission kernel describing the broadening of the emitted gluon during its formation. It corresponds to the evolution operator of a quantum particle immersed in the imaginary potential iv(x) in 2+1 dimensions and obeys the Schrödinger equation which resums multiple scatterings of the radiated gluon with the medium between the emission times t 1 and t 2 in the amplitude and its complex conjugate, respectively. For a general potential v(x, t) that includes the Coulomb tail at large momentum transfers, a closed form solution to Eq.  nevertheless be obtained for two special choices of the potential: vacuum and harmonic oscillator. In the vacuum case, setting v(x, t) = 0 leads to the following solution of Eq. (2.24) where ∆x = x − y and ∆t = t 2 − t 1 . Note that this contribution is explicitly removed in Eq. (2.23) so that the result is only sensitive to the purely medium-induced contribution.
In fact, we can also express the resummed propagator, given by the solution of Eq. (2.24), as a Dyson-like iterative equation that resums multiple interactions around the vacuum solution, namely (2.26) From the structure of the equation, we immediately see that, for a time independent rate v(x, t) = v(x), the function K only depends on τ ≡ t 2 −t 1 . This equation is equivalent to an expansion in medium opacity χ, defined as χ = L/ mfp . Computing the radiative spectrum -10 -by truncating the expansion in Eq. (2.26) at a fixed order in v(x, t), or χ, corresponds to the Opacity Expansion introduced in the previous section. Consistently, the n−th term in the OE scales as dI n /dω ∼ O χ n . The single scattering solution corresponds to the n = 1 truncation of the expansion and is often referred to as the GLV spectrum [22,26]. The other special case where Eq. (2.24) is analytically solvable is when v(x, t) = v HO (x, t) = 1 4q (t)x 2 ⊥ , that is, when the potential reduces to that of an harmonic oscillator. We recall that the IOE splits the leading logarithmic potential given in Eq. (2.15) as where Q is for now an undetermined matching scale, different from the one used for the broadening case. Yet, the effective jet quenching parameter isq(t) =q 0 (t) log Q 2 /µ 2 * . Thus, similarly to transverse momentum broadening discussed in the previous section, the solution to Eq. (2.24) with a quadratic potential, that we denote as K = K LO , corresponds to the leading order (LO) term in the Improved Opacity Expansion. It reads (2.28) Here, C(t 2 , t 1 ) and S(t 2 , t 1 ) are purely time dependent functions which are solutions to the initial condition problems [27] with the complex harmonic oscillator frequency Ω(t) given by More details on the properties of these functions can be found in Appendix A. Inserting Eq. (2.28) back into Eq. (2.23) and performing the time integrals, one obtains the spectrum at leading order in the IOE (or equivalently in the harmonic approximation). The final expression reads, and is often referred to as the BDMPS-Z spectrum [13,17]. The LO contribution to the IOE spectrum takes a particularly simple form in the case where the medium has an extension L with a constant density n; we refer to this simple medium model as the plasma brick model. In the brick model one can simply define the jet quenching parameter aŝ q(t) =q Θ(L − t), which allows one to write the C and S functions as S(t 2 , t 1 ) = 1 Ω sin Ω(t 2 − t 1 ) , and C(t 2 , t 1 ) = cos Ω(t 2 − t 1 ) . (2.32) In this case, the well-known behavior of the spectrum at asymptotically at low and high frequencies is for ω ω c , where the characteristic gluon energy ω c =qL 2 /2 corresponds to gluons with maximal formation time, i.e. t f = L. The behaviour in the soft limit highlights the Landau-Pomeranchuk-Migdal (LPM) interference [14,15] that occurs since the gluon is formed over timescales involving multiple interactions with the medium. The strong suppression at high gluon energies follows directly from the approximation of multiple soft interactions, implicit in the harmonic form. At these frequencies, i.e. ω > ω c , the contribution from a single, hard scattering can be shown to dominate, as we will discuss below.
Let us now construct the contributions to the IOE beyond the LO term. Adopting the decomposition provided by Eq. (2.27) that allows us to separate the harmonic part from the x dependent Coulomb logarithm, and in analogy to the resummation around the vacuum solution given by Eq. (2.26), the full kernel can be written as Truncating this relation at O(δv 2 ), it is easily seen that the LO kernel is given by K LO in Eq. (2.28). The NLO kernel reads which can be used in Eq. (2.23) to compute the NLO contribution to the IOE spectrum, as was done for the LO term. Like in the broadening case, we do not consider higher order terms since truncating the series at NLO is enough to reproduce the single hard and multiple soft regimes. At this order, the spectrum reads [51][52][53] and we defined the ratio Cot(t 2 , t 1 ) ≡ C(t 1 , t 2 )/S(t 2 , t 1 ). We now analyze the asymptotic forms of the spectrum by considering the brick model. In this case Eq. (2.37) reduces to At high frequencies, that is ω ω c or ΩL 1, one finds that Eq. (2.38) leads to k 2 (s) iω/(2s). The high-frequency behavior of the NLO term is given by It dominates the spectrum, given the quadratic ω suppression of the LO term, see Eq. (2.31). In this last equation, we recall that the medium opacity parameter is χ ≡q 0 L/µ 2 * ∼ L/ mfp and we introduced the high energy cut frequencyω c = 1 2 µ 2 * L. Higher-order terms are all suppressed by at least one additional power of 1/ω as well [53]. Thus, similar to the discussion for the broadening distribution P(k), one observes that the dominant term, given in Eq. (2.39), comes solely from the NLO contribution and it can be shown to exactly match the medium-induced spectrum obtained by considering a single hard scattering in the medium [26,51], i.e. n = 1 in the traditional Opacity Expansion. Furthermore, Eq. (2.39) is independent of the matching scale, analogous to what was observed for P(k) in Eq. (2.21).
At low frequencies, i.e. for ω ω c or ΩL 1, the NLO term, containing the single hard scattering physics, can be simplified by noticing that k 2 (s) −ωΩ, leading to [51][52][53] which is equivalent to the next-to-leading logarithmic result derived in Ref. [49]. Again, as observed for the broadening distribution, in the soft regime higher order terms in the IOE scale as the LO contribution, see Eq. (2.33), with increasing power suppression by [53]. This result implies that, in the soft limit, the full spectrum can be written in terms of the LO result with an effective jet quenching coefficientq eff that absorbs the additional logarithmic dependencies.
More importantly, Eq. (2.41) imposes further constraints on the matching scale Q 2 . To see this, let us first assume that the matching scale associated to the radiation spectrum Q 2 = Q 2 r is independent of ω. Then, in Eq. (2.41) the logarithms in the denominators would be frozen. However, the numerator logarithms would evolve quite rapidly for µ 2 * ω 2 (q 2 Q 4 r ) −1 , leading to a divergent series (notice that the LO contribution would be negligible in this case). Thus, one concludes that Q r = Q r (ω) in order for the spectrum to be free of unphysical divergences. In addition, one sees that the natural way to regulate the numerators is to take 8 Moreover, it can be shown [53] that this form follows directly from the fact that once all orders in the IOE are resummed the spectrum takes the functional form of the LO term. For the present paper and the following calculations, the main message is that Eq. (2.42) ensures that at low energies the spectrum is well behaved and non-physical divergences are absent. Also, and again in analogy to the broadening, at leading logarithmic order Q 2 r ∼ √q 0 ω, which using the above relations for the gluon formation time and the average accumulated momentum, can be translated into the typical momentum acquired by a gluon with frequency ω ω c . The solutions of Eq. (2.42) are discussed in Appendix B. Finally, we still need to ensure that Q 2 r µ 2 * in order to justify the expansion. Ignoring the logarithmic dependence in the matching scale, we observe that the IOE approach only works if where we defined the characteristic Bethe-Heitler (BH) frequency as ω BH = µ 4 * /q 0 . This condition means that the current scheme is not valid in the BH regime [68], see Ref. [34] for a similar conclusion and further discussion regarding the analytic treatment of the BH region. This regime is characterized by gluons with a formation time of the order of the mean free path in the medium, acquiring a momentum k 2 ⊥ ∼q 0 fmp ∼ µ 2 * and with a typical energy ω BH ∼ T of the order of the medium temperature. At this scale, nonlinear dissipation effects take place [69] such as gluon absorption. However, in the case of large or dense enough media (such that Q 2 m 2 D ) the BH regime is power suppressed and radiative energy loss is dominated by frequencies in the deep LPM regime in the calculation of inclusive jet observables [70].
In Fig. 2, we compute the medium-induced single gluon spectrum up to NLO in the IOE, comparing with a full numerical solution to Eq. (2.23) [33] and the GLV spectrum, corresponding to the limit of single scattering in the medium. The gray band indicates the region in which Eq. (2.42) does not have a valid solution, i.e. where the IOE approach is not valid. A similar numerical comparison was previously carried out in Ref. [34]. As discussed above, we numerically observe that in the soft sector, ω BH ω ω c , the difference between the full result and the LO contribution is small, and including the NLO provides a very good approximation. In addition, the IOE has no divergences since the matching scale is chosen for each ω by solving Eq. (2.42). At frequencies ω ω c , we observe that the LO is power suppressed, but the NLO term matches the full result, even faster than the GLV approximation. Overall, the agreement between the LO+NLO result and the full numerical solution is outstanding.

The medium-induced radiative kernel with the IOE
After having revised the building blocks of the IOE, we proceed to compute the fully differential medium-induced spectrum for a gluon with energy ω and transverse momentum k. We assume that the emitted gluon is soft, ω E, and collinear, θ 2 ∼ k 2 ⊥ /ω 2 1, with E being the energy of the emitter. The emitter follows an eikonal trajectory and its kinematics are frozen. Regarding the medium properties, which are encapsulated by the jet quenching parameterq, we assume that it has a smooth time profile almost everywhere and that, at large distances, the system reaches the vacuum sufficiently fast, i.e. lim t→∞q (t) = 0.
Then, we study a particular scenario where the medium has a simple time dependence: up to a distance L the jet quenching parameter is positive and constant, while for times larger than L,q = 0. This corresponds to the previously mentioned plasma brick model, where the medium is a slab with longitudinal size L, after which there is vacuum; mathematically it corresponds to defining the jet quenching parameter asq(t) =q Θ(L − t).
Under these assumptions, the purely medium-induced spectrum can be expressed as a convolution between the broadening probability distribution, P, and the splitting kernel, K, that we have introduced in the previous section. It reads, where t 1 and t 2 correspond to the gluon splitting light-cone times in amplitude and conjugate amplitude respectively, and span from the creation point inside the medium at t 1 = t 2 = 0 up to any possible in-vacuum or in-medium splitting time. In Eq. (3.1), we explicitly denote the starting time t 2 in the broadening distribution so, compared to our formulas in Section 2.1, P(x, L) ≡ P(x, L; 0). Also, in Eq. (3.1) we employ the adiabatic turn-off prescription [22], which prevents the emission of purely vacuum-like radiation at asymptotically large times, with the → 0 limit being implicit for the rest of the paper. The last term in the formula subtracts a contribution corresponding to purely vacuum -15 -radiation off the hard emitter given by (see Appendix C) Before proceeding further with the explicit analytic evaluation of Eq. (3.1), we anticipate a subtlety when carrying out the time integrals with the adiabatic turn-off prescription. Ignoring the e − t 1 suppression factor, the t 1 integral can be performed using Eq. (3.17). Keeping the prescription yields only an additional negative vacuum-like term, −(2π) 2 ω dI vac dωd 2 k , so that Eq. (3.1) can be expressed in a more convenient form as follows where now the for the t 1 integral prescription has been removed at the cost of a factor 2 multiplying the vacuum term. The limit → 0 has to be taken after the integral over t 2 . The details regarding the treatment of the adiabatic prescription are discussed in Appendix C.
In what follows, we will compute Eq. (3.1) in the IOE approach, including all terms up to O(δv) (NLO). For that, we extend Eq. (2.6) for a generic medium and express the broadening distribution P as Similarly, the emission kernel K can be expanded as in Eq. (2.34): Truncating these relations up to NLO accuracy and inserting them into Eq. (3.1), we obtain the spectrum which we write in the following way, To reiterate, the LO and NLO terms resum arbitrary number of soft medium interactions, encoded in v HO , and a fixed number (zero for LO, and one for NLO) number of hard interactions with the medium, through the potential δv.
While the vacuum spectrum is already given in Eq (3.2), the medium read as follows, and The LO term captures the physics associated with the production of gluon radiation due to multiple soft scattering in the medium, thus recovering the BDMPS-Z solution. The first term in the NLO term Eq. (3.8) includes the possibility of producing the gluon due to a hard scattering in the medium and when integrated over k gives the NLO contribution to the integrated spectrum studied in the previous section, Eq. (2.36) [51]. Finally, the last term in Eq. (3.8) arises from expanding the final state broadening distribution P.
Thus, it only affects the redistribution of the radiated gluon transverse momentum and it vanishes upon integration over k. In the following sections, we proceed to explicitly compute Eqs. (3.7) and (3.8).

Leading order contribution
The leading order contribution to the spectrum is captured by Eq. (3.7). The broadening distribution P LO , implicitly given in Eq. (3.4), reads where we define the bare saturation scale as a slight generalization of Eq. (2.17), reading 12) and the matching scale, Q b , satisfies (following Eq. (2.18)) Furthermore, the kernel K LO can be found in Eq. (2.28) but, for consistency, we also repeat it here in a slightly different form, namely where we recall that (see Appendix A for further details) Also, in Eq. (3.14) we have taken advantage of the anti-symmetry of S, i.e. S(t 2 , t 1 ) = −S(t 1 , t 2 ).
In all these functions, the value ofq enters as an argument and thus they are sensitive to the definition of the matching scale for radiation, Q r , that, as discussed in Section 2, is obtained by solving the transcendental equation (see Eq. (2.42)) At this point, an important remark is in order. The functional form of the matching scale is constrained by making the spectrum finite in the infrared. This leads to the ω-dependence in Eq. (3.16), which is constrained in this way up to an overall numerical coefficient. As was shown in Ref. [53], the dependence in such a factor is sub-leading for a fixed order calculation in the IOE. As such, and since all the time dependence of Q 2 r emerges fromq 0 , it is more convenient to define Q 2 r as a static scale. This simplifies the time integrations needed for the spectrum calculation without downgrading its accuracy.
In order to further simplify Eq. (3.7), we make use of a series of identities satisfied by the functions C(t 2 , t 1 ) and S(t 2 , t 1 ) that enter in the definition of K LO , see Eq. (2.29). In particular, we use Eq. (A.8) to obtain where in the second step we dropped an infinite phase which has already been accounted for in the vacuum subtraction term in Eq.
where the term proportional to ∂ x · x x 2 ∝ δ (2) (x) is purely imaginary and thus does not contribute given the prescription adopted in Eq. (3.3) (see Appendix C for more details).
The result obtained in Eq. (3.18), although compact, is somewhat obscure from a physical perspective. A more intuitive description of the in-medium emission can be achieved by using a momentum space representation In this form, we can interpret the first term in Eq. (3.19) as describing the emission of a gluon via some effective kernel at time t 2 followed by final state broadening. The second term corresponds to a vacuum-like subtraction contribution.
Furthermore, since P LO (p) is Gaussian, the remaining momentum integral can be performedˆp where we introduced the function 21) and the effective saturation scale is defined as Q 2 This form of the spectrum was first derived in Ref. [71]. Integrating over k and usinǵ k P(k) = 1, we recover the LO contribution to the energy spectrum discussed in the previous section [27,51], see Eq. (2.31).
As a sanity check, one can verify that Eq. (3.22) vanishes in the vacuum, i.e. in the limit Ω → 0, so that Thus, we obtain Plasma brick model. We proceed to evaluate the previous expressions for a concrete medium model, namely the brick whereq(t) =q Θ(L − t). Inside the medium, i.e. when both t 1 < L and t 2 < L, the C and S functions take simple forms (see Appendix A) and . On the other hand, for both t 1 > L and t 2 > L, the system evolves as in vacuum and the C and S are obtained by setting Ω → 0 in the previous equations such that This sharp separation of the problem into processes happening inside the medium and outside of it, suggests that an efficient way to evaluate Eq. (3.22) consists in splitting the time integral into two regions: one where 0 < t < L and a vacuum-like region where t > L. We refer to the first region as in-in since the gluon emission occurs inside the medium in both the amplitude and its conjugate. It can be easily obtained by replacing the upper -19 -limit in the integral in Eq. (3.22) by L and employing Eq. (3.25). Since the integral over t 2 has a finite extension, we can safely neglect the adiabatic suppression factor. This contribution to the total medium-induced LO spectrum reads where we have used that the saturation scale reduces to Q 2 s (∞, t) =q(L − t). The remaining region of phase space is obtained by imposing t > L in Eq. (3.22). In this case, there are two types of contributions: (i) a purely vacuum term corresponding to the scenario where the gluon is outside the medium both in the amplitude and its conjugate, in this situation the first and second terms in Eq. (3.22) cancel by construction, and (ii) an interference term where the amplitude gluon is emitted inside the medium while its conjugate counterpart is emitted in the vacuum (or vice-versa). The latter contribution, which we shall refer to as in-out, requires further manipulations.
We begin by constructing the C and S functions which have support inside and outside the medium. This is done by using the decomposition for the C and S, given in Eq. (A.9) [27] In addition, in Eq. (3.25) the broadening term (encapsulated in Q 2 s ) has support in (t, L) which is the vacuum region. Then, there is no final state broadening and one can take or, equivalently, Q 2 s = 0 in Eq. (3.25). Combining all these results, we find that where we have included the vacuum subtraction term, and where the implicit adiabatic prescription ∼ e − t in the first line allowed to drop the contribution from t → ∞.
-20 -This contribution, together with Eq. (3.27), constitute the medium-induced leading order spectrum, analogous to the BDMPS-Z result. The medium-induced spectrum at LO in the IOE reads then,

Next-to-leading order contribution
The computation of the next-to-leading order contribution to the spectrum can be done using similar manipulations to the ones performed for the LO term in the previous section.
The NLO spectrum is defined in Eq. (3.8). The first term corresponds to a genuine correction to the emission kernel, referred below to as the in contribution, while the second term, which we shall refer to as broad contribution, introduces the possibility of a hard scattering in the final state broadening process. Also, the vacuum-like subtraction terms that appeared in the LO contribution are absent at O(δv).
To summarize, the two contributions to the NLO spectrum read where we have rearranged the time integrations. 9 Let us begin by considering Eq. (3.34). Note, that both the t 1 and s integrals have support only inside the medium, hence the naming as the in contribution. By using Eq. (3.17), we can directly perform the s-integral such that The remaining derivative operator gives so that the spectrum further simplifies to (adopting hereafter the more compact notation Cot 21 ≡ Cot(t 2 , t 1 ), C 12 ≡ C(t 1 , t 2 ) and S 12 ≡ S(t 2 , t 1 )) The remaining integration in x is Gaussian, and can be executed to obtain Let us re-emphasize that, in the previous expression, the matching scale associated with Q s inP 21 is Q b . Replacing δv by its explicit definition, that depends on Q r , the in contribution reads Notice that we tend to write the largest time as the first argument in the functions. Nonetheless, this is not always possible since in general the C function has no definite parity under the exchange of the arguments, unlike the S function which is always odd. Comparing Eq. (3.40) to the LO contribution given by Eq. (3.18), we observe an additional transverse integral in the z variable which is no longer Gaussian due to the logarithmic dependence in δv. Nevertheless, this integration can be performed analytically too. The angular part can be performed by recalling the definitions of the Bessel functions of the first kind, that lead to (3.43) The more challenging radial integral can be solved by using a convenient decomposition of the logarithmic function. Namely, the relation allows us to transform the original integral into a sum of Gaussian integrations that can be readily performed. In particular, the z-integrals in Eq. (3.42) can be compactly expressed as and Then, replacing the logarithms in the previous equations by the decomposition in Eq. (3.44) allows one to write I a and I b in terms of the exponential integral function Ei, leading to and Taking advantage of these further simplifications, the in contribution to the NLO gluon spectrum can be compactly written as Hence, we have managed to reduce the number of integrals over transverse positions and times down to two time-integrations. Turning now to the broad contribution, given in Eq. (3.35), we can perform the derivatives on K LO and integrate over t 1 using Eq. (3.17). Then it reads whereP 2 20 ≡P 2 (t 2 , 0) given by Eq. (3.21).
Plasma brick model. So far, we have made no approximation regarding the time profile of the medium. As for the LO contribution, we now assume the plasma brick model, i.e. q(t) =q Θ(L − t). As in the previous case, this simple model allows one to further simplify Eqs. (3.49) and (3.51) by splitting the time integrations appropriately. We start with the broad term since it has only one time integral left. In addition, the spectrum is proportional to Q 2 s0 , and thus this term only has support inside the medium t 2 < L. As a consequence, the C and S functions are given directly by Eq. (3.25), and Eq. (3.51) reduces to Next, let us analyze the in contribution given by Eq. (3.49). It has support both inside and outside of the medium and thus two contributions appear. One option is that the gluon is emitted inside the medium in the amplitude and its conjugate, which we identify as the in-in contribution. The second term refers to the case in which one of the emissions happens outside of the medium, that we denote as in-out. The in-in term obeys the time ordering´L 0 dt 2´t 2 0 dt 1 in Eq. (3.49), with the C and S having support only inside the medium and thus are given by Eq. (3.25). Therefore, this contribution reads , and the auxiliary functions reduce tô Ωcot(Ω(t 2 − t 1 )) + Ωcot(Ωt 1 ) + 2iω The in-out contribution can be further simplified following similar steps as in the LO case. Now, the time integrals read´∞ L dt 2´L 0 dt 1 and one can set Q 2 s = 0 everywhere in Eq. (3.49) since this term only has support outside of the medium. Then, the spectrum reads where we recall that the C and S functions are distinct from the ones used in the in-in term, since they have support both inside and outside of the medium. Nonetheless, as for the LO case, they can be written in terms of the purely in-medium and in-vacuum C and S functions by using Eq. (A.9). Taking t 0 = L, we find that for the above time ordering ) Ω , such that . (3.58) For the reversed time ordering, we find leading to Combining all these results, the auxiliary functions now read Inserting these expressions into Eq. (3.56), one realizes that the remaining t 2 integral can be carried out As a consequence, the in-out contribution to the NLO spectrum can be finally written as 2ωΩ tan(Ω(L−t 1 )) .

Final formulas
At this point, we summarize the main results obtained in the two previous sections. Our aim is to provide a set of compact equations which can be directly used in phenomenological studies or implemented in jet quenching Monte-Carlo codes. In what follows, we first present the results for a generic medium profile and then take the brick limit.

Spectrum at LO+NLO for a generic medium profile
Up to NLO in the IOE, the purely medium-induced gluon spectrum, i.e. after subtracting vacuum radiation, can be written as The LO term reads In addition, the NLO contributions are given by (3.69) In the above equations we introduced . The remaining functions are defined as follows, The matching scale Q b , that enters everywhere into the broad and only into the Q s definition for the in case, is obtained by solving the transcendental equation withL some effective medium length. The exact value ofL is not important, as long as it is taken such that the relevant support for the integration ofq is covered. The radiative matching scale, Q r , appears in all other terms related to the kernel expansion and is the solution of (3.73)

Spectrum at LO+NLO for the brick model
The previous results are simplified when the medium is modelled as a plasma brick of length L withq (t) =q Θ(L − t) . (3.74) In this case, the full medium-induced spectrum at NLO in the IOE can be written as The leading order terms read Here we used The NLO in-in term can be written as where Q 2 s (L, t 2 ) is defined as above, the C and S functions are given in Eq. (3.25), and Ωcot(Ω(t 2 − t 1 )) + Ωcot(Ωt 1 ) + 2iω (3.81) The NLO in-out piece is Finally, the NLO broad contribution is given by

Asymptotic behavior
The complete expressions for the in-medium branching kernel, that we have summarized in the previous section, are written in terms of a few integrations that we did not manage to solve analytically. Before presenting their numerical implementation, we would like to give further analytical insight into the discussion. To that end, we analyze the behavior of the IOE spectrum up to NLO in two physically relevant asymptotic regimes: when the emitted gluon is either (i) soft (ω ω c ) and collinear (k 2 ⊥ qL) or (ii) hard (ω ω c ) and wide angled (k 2 ⊥ qL). That is, the regime of validity of BDMPS-Z and GLV approaches, respectively. Our results below are obtained by taking the brick limit, although similar conclusions are obtained for other choices of medium profile. In addition, we neglect the purely vacuum radiation as it is completely irrelevant for this discussion.

Multiple soft scattering regime
We begin by analyzing the regime in which the emitted gluon is soft, i.e. ω ω c , and its typical formation time is much shorter than the medium length t f ∼ ω q L. The latter condition can be translated into a constraint on the transverse momentum off the emission, q 2 ⊥ ∼qt f ∼ √q ω qL, 10 which implies that short formation time gluons typically acquire most of their transverse momentum due to final state broadening. Under these conditions, the IOE spectrum simplifies significantly. The general formula for the medium-induced spectrum given by Eq. (3.1) can be re-written in momentum space as 11 where we have exploited that K and P are invariant under time translations, i.e. depend only on time differences, when the plasma is homogenous. Eq. (3.86) can be simplified noting that τ ∼ t f t 2 and thus one can set t 2 → ∞ in the τ integration upper limit. That is, the two time integrations decouple. In addition, we note that q ⊥ ∼ 1/x ⊥ corresponds to the transverse momentum acquired in the branching process, q ⊥ ∼ √q ω, that, as we have anticipated, is small with respect to the characteristic broadening momentum, i.e. q ⊥ qL. As a consequence, we neglect q with respect to k inside P. Then, the q integral acts solely on K and the x and q integrals yield

(3.87)
A familiar element in the previous equation is the medium induced rate defined as Then, Eq. (3.87) can be finally written as That is, in the soft and collinear limit, the spectrum is given by the product of the time integral of the broadening distribution and the medium induced energy rate. This result is not tied with the IOE approach and has been previously obtained in the literature in the context of BDMPS-Z calculations [19] and exploited in Monte-Carlo simulations [55,57]. We proceed to compute Eq.(3.89) in the IOE. 10 Notice that this condition refers to the momentum of the in-medium vertex, rather than the final momentum of the gluon. Even for soft gluon emissions, final state broadening can lead to a final momentum k 2 ⊥ ∼qL. 11 Strictly speaking, the upper bound of the integrals should be L. We take L → ∞ to facilitate analytical manipulations.
-29 - The soft limit of the IOE energy spectrum at all orders was computed in Refs. [51,53] and reads where ω dI LO dω =ᾱ qL 2 ω corresponds to the well known BDMPS-Z result. The effective jet quenching parameter is given at leading-logarithmic order by [53] 12 Eqs. (3.90) and (3.91) show that the IOE energy spectrum is governed by the LO result with higher orders suppressed by a logarithmic power which can be written in terms of the ratioq 0 q . A similar conclusion can be reached regarding the broadening distribution in the kinematical limit k 2 ⊥ qL, see for example Eq. (2.22) for the result up to NLO. Combining these two results, one concludes that the soft and collinear limit of the fully differential spectrum in Eq. (3.89) also obeys this functional form. Note that the spectrum will consist of terms where the matching scale is given by Q r and others where it is Q b , depending if the terms come from the expansion of the kernel or of the broadening distribution.
In Appendix D we explicitly show that the in-out NLO contribution in the IOE spectrum scales as the LO term multiplied by a logarithm that arises from the ratioq 0 q .

Rare, hard scattering regime
Let us now consider the orthogonal regime with respect to the previous section. Here, the gluon is hard, ω ω c , and carries a large transverse momentum k 2 ⊥ qL, i.e. the intrinsic momentum of the gluon is significantly larger than what it typically acquires through broadening in the medium. In this case, the multiple soft scattering contribution is suppressed by the LPM effect and the emission spectrum is dominated by single hard scattering in the medium. This corresponds to the truncation of the opacity expansion, considered by GLV, at first order [22,26], leading to an emission spectrum reading where u = L 2ω k 2 ⊥ and γ = µ 2 L 2ω . The medium potential was taken to be the GW model and thus µ is its infrared regulator. It can be related to the universal physical mass µ * , as mentioned in Section 2.1 and detailed in Refs. [50,53].
Note that the conditions ω ω c and k 2 ⊥ qL correspond to u 1 γ in the previous equation. To take this limit in Eq. (3.92), let us consider the integral with u held constant. The contribution from the first region can be easily computed to leading-logarithmic accuracy where we introduced 13 (3.96) In the case of I > , we first notice that x 1 so we can drop the sin(x) term. Defining a ≡ γ/u 1 we have

(3.97)
Combining the results from the two different regions, the I integral gives Inserting this result into the GLV spectrum given in Eq. (3.92) yields The expected 1/k 4 ⊥ power tail naturally arises from a Coulomb-like single hard scattering in the medium. Counterintuitively, even though we are considering here the high energy limit, the spectrum is still sensitive to the infrared details of the in-medium scattering potential via the thermal mass µ. For the sake of comparing with the IOE in what follows, a couple of manipulations are required. First, we re-write the resulting logarithm as and then replace the GW mass by the universal infrared scale µ * through the leadinglogarithmic prescription 4µ 2 * = µ 2 e −1+2γ E [50,53]. This leads to our final expression for the GLV spectrum: Adding the two components results into a vanishing spectrum at this order in k and indicates the need to go to higher orders that, as we will see, affect rather differently the in-in and in-out terms. The latter is exponentially suppressed when including higher orders as can be derived from Eq. (3.32). In turn, the in-in term follows a power-law suppression. To see this, we perform a second order gradient expansion of the broadening distribution in the Q 2 Now, because k is large the phase oscillates rapidly unless t 2 is small enough. To estimate the support of the t 2 integral, we exploit the fact that in the high energy regime ΩL 1 and thus the dominant contribution to the integral comes from the region where t 2 2ω k 2 ⊥ Ω −1 . As a consequence, one can replace the integration limit L → ∞ and linearize the tangent, to obtain the leading asymptotic behavior of the spectrum Notice that the logarithm depends on the broadening matching scale since it originates from the last line in Eq. (3.104). Interestingly, the leading order contribution exhibits a 1/k 4 ⊥ tail, physically corresponding to early time hard emissions, which then suffer multiple scatterings in the medium, acquiring a momentumqL much smaller than the momentum off the emission vertex. Compared to the vacuum like emission in Eq. (3.102), we observe that although final state broadening does not change the power-law dependence on the transverse momentum, it power suppresses this second order by aq L k 2 ⊥ 1 factor. At NLO, we need to analyze individually each of the three identified terms in this asymptotic regime. Starting with the broad term, we note that in the high energy limit It is convenient to write the remaining integral aŝ (3.110) Truncating the previous exact result to leading order in x max , we obtain Plugging this last result into Eq. (3.108) yields We restrain ourselves from doing any physical interpretation of this result at this stage and proceed to compute the in-in contribution. In this case, the lack of a vacuum-like (dt/t) divergence simplifies the calculation. First, we use the asymptotic form of the exponential integral function, in Eqs. (3.47) and (3.48), to obtain the leading asymptotic forms of I a and I b Moreover, since ΩL 1 we can simplify the auxiliary functions given by Eq. (3.55) down toR .
Consequently, the in-in spectrum reduces to Finally, for the in-out term we obtain which is power-suppressed with respect to the broad and in-in terms and can be ignored. Then, the NLO contribution to the IOE spectrum is given by Finally, combining the LO and NLO results we obtain that the IOE spectrum at high energies reduces to A few important remarks are in order at this point. The most obvious one is that the LO+NLO exactly matches the GLV result in the large ω, large k regime. This was somehow expected but not trivial to confirm explicitly. Among other technicalities, a second order gradient expansion in transverse momentum for the LO term was essential. Another remarkable and related feature of the final result is that both the LO and NLO depend on Q b , while their sum does not. This was already encountered in the energy spectrum calculation as discussed in Ref. [53] and constitutes a sanity check of the Improved Opacity Expansion framework. We note again that unlike the soft limit considered before, this regime, although having a non-trivial cancellation of the matching scale dependence between different orders, does not provide any constraint on the functional form of Q b , as can be observed from Eq. (3.119). This is unlike, for example, the result in Eq. (D.7), that forbids the matching from being a numerical constant. From a more pragmatic point of view, the exact matching between the IOE and GLV provides a non-trivial check on the computations performed in the previous sections.

Numerical results
In this section, we numerically explore the IOE spectrum in the brick model. We compare our results of the IOE spectrum truncated at LO and LO+NLO (see Section 3.3.2) to (i) the single, hard scattering limit encompassed in the GLV spectrum (see Eq. (3.92)) and (ii) an all-orders resummation of the spectrum presented in [33]. Notice that the LO result can be considered as the BDMPS-Z solution with ultraviolet regulator taken to be the radiative matching scale Q r . These comparisons should be regarded, at this point, as a merely theoretical exercise. However, we choose the medium parameters to be in the ballpark of LHC conditions. More concretely, the LHC-inspired medium has:q 0 = 0.156 GeV 3 , length L = 6 fm and infrared regulator µ * = 0.355 GeV. Further, we take a fixed value of the strong coupling constant α s = 0.28 and consider radiation off a hard quark such that C R = C F = 4/3. These set of parameters lead to a critical frequency scale ω c0 ≡q 0 L 2 = 140 GeV and the saturation density Before comparing our result to other approaches, we first address a natural question: what is the dependence of the IOE spectrum on the matching scales Q r and Q b ? In Fig. 3 we plot the medium-induced spectrum as a function of k ⊥ /Q s0 for two different gluon frequencies, ω = 5, 100 GeV, truncating the spectrum at LO (left) or at LO+NLO (right). The central curves are obtained by solving Eqs. (2.18) and (2.42). Then, we perform an independent variation of the matching scales by a factor of 2 (1/2) that lead to the uncertainty bands around the mid value. We recall that this variation is associated to the uncertainty in the definition of such scales, which are constrained up to an overall constant factor.
Let us discuss first the large ω, large k ⊥ regime, i.e. the inset in the bottom row plots. Analytically, we have shown that, in the asymptotic kinematical region, the dependence on the matching scale vanishes when one considers LO+NLO contribution, see Eq. (3.119). This is exactly what we observe numerically. Although this conclusion was reached for highly energetic gluons, the numerical results indicate that it holds reasonably well in the case of soft gluons too. Notice that when analyzing only the LO, a bigger, but still weak dependence on the matching scale Q b is observed. This corresponds to Eq. (3.106), where a logarithmic dependence on Q b appears. Then, we reemphasize that only when considering -35 -LO+NLO the dependence on the matching scale is residual as due to the cancellation occurring between these two terms.
The small ω, small k scenario is represented by the top row plots in Fig. 3. In this region, we argued in the previous section that all orders scale as the LO term, with logarithmic power corrections as one goes higher in the IOE, see Eqs. (3.89) and (3.90). Numerically, we observe that there is a large uncertainty due to the variation of the matching scales at LO, mainly from Q r . However, if one includes the NLO term (top right) the dependence on the matching scales almost disappears. Regarding Q b , we note that higher orders in the broadening factor P also enter through logarithmic power corrections on top of the LO term. Thus, the weaker dependence of LO+NLO on Q b as compared to LO is to be expected. These findings are inline with [53], where it was observed that for the energy spectrum the variations of the matching scale Q r although could drastically change the LO and the NLO terms separately, the sum LO+NLO was only sensitive to these variations at NNLO. This is a consequence of the fact once all orders in Eq. (3.91) are considered, the spectrum becomes independent of Q r . Since in Fig. 3 the largest uncertainty comes from the scale Q r , we argue that the result obtained here is a manifestation of the findings of [53].     Fig. 4 we present the final comparison between the IOE spectrum and the other approaches mentioned above. We considered two gluon frequencies, ω = 0.05 ω c0 and ω = 2 ω c0 , corresponding to the cases of soft and hard gluon and we use the GW mass µ 2 = 0.43 GeV 2 . Again, we vary the matching scales of the IOE spectrum up and down by a factor of two, independently and then take the envelope to build the uncertainty band.
The overall conclusion from the two plots in Fig. 4, and the most important result of this paper, is that the IOE spectrum, up to NLO, already does a reasonable job at capturing the full solution (less than 25% deviations), with the advantage that it requires considerably less computational power. The observed deviation from the full numerical result reflects the sensitivity of the transverse momentum distribution to the infrared. A wider separation between µ 2 and Q 2 b or Q 2 r would yield a better agreement. Let us split the discussion of Fig. 4 into small and large transverse momentum.
The small k ⊥ and small ω regime is dominated by multiple scattering contributions. Then, it is natural that the GLV spectrum fails to capture the full solution. The LO term of the IOE (related to the BDMPS-Z solution) already does a good job at describing the full result. Nonetheless, including the NLO term, improves not only the overall agreement with the full result, but also reduces the uncertainty band associated with the variation of the matching scales, as discussed in the previous section. When increasing the gluon's frequency, and still at small k ⊥ , we observe that the LO+NLO result remarkably captures the full solution up to 5% deviations. Regarding GLV, its agreement with the full solution is improved with respect to the small frequency case and, curiously, coincides with the LO term. We do not expect this to be a systematic result for other choices of the medium parameters.
The large k ⊥ tail is generated by rare, hard scatterings in the medium. It is well known that the BDMPS-Z approximation does not correctly capture such contributions and, therefore, fails to describe the full solution. At small frequencies, we observe that the LO+NLO result approaches much faster the full solution than compared to GLV. This is to be expected, as at large, but not infinite k ⊥ , multiple scatterings still play a role, despite being sub-leading. GLV lacks those effects and then needs an asymptotically large value of k ⊥ to reproduce the full solution, while the LO+NLO works even far from the asymptotic regime. At large frequencies, the spectrum is really dominated by a single hard scattering in the medium and thus the LO+NLO, full and GLV results rapidly converge.

Conclusion and Outlook
This work constitutes a natural extension of the recent studies of the medium-induced energy spectrum and broadening distribution using the Improved Opacity Expansion [50][51][52][53]. We have computed the in-medium radiative kernel using the IOE up to next-to-leading order accuracy, in the soft gluon approximation for a generic medium profile as well as for a brick plasma for which we have performed numerical computations that we compare to full numerical results from [33]. We observe a very good agreement for a LHC-motivated choice of medium parameters. , the IOE at LO+NLO (solid, red) and the all-order spectrum (solid, navy) as computed in Ref. [33] for two gluon frequencies: ω = 0.05 ω c0 (left) and ω = 2 ω c0 (right). The ratio to the full solution is presented in the bottom panels. The uncertainty band arises from variations in the matching scale.
From a theoretical viewpoint, the differential spectrum calculation highlights the role played by the matching scales that enters the definition ofq, given that it convolutes contributions due to final state broadening with in-medium radiative terms. Each of these physical processes enters the IOE expansion with its own matching scale that we denote by Q b , associated to final state broadening terms, and Q r , related to the radiative kernel. These two scales are obtained by solving their corresponding transcendental equations that are given, in full generality, by Eqs. (3.72) and (3.73). We emphasize that these two scales have to be treated separately in order for the expansion to be consistent and well defined. Taking this into account, we derive the medium-induced spectrum for a smooth medium time profile and also in the brick limit. The final formulas are given in Sec. 3.3.
Besides the master formulas, we analytically study the IOE for the plasma brick model in two asymptotic kinematical regimes. Firstly, we consider the regime where multiple soft exchanges with the medium constitute the dominant contribution, i.e. ω ω c , and with the further assumption that the gluon is collinear, i.e. k 2 ⊥ qL. In this case, we recover the well-known factorization formula given by Eq. (3.89), often used in jet quenching phenomenology [19,[55][56][57]. In particular, this result implies that in the soft limit the differential spectrum can be written as the product between the LO term and powers ofq 0 q that correspond to higher order contributions. This result agrees with what was observed in the energy spectrum calculation, as detailed in [53]. Secondly, we study the physically opposite regime where a single hard scattering with the medium governs the dynamics of the medium-induced spectrum. This corresponds to a region of phase space where ω ω c and final momentum of the gluon satisfies k 2 ⊥ qL. In this regime, it is expected that the exact spectrum is given by the GLV result, and thus should be reproduced by the IOE approach. Indeed, we confirm that after considering both the LO and NLO contributions, non-trivial cancellations between these two orders occur such that one recovers the GLV spectrum. Not only the GLV result is reproduced, but also the dependence of the spectrum on the matching scales disappears order by order in 1/k 2 ⊥ . Again, these cancellations resemble the situation in the energy spectrum calculation [53]. Additionally, the explicit and detailed calculation carried out in order to check that the IOE recovers the GLV result provides a non-trivial cross-check on the main formulas derived in this paper.
Regarding the final numerical evaluation we find that, for the plasma brick model with LHC-inspired parameters, the computing time is in the ballpark of the LO/BDMPS-Z result [35]. More concretely, we have evaluated the code's performance and encountered that the small ω and small k regime requires more computational power due to the oscillating phases in the integrands. We provide ancillary Python files with the IOE spectrum together with the GLV expressions in Ref. [62]. The comparison with an all-orders resumed spectrum reveals a globally good agreement between the NLO spectrum from the Improved Opacity Expansion and the full solution, for this set of medium parameters. In particular, the agreement improves at high-frequencies, where the deviations between the two approaches are below 10%. This is a remarkable result given the relative simplicity of the approach presented in this paper as compared to larger computational cost needed to resum the spectrum to all orders. It is indicative of the power of the IOE approach to cap--39 -ture the correct dynamics at small and large frequencies simultaneously. A more thorough comparison with the full numerical result including a scan of the parameters space is left for future work. We expect the agreement to systematically improve for denser or larger plasma for which the scale separation between Q r or Q b and the infrared scale µ 2 is larger. Obviously, the IOE scheme is exact asymptotically.
Our results provide for the first time a unified analytic framework for the fully differential medium-induced radiation spectrum that accounts for both the GLV and BDMPS limits. We expect that adopting our radiative kernel in future phenomenological studies would substantially reduce model dependence of jet quenching observables as well as theoretical uncertainties on the extraction of medium transport properties such as the jet quenching parameter,q that is a function of the typical scale of the process. Two phenomenological applications have been already proposed in the literature that concern quenching effects on the jet spectrum [36,37]. A natural continuation of this work is to use the in-medium radiative kernel that we have derived in this paper to analytically compute observables where the gluon transverse momentum information is not integrated out, namely jet substructure observables. In particular, on-going measurements of the k ⊥ -distribution of the hardest splitting [72] would benefit from a theoretical calculation in which both multiple soft scatterings and hard momentum exchanges are correctly incorporated. This study would open up a new theoretical window onto extending the IOE framework to describe the energy loss of a quark-antiquark antenna. In parallel, we would like to implement the formulas that we have derived in this paper in a suitable Monte Carlo framework such as Ref. [55,57].

Note
While this manuscript was being produced an independent derivation of the differential spectrum (for a massive quark) using the IOE approach was presented in Ref. [73]. Although performing a numerical comparison is beyond the scope of this work we would like to point out a couple of differences. Firstly, in comparison with [73] we have presented results for a generic medium profile for which we were able to reduce further the number of integration variables. Secondly, in Ref. [73] a single matching scale was used for the radiative and broadening parts, i.e., Q b = Q r which we found leads to an incorrect description of the spectrum.
where we recall that the C and S functions satisfy + Ω 2 (t) S(t, t 0 ) = 0 , S(t 0 , t 0 ) = 0 , ∂ t S(t, t 0 ) t=t 0 = 1 , For any Ω, i.e. for a generic time profile of the medium, one can derive certain identities relating the C and S functions that were employed in the main text to simplify the emission spectrum. Firstly, these solutions are related by C(t 1 , t 2 ) = ∂ t 2 S(t 2 , t 1 ) and by the associated Wronskian (W ), which reads where we used the above the initial conditions. The condition W = 1 can be used to show that , which is used to derive Eq. (3.17). In addition, W = 1 implies that C and S are linearly independent solutions and thus any other solution to the above ordinary differential equation can be written as a linear combination of them. Using this fact, and for a time ordering t 2 > t 1 > t 0 , any solution in (t 2 , t 1 ) can be written as [27] S(t 2 , where it is easy to verify that these equations satisfy the above initial conditions and that S(t 2 , t 1 ) = −S(t 1 , t 2 ). This decomposition of the C and S is extensively used in the main text; here we give a simple application to derive another useful identity. Let us consider the brick model introduced in the main text, with a medium of extension L such thatq(t ≥ L) = 0, but still allowing the jet quenching parameter not to be constant in time inside the medium. In the vacuum, the solutions to the C and S functions are trivially found S(t 2 , t 1 ) = t 2 − t 1 , C(t 2 , t 1 ) = 1 , (A. 10) and indeed when introduced back in Eq. (A.5) give Eq. (A.4). Using Eq. (A.9) with t 2 = +∞, t 1 > L and t 0 = L, combined with the explicit forms for vacuum C and S solutions, one observes that terms proportional to S(t 2 , t 1 ) dominate, leading to the handy formula where the last equality holds if C is even in its arguments. Finally, in the case of the brick model with a time independentq inside the medium, one finds that S(t 2 , t 1 ) = sin(Ω(t 2 − t 1 )) Ω , C(t 2 , t 1 ) = cos(Ω(t 2 − t 1 )) . (A.12) Note that for the special cases of vacuum and static medium profiles, the C function is even in its arguments.

B Solutions to the IOE matching scale
In this appendix we pin down some of the basic properties of the real solutions to Eq. (2.42). 14 It is convenient to rewrite Eq. (2.42) as a fixed point equation. Introducing which implies that the roots are in general not unique: either there is no real solution, there is a single solution or at most two real solutions. For a real solution to this equation to exist one must require that which can be derived from the limiting case where there is a single solution [34]. In particular we find that so that in the regime of interest where α 1, one always has 2 real solutions (Q − and Q + ) to Eq. (2.42), which are numerically obtained depending on the initial condition used in solving the recursion relation. Since one typically uses as an initial condition Q 0 µ * , the solution obtained is always the largest root Q + , while as α → ∞, Q − → 1. These observations are illustrated in Fig. 5, where we numerically confirmed that in the limit α 1, there are always two roots (i.e. two real solutions for Q) with the smaller one asymptotically approaching 1. In the opposite case, when α 1 there are no real solutions, as expected in the unphysical region where µ 2 * √q 0 ω.

C Details on the adiabatic prescription and the vacuum contribution
In this appendix, we discuss the vacuum limit of the full spectrum and the importance of the adiabatic prescription, ensuring that interactions are properly turned off at asymptotically large times. We also give a justification for the slightly modified adiabatic prescription that allows us to carry out time integrations for a general medium profile while correctly accounting for the vacuum contribution.
14 As noted in Ref. [37], the transcendental equation can be written in terms of special functions, i.e.

C.1 Revisiting the pure vacuum spectrum
Before turning to the in-medium calculation, let us revisit the calculation of gluon emission in the vacuum, with the goal of understanding how to tame the emission of vacuum-like radiation at asymptotically large times. The matrix element which encodes the production of soft radiation with energy ω and transverse momentum k off a hard parton with lightcone energy p + , reads where we have used k − = k 2 ⊥ /2ω and anticipated that only the transverse component of the gluon will contribute at leading order. Squaring the amplitude and inserting the necessary phase space, spin and color factors, one recovers the vacuum spectrum introduced in Eq. (3.2) To make contact with the main text calculation, we use the typical trick to write the denominator in an integral representation Going back to Eq. (C.1), we can express the square amplitude in terms of a double integral Although this representation seems unnecessarily more complicated than the one in Eq. (C.2), it is of the form of the medium-induced spectrum in Eq. (3.1). In addition, we observe that the adiabatic prescription we introduced in the main text is related to the Feynman prescription for the pole structure of the propagators. The integral over t 1 is easy to carry out and gives At this point, one could just perform the remaining t 2 integral and recover Eq. (C.2). However, since we know the correct result and this example is sufficiently simple, let us swap the implicit limit → 0 with the integral in Eq. (C.2) first. Doing this we would obtain for the vacuum spectrum which is twice the result in Eq. (C.2), thus showing that the limit does not commute with the integral. As a consequence, one must always perform first the integrations and only take the limit → 0 at the final step. Performing first the t 2 integral in Eq. (C.5), we obtain (C.7) Taking the limit → 0 in the second term we obtain − i 2 ( The divergent term, ∼ −1 , is purely imaginary and does not contribute to the squared amplitude. Keeping only the real term and combining with the first term in Eq. (C.7), we obtain matching Eq. (C.2) after including the necessary overall factors, as expected.
In the next section, we will detail how to deal with the adiabatic prescription in the case of in-medium emission. The strategy followed will be the same as the one detailed here, i.e. only take the → 0 limit at the end. However, the time integrals are not as simple to perform as in the vacuum case, so slightly more elaborate techniques are necessary, which we follow to describe. dt 1ˆx e − (t 2 +t 1 ) e −ik·x P(x, ∞; t 2 ) × ∂ x · ∂ y K(x, t 2 ; y, t 1 ) y=0 , (C.10) where we have introduced the adiabatic regulator , and the limit → 0 has to be taken after the integrations have been performed [22,75]. We have shown in the previous section that the t 1 and t 2 time integrals can be performed for pure vacuum radiation albeit with a careful treatment of interactions at infinity. In the medium, it is in general not possible to perform analytically both time integrals. However, we may use the fact that the primitive of ∂ y K(x, t 2 ; y, t 1 ) with respect to t 1 is known and carry out the t 1 integral, see Eq. (3.17). Note that the e − t 1 suppression factor stands in the way of the immediate application of Eq. (3.17). However, this is can be dealt with in two ways which we shall present in what follows. Let us first briefly outline the discussion presented in Ref. [71]. Defining the function f (t 2 , t 1 ) = 1 k 2 ⊥ˆx e −ik·x P(x, ∞; t 2 )∂ x · ∂ y K(x, t 2 ; y, t 1 ) y=0 , (C. 11) it is straightforward to show that at late times t 2 > t 1 > T , where P(x, ∞; t 2 ) = 1 and K(x, y) ≈ K 0 (x, y), the above function reduces to with T an arbitrary time chosen to be much larger than the typical medium extension, i.e. T L. Using this fact that at late times f can be written as pure phase e iφ(t 2 −t 1 ) , one can show that [71]  . Applying this to the full emission spectrum above, the subtraction term becomes the vacuum spectrum (2π) 2 ω dI vac dω d 2 k = 4ᾱπ k 2 ⊥ , (C.14) leading to the final spectrum for the medium-induced spectrum considered in Eq. (3.3). An alternative strategy consists in absorbing the e − t 1 prescription in K while preserving its integrability by following the chain of operations: . This shows that instead of using the e − t 1 adiabatic phase, one can do a slight rotation of the frequency ω in the complex plane. The generalization of the above discussion to the full K is straightforward since the prescription is only relevant at asymptotically large t 2 and t 1 . Let us show how this prescription applies to Eq. (3.1). The integral over t 1 readŝ t 2 0 dt 1 e − t 1 ∂ y K LO (x, t 2 ; y = 0, t 1 |ω) →ˆt (C.17) When the differential operator acts on the exponential factor one recovers the medium induced contribution given for example in Eq. (3.18); in this appendix we are interested in the first term, which was overlooked before. It is easy to show that is a valid representation for the Dirac delta function. Using this result, the x integral in Eq. (C.17) becomes trivial. After some simple manipulations, we obtain that the net contribution to the spectrum is simply (C. 19) This contribution exactly matches the term subtracted in Eq. (3.3). Also notice that as a consequence, in the main text calculation, one can ignore the action of any differential operator on x x 2 , since its contribution has already been taken into account. Note that the prescription is only relevant for the purely vacuum term emerging from Eq. (C.17). In the term considered in the main text, one can safely set ω → ω and → , owing to the fact that there is only one remaining time integral t 2 for which the details of the regularization, such as multiplying by an arbitrary factor, do not matter. In addition, one can show that the adiabatic prescription can be ignored for higher order contributions in the IOE. A simple way to see this, at NLO, is to introduce Eq. (3.9) and Eq. (3.10) in Eq. (C.16), after adjusting the time integration limits and inserting extra K LO factors inherited from the full NLO spectrum; see Eq. (3.8). Also, one must recall that the NLO contributions are proportional to δv. The calculation follows as for the above example, where the terms coming from ∂ x x x 2 , will be proportional to δv(x)δ (2) (x) = 0, showing that the adiabatic regularization can be overlooked at higher orders in the IOE.
-47 -D The in-out contribution to the IOE spectrum in the multiple soft scattering regime In this appendix, we explicitly show that in the multiple soft regime considered in Section 3.4.1, the NLO contribution in the IOE scales as the LO term multiplied by a logarithm that arises from the ratioq 0 q . For the sake of the argument, we focus on the in-out contribution, which although giving being sub-leading, still exhibits this scaling and it is straightforward to compute. In this kinematic regime ΩL 1, and we can use (1 + i) , (D.6) 15 In the first integral one could further approximate the tangent inside the logarithm and obtain instead The logarithm inside the brackets is small [53] since at leading-logarithmic order Q 2 r = √q ω, as discussed in Section 2.2. As detailed above, indeed we explicitly find that the NLO result scales as the LO term times a logarithmic factor encapsulated byq 0 q .