Improved opacity expansion for medium-induced parton splitting

We present a new expansion scheme to compute the rate for parton splittings in dense and finite QCD media. In contrast to the standard opacity expansion, our expansion is performed around the harmonic oscillator whose characteristic frequency depends on the typical transverse momentum scale generated in the splitting. The first two orders account for the high frequency regime that is dominated by single hard scatterings together with the regime of multiple soft scatterings at low frequency. This work generalizes the findings of Ref. \cite{Mehtar-Tani:2019tvy} beyond the leading logarithmic approximation allowing to account also for the Bethe-Heitler regime and compare to the full numerical results from Ref. \cite{CaronHuot:2010bp}. We investigate the sensitivity of our results to varying the separation scale that defines the leading order. Finally, the application to Monte Carlo event generators is discussed.


Introduction
Measurements of significant modifications of hard probes observables, in particular jets, in heavy ion collisions as compared to proton-proton collisions at RHIC and LHC have firmly established the prominent role of final-state interactions in the dense nuclear medium created in heavy-ion collisions. For large systems, radiative processes are the main mechanism responsible for the observed quenching effects [3] (for recent reviews see [4,5]). Therefore a precise description of these processes is essential for a quantitative understanding of the mechanisms driving in-medium jet modification and probing non-equilibrium dynamics of the quark-gluon plasma (QGP).
For energetic particles propagating close to the light-cone through a QCD medium, the problem reduces to describing the (non-relativistic) dynamics in the transverse plane perpendicular to the direction of propagation. A formalism for dealing with multiple scattering in a QCD medium was developed by Baier-Dokshitzer-Mueller-Peigné-Schiff [6,7] and Zakharov [8,9], usually referred to as the BDMPS-Z formalism. 1 This resummation [11] can also be cast as an expansion in medium opacity [12][13][14][15]. 2 The regime of strong interactions can be approximated by diffusive broadening of the transverse momentum, governed by the diffusion coefficientq. In this regime, characterized by the formation time of the radiation, t f , being larger than the mean free path mfp , i.e. mfp t f L, interference effects between subsequent scattering centers have to be taking into account leading to the Landau-Pomeranchuk-Migdal (LPM) suppression. In this regime, the transport coefficientq is proportional to the Coulomb logarithm which must be regulated in the UV by the typical transverse momentum acquired by many soft scatterings, q 2 med , see e.g. [18]. This approximation ceases to be valid for relatively hard emissions for which the dominant effect comes from single scattering with the medium [19,20]. Apart from a full numerical solution of the propagator [2,21,22], it is currently unclear how to account for both regimes in a semi-analytic fashion.
The sensitivity ofq to a high-energy cut-off scale comes from the underlying assumption that the medium in heavy-ion collisions consists of dressed quasi-particles whose interaction cross-section is given by a Coulomb-like power-law. A first step towards the unification of the two limits described above was undertaken in [1]. The main idea, inspired by the Moliére theory of scattering [23] (see also [24] for a more recent application in the context of momentum broadening in high energy proton-nucleus collisions), is to treat the leading logarithm generated by the Coulomb tail in the medium potential to all orders in opacity since in this case the problem simplifies to solving for a harmonic oscillator potential and the remainder is treated a perturbation. The expansion parameter will thus be ln −1 (q 2 med /µ 2 ), where q med the typical transverse momentum transfer with the medium and µ an IR cut-off such as the Debye mass.
In this work, we generalize the approach in [1] to go beyond the leading-logarithmic corrections toq. This is achieved by expanding the scattering kernel around a harmonic oscillator approximation to incorporate the effects of hard scattering on top of multiple soft interactions. This is equivalent to a shift of the conventional opacity expansion around vacuum propagation to a solution that directly accounts for multiple soft scattering in the medium. Since the procedure in principle encodes the full information about the power-law behavior of Coulomb interactions, it also describes the Bethe-Heitler regime of very soft gluon emissions, i.e. t f L, which is especially relevant for dilute media, mfp L. Our result for the spectrum of medium-induced splittings is therefore valid for arbitrary energies and medium sizes, interpolating between three regimes: 1) Bethe-Heitler (t f < mfp ), 2) LPM ( mfp < t f < L) and 3) single, hard scattering (L < t f ). Although we implicitly assume that the jet is energetic enough such that E/q > L (the "thin" medium limit, according to [20]), our improved formula also accounts for "thick" media where L > E/q.
The manuscript is organized as follows. In Sec. 2, we recall the general expression for leading order splitting distribution in the presence of a dense QCD medium. In Sec. 3, we evaluate the spectrum by expanding close to the Harmonic Oscillator. We carry out the analytic calculations for the first two terms that are sufficient to interpolate between multiple soft and single hard approximations. In Sec. 4, our results for the spectrum and the rate are plotted. The latter is compared to the full numerical results from the McGill group which was first presented in [2].

