A complete set of in-medium splitting functions to any order in opacity

In this Letter we report the first calculation of all ${\cal O}(\alpha_s)$ medium-induced branching processes to any order in opacity. Our splitting functions results are presented as iterative solutions to matrix equations with initial conditions set by the leading order branchings in the vacuum. The flavor and quark mass dependence of the in-medium $q \rightarrow qg$, $g\rightarrow gg$, $q \rightarrow g q$, $g \rightarrow q\bar{q}$ processes is fully captured by the light-front wavefunction formalism and the color representation of the parent and daughter partons. We include the explicit solutions to second order in opacity as supplementary material and present numerical results in a realistic strongly-interacting medium produced in high center-of-mass energy heavy ion collisions at the Large Hadron Collider. Our numerical simulations show that the second order in opacity corrections can change the energy dependence of the in-medium shower intensity. We further find corrections to the longitudinal and angular distributions of the in-medium splitting kernels that may have important implications for jet substructure phenomenology.


I. INTRODUCTION
In Quantum Chromodynamics (QCD), the probability of one parton to split into a two-parton system is governed by the well-known Altarelli-Parisi splitting functions in vacuum [1]. These probabilities are the key ingredients in all modern high-precision calculations in the theory of strong interactions and in Monte-Carlo event generators for high energy physics. For jet physics, quark and gluon splitting kernels play an essential role in understanding the radius dependence of inclusive cross sections and of jet substructure. These can be directly measured using novel observables, such as the groomed soft-dropped momentum sharing distribution of the two leading subjets inside a reconstructed jet [2]. For jets propagating through hot or cold QCD matter, the vacuum splitting functions receive medium-induced contributions [3][4][5][6][7][8][9][10] where non-local destructive interference, known as the non-Abelian Landau-Pomeranchuk-Migdal (LPM) effect [11][12][13][14][15][16], plays an important role. These modifications can then be calculated order by order in powers of the opacity χ, the average number of scatterings in the medium. The basis of the opacity expansion and the current state of the field are reviewed in Ref. [17] and the references therein.
So far, a realistic numerical evaluation of higher order-in-opacity corrections has been missing from the literature. Early works investigated the convergence of the opacity series to second and third order, utilizing the soft gluon emission limit available at the time [14,18]. These studies used models of the QCD medium density carefully chosen to smooth out the non-Abelian LPM interference pattern, making certain integrals analytically calculable. It was found that for opacities as large as five the series converges quickly and the corrections to parton energy loss can be of the order of a few percent, especially at typical jet energies available at the LHC. The energy spectrum of radiated gluons receives larger corrections, of order 10 -30%, and the fully differential angular distribution of the radiation can be affected even more. Similar conclusions were reached in Ref. [19], where corrections up to ninth order in opacity were evaluated, albeit with much reduced numerical accuracy. In the high-twist approach second order corrections were calculated in Ref. [20]. The possibility of a hard scattering was added to the infinite medium soft scattering soft scattering approximation [21]. Importantly, it was recently shown in the work of Feal and Vasquez [22,23] that analytic approximations to in-medium splitting processes are not a good substitute for their full numerical evaluation.
In Ref. [17] we constructed a new framework to compute these splitting functions to higher orders in the opacity expansion, which was previously limited to either the soft-gluon approximation x 1 or to first order in opacity. With the derivation of the full in-medium splitting functions and theoretical formalisms that bridge the gap between high energy and heavy ion physics, the question of corrections from higher orders in opacity must now be revisited.
We are also motivated by the more precise experimental measurements of jets and jet substructure that have become available or are expected in the near future, for recent examples see [24][25][26][27][28][29][30][31]. In-medium splitting kernels enter directly into the fixed order and resummed calculation of such observables [32][33][34][35][36]. In this Letter, we present new results which extend the formalism constructed in Ref. [17] to allow the calculation of all partonic splitting channels beyond the soft-gluon approximation to any order in opacity, and we evaluate them numerically at second order in opacity for the first time. Unlike the early works which used medium models chosen for theoretical convenience, here we employ realistic state-of-the-art hydrodynamic models of the expanding quark-gluon plasma to enable precision jet phenomenology.
For the purposes of this Letter we will limit ourselves to the final-state interactions of a jet produced in the medium, and we will work in the limit where the energy of the jet is much larger than any medium scales such that collisional energy loss is suppressed compared to radiative branching. Throughout this Letter we will use light-front coordinates We work in a general frame in which the jet moves along the +z axis with longitudinal momentum p + ; all the results are manifestly boost invariant, and the medium rest frame is just the special case for which p + = √ 2E with E the jet energy. The rest of our Letter is organized as follows: in Section II we outline the analytic derivation of all medium-induced splitting functions. Their numerical evaluation is discussed in Section III, and explicit new results to second order in opacity are given. We give our conclusions in Section IV.

II. UNIVERSAL FEATURES OF MEDIUM-INDUCED PARTON BRANCHING
For the four different branching channels in QCD to lowest order, we will write a → bc with partons a, b, c being either (anti)quarks or gluons. We take the parent parton a to have large longitudinal momentum p + along the jet axis, with the daughter partons b and c carrying fractional momenta (1 − x)p + and xp + , respectively. We also denote by k the relative transverse momentum given to parton c by the branching. As was shown for the q → qg branching in Ref. [17], the splitting functions can be expressed in terms of the following quantities: the light-front wave function (LFWF) of the branching ψ [37,38], the corresponding energy denominator ∆E − , the cross-section dσ el d 2 q to interact with a scattering center in the medium, and the scattering length λ of a parton. These features are universal to all partonic branching channels a → bc and to any QCD medium; the LFWF themselves will differ for different channels, but the way in which they are dressed by the medium is the same.
The LFWF (squared and averaged over quantum numbers) and energy denominator for any channel can be written in the form with functions f (x), g(x) describing the distribution of longitudinal momentum in the branching and a coefficient ν accompanying the mass dependence. These functions and coefficients will be fully listed in Table I below. The LFWF describe the branching in vacuum, that is, at 0 th order in opacity O χ 0 : where C 0 is the color factor for the vacuum splitting and p + dN0 d 2 p dp + is the distribution of the parent parton with momentum (p + , p) created by the hard scattering in the medium. As in Ref. [17], each additional order in opacity corresponds to an additional correlated rescattering in the medium and generates a corresponding integral over both the position z + of and momentum exchange q with the scattering center: R is the mean free path of a parton in color representation R (λ + q for a quark or λ + G for a gluon) and σ el = d 2 q dσ el d 2 q is the elastic scattering cross section on a scattering center. Local color neutrality of the scattering centers implies that the leading-order scattering process involves two-gluon exchange at the level of the amplitude squared, corresponding to 9 "direct" scattering diagrams and 8 "virtual" scattering diagrams, labeled D 1−9 and V 1−8 respectively in the notation of Ref. [17]. The color factors for each of these diagrams are specified by coefficients d 1−6 and v 1,3,5,7 that differ for each channel, but they are preserved and iterated at successive orders in opacity. 1 Each partonic splitting channel is then fully specified by the quantities summarized in Table I. The choice of representation R used in the mean free path is arbitrary, but this choice rescales the scattering color factors d 1−6 and v 1,3,5,7 ; here we adopt the convention of using the representation of the parent parton for each channel. Note that the q → gq channel is obtained from the q → qg channel under the mapping x → (1 − x) , k → −k. This transformation leaves the splitting function invariant as long as the cross section dσ el d 2 q is even under q → −q. For finite kinematics, we can expect some deviations from this symmetry.
As derived in Ref. [17], the contribution from N th order in opacity to the distribution of partons within a jet can be expressed as The function f F/F is computed using the generalized recursion relations (referred to as the reaction operator) These equations can also be conveniently written in an upper-triangular matrix form [17]. Though cumbersome, they are straightforward to implement order by order in the opacity expansion in algebraic manipulation packages such as Mathematica [39]. The reaction operator Eq. (4) encodes both the induced splitting and the momentum broadening of the total jet transverse momentum to a final observed value p from some initial value p − q tot . In the broad source approximation, one neglects the shift q tot compared to the hard scale in the initial distribution p + dN0 d 2 (p−qtot) dp + of hard partons produced in the medium. Then, the medium-induced splitting function can be considered to depend only on x and k through the ratio: dN d 2 k dx ≡ p + dN d 2 k dx d 2 p dp + /p + dN0 d 2 p dp + . We have developed a Mathematica code to solve the recursion relations Eqs. (4) with initial conditions Eqs. (5) to any order in opacity for all of the four lowest order splitting kernels. We checked that our results recover the known expressions in the literature for both massless and massive quarks, and gluons to first order in opacity [5,10]. In the soft gluon emission limit we further verified the that our results for the diagonal splitting functions to second order in opacity coincide with the solutions reported in [14,18]. The interested reader will find the full second order in opacity corrections to the medium-induced splitting kernels that we calculate here for the first time as supplementary material.