Spectrum of medium-induced gluons
The probability for a high energy parton a, of energy E, to split into a two partons b and c carrying a fraction z and 1 − z of its energy, respectively, due to multiple scatterings in a dense QCD medium is given by (2.1) where P b(c)←a (z) stands for the Altarelli-Parisi splitting functions that read for the various relevant branching processes (3). N f is the active number of quark flavors which we fix to 3 for our applications throughout. Note that we have adopted the instant form notations for simplicity. The time variables stand in fact for light-cone time, i.e., t = x + and the the mass (energy) E corresponds the the longitudinal light-cone momentum p + .
The Green's function K(x, t 2 ; y, t 1 ) accounts for the interactions with the medium taking place during the formation time of the splitting. It implicitly depends on the color representation of the in-coming and out-going partons and obeys the following Schödinger equation, where medium interactions are incorporated via the interaction Hamiltonian H int = iv(x, t).
The free part of the non-relativistic Hamiltonian in 2+1 dimensions is given by H 0 = i∂ t + ∂ 2 /(2ω) with a "mass" parameter given by ω ≡ z(1 − z)E is solved by the following vacuum Green's function, Naturally, for z < 1 we can identify ω with the energy of the emitted, soft daughter particle, ω ≈ zE. Thus, Eq. (2.6) describes the propagation through, and subsequent transverse broadening in a medium described by an imaginary three-body potential iv ba (x, t), which is given by where C a , C a and C c are the color factors associated with partons in representations a, b and c, respectively. Explicitly, we have The elastic scattering potential with the medium can be extracted from thermal field theory in a weakly-coupled medium [25][26][27], but is often modeled as (2.14) also referred to the Gyulassy-Wang model [28] (here, n = n(t) ∼ T −3 corresponds to the density of scattering centers in the medium). The potential is screened at the scale µ that is related to the Debye mass in a thermal medium, i.e., µ 2 ∼ m 2 In the case of interest many scattering centers contribute during the branching process and as a result the typical transverse momentum acquired by the three body system (a, b, c) is much larger than the Debye mass k ⊥ ∼ x −1 ⊥ µ. As a result,ṽ(x, t) is dominated by a large Coulomb logarithm ln(x ⊥ µ) −1 . Using the GW model for the potential given by Eq. (2.14), we find up to sub-leading corrections of order ∼ x 4 . Here, K 1 (x) is the modified Bessel function of the second kind and γ E ≈ 0.577 . . . the Euler-Mascheroni constant. For later convenience we have introduced the "bare" quenching parameter, stripped of the Coulomb logarithm For a thermal medium we simply have n = (3/2)T 3 . For T = 0.4 GeV and g = 1.94 (so that α s ≈ 0.3) we findq 0 (t) = 1.83 GeV 2 /fm and m D = 0.9 GeV.

Expanding around the harmonic oscillator
It is in general difficult to solve Eq. (2.6) exactly besides using numerical methods [21,22,27]. A first strategy consists in a plain expansion in v(x, t), which stands for the standard opacity expansion, where opacity is defined as a ratio between the medium length and the mean free path ∼ mfp /L. The first order in the latter approach is typically referred to as the Gyulassy-Levai-Vitev (GLV) approximation [12], see also [11]. However, this approximation breaks down in a dense medium at frequencies ω < ω c ∼qL 2 , which for realistic values such asq = 2 GeV 2 /m and L = 5 fm for instance yields a large value ω c = 250 GeV. In this case, the transverse momentum accumulated during the branching is determined by multiple soft scattering, i.e. k 2 ⊥ ∼ √q ω. In this non-perturbative regime the potential can be approximated by v(x, t) ∼ x 2 by neglecting the variation of the Coulomb logarithm (cf. Eq. (2.15)). In this case, the equation of motion is identical to that of a harmonic oscillator with complex frequency, hence this scheme is often referred to as the "harmonic oscillator" (HO) approximation. Of course, in order to obtain a quantitively sound result one needs to estimate the argument of the logarithm which introduces an uncertainty which is of order the inverse of the Coulomb logarithm.
Our strategy in what follows is to shift the expansion point to the "harmonic oscillator" solution as follows , will be treated as a perturbation. The HO potential is given by with the effectiveq parameter chosen to bê whereq 0 (t) is given in Eq. (2.17) and a = 4e −2γ E +1 . Here, we have limited our discussion to a parton in representation R = q, g that radiates a gluon of energy fraction z. The constant a accounts for the constant terms in Eq. (2.15). The jet quenching coefficientq eff is logarithmically dependent on a subtraction scale which is only a function of ω, i.e. Q sub = Q sub (ω). This scale must be chosen to be the typical transverse momentum generated during the splitting, that is, in the HO approximation. The advantage of choosing this expansion point, is that in the HO approximation the Schödinger equation, admits the known analytic solution [29], Here, the functions S(t, t 0 ) and C(t, t 0 ) represent the two independent solution of the equation where the frequency is given by Ω(t) = q(t)/(2iω), with boundary conditions S(t 0 , t 0 ) = 0 and ∂ t S(t, t 0 )| t=t 0 = 1, and C(t 0 , t 0 ) = 1 and ∂ t C(t, t 0 )| t=t 0 = 0, respectively. They are related through a constant Wronskian, In what follows we shall solve the latter iteratively for the first two orders: the leading-order (LO) term reads K LO = K (0) , where The next-to-leading (NLO) correction is given by so that the full NLO term is simply K NLO = K (0) + K (1) . Higher orders are found in a analogous manner. Our results below show that the two first terms already provide a reasonable description of the spectrum and the related splitting rate.