III. RESULTS
To evaluate the recursion relations Eqs. (4) for the medium-induced splitting functions, we employ the broad source approximation in a frame such that the initial jet momentum relative to the direction of propagation is p T = 0 and consider a particular model for the medium, assumed here to be a quark-gluon plasma with partial densities for relativistic massless quarks and gluons given by Here, T (x, t) is the temperature that depends on the time t and the position in the transverse plane of the collision x, N g DOF = 16, and for two light quark flavors N f = 2 we get N q DOF = 24. With the Debye screening scale µ 2 D (x, t) = g 2 T 2 (x, t)(1 + N f /6), the jet on medium quasiparticle scattering cross sections [40] We use a typical value of the in-medium coupling constant g = 2 for illustration. The upper limits of the q T integrals in Eqs. (4) are given by q max iT = √ µ D i E and effective thermal mass of partons is included by replacing Here, the index i corresponds to position z i in the nested path length integrals along the direction of jet propagation, N of them to O(χ N ). We also go from light-front to Cartesian coordinates, for example z + /E + = z/E and the dz + /λ + R = dz/λ R . For the upper limit of these line integrals we take a length much bigger than the size to the QCD medium, such that the contributions outside of the quark-gluon plasma naturally vanish by ρ → 0, λ R → ∞. The scattering length λ R we obtain as 1/λ R = σ el Rq→Rq ρ q + σ el Rg→Rg ρ g , where R can be a quark or a gluon.
The only further input needed to describe the thermal medium and compute the in-medium splitting functions is the spacetime temperature profile along the jet propagation axis. For an expanding quark-gluon plasma, this profile is set by using the 2+1D event-by-event viscous hydrodynamics code iEBE-VISHNU [41] with Monte-Carlo-Glauber initial conditions at an initial time of τ 0 = 0.6 fm. The hydrodynamical evolution uses the S95-1 Equation of State [42] to generate the temperature grids which are used to evaluate Eqs. (7) and the relevant medium quantities. The jets themselves are initiated randomly according to the number density of binary collisions.
Numerical values of the splitting functions are then obtained by evaluating the integrals given in Eq. (4) in the broad source approximation using the VEGAS algorithm [43]. Compared to the leading order calculations, the second order calculations require integration of additional three dimensions for z 2 and q 2 . As the integration dimension becomes larger, the convergence of the VEGAS algorithm that approximates the probability distribution function of the integrand becomes slower. Therefore, in the calculation of the second order corrections we use a 10 times larger number of random samples than that used in the leading order calculation, in order to make the statistical uncertainties of the correction terms similar to those of the leading order results. In addition, the integrand of the second order correction term requires about 10 times larger number of floating point operations than the leading order calculation.
We are now ready to present our numerical results and first define the integral of the splitting kernel weighed by the energy fraction of the emitted parton over a finite range x min ≤ x ≤ x max as follows with the transverse momentum of the splitting Λ 2 QCD ≤ k 2 T ≤ 4E 2 x(1−x). If we restrict the energy fraction Eq. (8) to only the diagonal channels q → gq and g → gg and the soft gluon limit x 1, it can be interpreted as the fractional The one-dimensional differential splitting functions dN dx for a 100 GeV jet as a function of the splitting fraction x. The symbol and color scheme is the same as in Fig. 1, and so is the underlying QGP medium. Computations for both the diagonal and off-diagonal branchings are presented. radiative energy loss ∆E E of the parent parton due to the emission of soft gluons. In the general case considered above, the concept of "parton energy loss" is ill-defined but I xmax xmin gives us a sense of the corrections due to second order in opacity to inclusive observables, such as inclusive hadron or inclusive jet suppression in nucleus-nucleus collisions. Results are shown for lead-lead (Pb+Pb) collisions at √ s N N = 2.76 TeV in Fig. 1 for the diagonal spitting channels mentioned above. The integral Eq. (8) is plotted in Fig. 1 for the 0 − 10% centrality bin, using cutoffs For light quark jets and gluon jets, a smoothing of the energy dependence is seen, with reduced in-medium shower intensity at lower jet energies and increased in-medium shower intensity at higher jet energies. For jets with energies around ∼ 100 GeV, the corrections from second order in opacity appear to be modest, of order ∼ O(10%), while the corrections at both higher and lower jet energies are larger. Based on the calculations of higher orders in opacities in the soft gluon approximation [14,18,19], the corrections tend to alternate in sign, such that the first and second orders in opacity bound the resummed result both above and below, approximating an error band associated with the opacity expansion. We show by a gray line the average of the 1 st and 1 st +2 nd orders in opacity, which we expect will be suitable for phenomenological implementations. For b-quarks the most pronounced feature remains the reduction of the in-medium shower intensity for small jet energies due to the heavy quark mass effects, already seen in the soft gluon approximation [44,45].
The one-dimensional spectrum of the parton splitting dN dx is shown in Fig. 2 for a 100 GeV jet. For light partons the corrections from second order in opacity in the small x and large x regions are all negative and sizeable, of order 30% − 70% depending on the channel. Interestingly, however, in all channels in the region of democratic branching x ∼ 0.5 the second order correction is positive, enhancing the branching probability and thereby producing a smoother x dependence. We have checked that at higher jet energies this region grows considerably, covers most of the phase space in x and is entirely responsible for enhanced in-medium shower intensity relative to the lowest order in opacity. For heavy quarks 2 nd order in opacity corrections are noticeably smaller in most parts of phase space.
Finally, the transverse momentum distribution of the medium-induced radiation is illustrated in Fig. 3 for a 100 GeV jet at x =0.3. We have performed simulations for several centralities and choose the 10 − 20% centrality bin as an example. Here, the medium-induced spectrum dN med is expressed to O χ N as a ratio to the vacuum radiation Eq. (2) to highlight the changes induced by the medium. Note that the medium-induced part can be negative, since it excludes the vacuum contribution O χ 0 . Shown in Fig. 3 are the spectra at 1 st and 1 st + 2 nd orders in opacity, along with the average of the two to approximate the resummed curve due to the oscillatory nature of the opacity expansion as discussed above. At small k T the 2 nd order corrections can alter the angular spectrum considerably. At larger k T the second-order corrections are milder, reducing the peak region from first order and gradually dying off around k T ≥ 5 GeV. All corrections for the g → qq channel are much smaller and do not change the sign of the spectrum.

IV. CONCLUSIONS
In this Letter we presented solutions for the q → qg, g → gg, q → gq, g → qq in-medium splitting kernels beyond the soft gluon approximation to O(α s ) and to any desired order in opacity. Our formalism relies on the universal features that such branching processes exhibit when they are expressed as solutions to iterative equations that build correlations between multiple scattering centers in the medium. The differences between the splitting channels are conveniently captured by the light-front wavefunctions and coefficients that depend on the color representation of the partons that interact with the medium. As the first order in opacity branchings are well-documented in the literature [5,10], we only included second order in opacity corrections explicitly as supplementary material, together with the process-dependent coefficients and functions in Table I. The interested reader can recover our results and obtain even higher order in opacity corrections if desired.
We also demonstrated that the full process-dependent in-medium splitting kernels can be evaluated numerically to higher orders in opacity in a realistic nuclear medium. In this Letter such a medium was exemplified by a quark-gluon plasma produced in high energy Pb+Pb collisions at the LHC. We presented new numerical results which explore for the first time the corrections due to second order in opacity with exact kinematics for all partonic branching channels in both the massive and massless quark limits. As anticipated, we found the largest corrections for the most differential quantity, the transverse momentum spectrum of the medium-induced splitting functions, with these corrections being smoothed over as we move to the more integrated quantities.
The derivation and numerical evaluation of higher order in opacity corrections is important and timely, as their real application is in medium-induced splitting kernels that enter the calculation of fixed order and resummed observables involving reconstructed jets. Corresponding to the changes we see in the more differential splitting functions, we anticipate that the greatest phenomenological impacts will be on the interpretation of jet substructure measurements in heavy ion collisions [24][25][26][27][28][29][30][31]. Even in more inclusive observables, such as the suppression of jet production denoted R AA , we anticipate the second order in opacity corrections to play a role, albeit milder. Specifically, they may be relevant for the theoretical understanding of the weak transverse momentum dependence of the several-hundred GeV jets R AA observed by the ATLAS collaboration [30].
The splitting functions computed here are now ready for broad phenomenological applications. The very general formulation presented here is also amenable to scanning different collision systems across system size, deformation, and other parameters [46]. In future work, we will pursue this new program of jet phenomenology at second order in opacity to make contact with the realistic values of opacity achieved in typical heavy-ion collisions. A future extension of this numerical analysis to the third order in opacity would also be a useful validation of the oscillatory nature of the opacity series and the use of the first and second orders as upper and lower bounds for the medium-induced splitting functions.