Leading order: the harmonic oscillator approximation
Let us first consider the leading order that corresponds the BDMPS approximation. Inserting Eq. (3.6) in Eq. (2.1) yields where we have left the indices indicating the parton splitting to be implicit. Using the following property [30] ∂ t C(t, s) S(t, s) the t 2 integration can be carried out and reads . (3.14) The first term in Eq. (3.14) cancels against the vacuum piece, i.e. the second term in Eq. (3.12), while the second one can be integrated further over t 1 , where we have used the decomposition of C(∞, s) and S(∞, s) as a superposition of other solution to the wave equation [30], Recall that q eff is a function of z and depends on the flavor of the partonic configuration (cf. Eq. (3.3)).

Next-to-leading order correction to the harmonic oscillator
Let us turn now to evaluating the next-to-leading term, given as a sum of Eqs. (3.10) and (3.11). The physical meaning of the NLO correction is that one "soft" scattering, described purely via diffusive transverse momentum broadening, is replaced by a "hard" scattering, i.e. an interaction with the medium described by the Coulomb potential. This opens for the possibility that a single, hard kick from the medium can dominate the total transverse momentum transversed from the medium during the formation time. The correction to the splitting distribution reads, see [1] for more details, where the indices are suppressed and In particular, for a medium with constant density n(s) = nΘ(L − s), we find For simplicity and without loss of generality, let us focus on the case where a gluon of energy fraction z is emmitted of a parton R. Then, In order to integrate over u it is convenient to use the Fourier representation of the dipole cross-section. Consider for instance the contribution from the second term in Eq. (3.24), The contribution a relates to the compete dipole cross-section from which the HO part, i.e. contribution b, must be subtracted. Strikingly, the integrations can be performed analytically by noticing that the u integral yields Furthermore, by changing variables to y = z 2 q 2 /(4k 2 (s)) we find where the first term appears due to the full elastic cross section, cf. Eq. (3.26), while the second term is the subtraction of the harmonic oscillator term, cf. Eq. (3.27). The first and third remaining terms can be found by simply substituting z → 1 and z → 1 − z, respectively. Explicitly, the full NLO correction then takes the form where we introduced the shorthand The full spectrum at NLO therefore becomes the sum of Eqs. (3.18) and (3.31). Let us investigate the two limits of the complex function k 2 (s), given in Eq. (3.23). In the limit ω qL 2 , it follows that Ω 1, leading to tan(Ω(L−s)) ≈ 0 and cot Ωs ≈ (Ωs) −1 . It follows that as in the vacuum. Then, expanding for small x = µ 2 /[4k 2 (s)], we find Hence, the spectrum in the high frequency regime ω qL 2 becomes which is a well-known limit of the GLV spectrum.
Turning to the small frequency regime, note that k 2 (s) also becomes vacuum-like since Q 2 sub ≈ µ 2 for ω < ω BH and thereforeq → 0. We can now expand for large x, to find This is the characteristic behavior of the Bethe-Heitler regime, that is also contained in the GLV spectrum.
It is worth pointing out that these limits are universal and do not depend on the choice of matching scale Q sub inq. The approach to these values is however affected by the exact value of the matching scale, which we shall explore in the next section.

Numerics
Let us first return to the subtraction scale Q sub introduced earlier in the context the transport coefficientq =q 0 ln aQ 2 sub /µ 2 . As discussed in [1,27], it is natural to define it in relation to the characteristic transverse momentum of the medium-induced emission process, i.e. Q 2 sub ∼ k 2 ⊥ . For radiative processes, we have that k 2 ⊥ ∼ (ωq) 1/2 ∼ (ωq 0 ln a √ ωq 0 /µ 2 ) 1/2 , sinceq itself is running with the subtraction scale. 4 Crucially, to ensure the matching with the Bethe-Heitler spectrum at small gluon frequencies the latter logarithm should vanish for Q 2 < µ 2 . This implies thatq → 0 when √ ωq < µ 2 or ω < ω BH ≡ µ 2 mfp (where we used thatq ∼ µ 2 / mfp and we are left with a vacuum spectrum in this regime. To do so we use the following interpolating form This choice guarantees that Q sub smoothly goes to µ. Finally, in order to test the sensitivity to the chosen matching, we will multiply the numerical factor a by 1/2 and 2. Let us first turn to the eikonal limit, where assume x 1. We compare the results from our formula in the eikonal limit with the expectation from first order in opacity (for short, labelled "GLV") and the multiple soft scattering approximation (labelled "BDMPS") in Figure 1. The spectrum interpolates well between the three physical regimes for inmedium QCD bremsstrahlung. For the choice of parameters here, corresponding to mfp < L < E/q, it clearly demonstrates that LPM interference effects are suppressing the spectrum over a large range of gluon energies. It is worth keeping in mind that we have only included the first order correction to the standard HO baseline, which leaves further room to improve on the matching by adding higher orders. The band around the "LO+NLO" curve corresponds to varying the matching parameter a by a factor 2 up and down, which describes the inherent ambiguity in defining the LPM suppressed regime. Most importantly, our results reproduce the universal expectations at low-(corresponding to the Bethe-Heitler limit) and the high-frequency (corresponding to the "GLV" limit) and, all in all, the result over a wide frequency range is quite good.
The results for the improved spectrum with full x dependence, calculated for a gluon jet with for three different energies E = {1000, 250, 62.5} GeV, are given in Figure 2. For the upper energy, the jet traverses a "thin" medium where E/q > L while for the lower energy, the medium is "thick", i.e. E/q < L. This affects the range where the LPM is at play. The turning over of the curves at small-x correspond to the BH regime which takes place at a fixed gluon energy ω BH , not a fixed x for different jet energies.
We have also computed numerically the rate of emissions, defined as dI/(dω dt), in Figure 3. As expected, at small times the rate grows linearly with time like the GLV spectrum, i.e. ∝ t. At later times the rate satures in the LPM regime. We observe the uncertainty related to the implementation of the LPM regime as a band that predominantly  appears whenever the rate saturates. For the same reason, we have therefore chosen not to plot the expectation from the "bare" BDMPS spectrum since it would differ from our  curve by a constant offset related to choice of scale inq. Finally, in order to compare our compact, analytical formula with the full, all-order in opacity solution of the spectrum, solved numerical in [2], we plot in Figure 4 the rate for two choices of medium parameters and jet kinematics. Since our result only includes the first corrections from hard scattering, the agreement is reasonable ( 30%). Note that some of the discrepancy may be attributed to the different choices of the potential and the lack of thermal masses in our approach.

Conclusions
In this work, we revisit the calculation of medium-induced parton splitting and present a new analytic approach that allows for the first time to account for the various known limits. Our method, dubbed "Improved Opacity expansion" resums multiple soft scatterings to all orders while treating single hard scattering as a perturbation. To do so, following [1], we have suggested a shift of the expansion of the in-medium propagators around the so-called "harmonic oscillator" solution which takes into account diffusive momentum broadening. Perturbations around this solution correspond to hard, transverse kicks that reveal the quasi-particle structure of the underlying medium. In Ref. [1], the radiative spectrum was calculated in the leading logarithmic approximation and therefore is applicable so long as x ⊥ µ −1 , which translates into ω ω BH . In the present work we go beyond by accounted for the full imaginary potential order by order.
We have demonstrated the validity of the framework by computing the in-medium radiative spectrum. The NLO term, with a suitable choice of subtraction scale, allows to properly link the LPM regime, appropriate for dense media, to the regimes where singlescattering in the medium dominates, including the Bethe-Heitler regime at low frequencies, where formation time of the radiation probes the scale of the medium mean-free-path, and the GLV regime at high frequencies, where the formation time of the bremsstrahlung  exceeds the length of the medium. This demonstrates that the approach, albeit formally suitable for large and dense media, where diffusive broadening dominate at small angles. The compactness of our main result and the good agreement with the results from an allorder resummation [2], which is implemented in the MARTINI event generator [31,32], makes it amenable for implementation in a fast Monte-Carlo event generator for jet quenching.