Moli\`ere Scattering in Quark-Gluon Plasma: Finding Point-Like Scatterers in a Liquid

By finding rare (but not exponentially rare) large-angle deflections of partons within a jet produced in a heavy ion collision, or of such a jet itself, experimentalists can find the weakly coupled short-distance quark and gluon particles (scatterers) within the strongly coupled liquid quark-gluon plasma (QGP) produced in heavy ion collisions. This is the closest one can come to probing QGP via a scattering experiment and ultimately learning how a strongly coupled liquid emerges from an asymptotically free gauge theory. The short-distance, particulate, structure of liquid QGP can be revealed in events in which a jet parton resolves, and scatters off, a parton from the droplet of QGP. The probability for picking up significant transverse momentum via a single scattering was calculated previously, but only in the limit of infinite parton energy which means zero angle scattering. Here, we provide a leading order perturbative QCD calculation of the Moli\`ere scattering probability for incident partons with finite energy, scattering at a large angle. We set up a thought experiment in which an incident parton with a finite energy scatters off a parton constituent within a"brick"of QGP, which we treat as if it were weakly coupled, as appropriate for scattering with large momentum transfer, and compute the probability for a parton to show up at a nonzero angle with some energy. We include all relevant channels, including those in which the parton that shows up at a large angle was kicked out of the medium as well as the Rutherford-like channel in which what is seen is the scattered incident parton. The results that we obtain will serve as inputs to future jet Monte Carlo calculations and can provide qualitative guidance for how to use future precise, high statistics, suitably differential measurements of jet modification in heavy ion collisions to find the scatterers within the QGP liquid.


Introduction
When the short-distance structure of quark-gluon plasma is resolved, it must consist of weakly coupled quarks and gluons because QCD is asymptotically free. And yet, at length scales of order its inverse temperature 1/T and longer, these quarks and gluons become so strongly correlated as to form a liquid. Heavy ion collisions at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) produce droplets of this liquid QGP whose expansion and cooling is well described by relativistic viscous hydrodynamics with an unusually small viscosity relative to its entropy density. (For reviews, see Refs. [1][2][3].) This discovery poses a question: how does this strongly coupled liquid emerge (as a function of coarsening resolution scale) from an asymptotically free gauge theory? In other contexts, the path to addressing a question like this about some newly discovered complex strongly correlated form of matter would begin with doing scattering experiments, and in particular would begin with doing scattering experiments in which the momentum transfer is large enough that the microscopic constituents (in our case, weakly coupled at short distance scales) are resolved. Some analogue of such high resolution scattering experiments are a necessary first step toward understanding the microscopic structure and inner workings of QGP. Since the droplets of QGP produced in heavy ion collisions rapidly cool and turn into an explosion of ordinary hadrons, the closest that anyone can come to doing scattering experiments off QGP is to look for the scattering of energetic "incident" partons that are produced in the same collision as the droplet of QGP itself. Since such energetic partons shower to become jets, this provides one of the motivations for analyzing how jets produced in heavy ion collisions are modified via their passage through QGP. Pursuing such measurements with the goal of understanding the microscopic workings of QGP has been identified [4][5][6] as a central goal for the field once higher statistics jet data anticipated in the 2020s, at RHIC from the coming sPHENIX detector [7] and at the LHC from higher luminosity running, are in hand. The short-distance, particulate, structure of liquid QGP can be revealed by seeing events in which a jet parton resolves, and scatters off, a parton from the droplet of QGP. If the QGP were a liquid at all length scales, with no particulate microscopic constituents at all, as for example is the case in the infinitely strongly coupled conformal plasma of N = 4 supersymmetric Yang-Mills (SYM) theory, then the probability for an energetic parton plowing through it to pick up some momentum q ⊥ transverse to its original direction is Gaussian distributed in q ⊥ [8][9][10], meaning that large-angle, large momentum transfer, scattering is exponentially (maybe better to say "Gaussianly") rare. The q ⊥ distribution should similarly be Gaussian for the case of an energetic parton plowing through the QGP of QCD -as long as q ⊥ is not too large. One way to see this is to realize that as long as q ⊥ is small enough the energetic parton probes the QGP on long enough wavelengths and "sees" it as a liquid. Another way to reach the same conclusion is to imagine the nottoo-large q ⊥ as being built up by multiple soft (low momentum transfer; strongly coupled) interactions with the QGP. The key point, though, is that in QCD, unlike in N = 4 SYM theory, this cannot be the full story: real-world QGP must be particulate when its short-distance structure is resolved. This means that large-angle, high momentum transfer, scattering may be rare but is not Gaussianly rare, as Rutherford would have understood. So, if experimentalists can detect rare (but not Gaussianly rare) large-angle deflections of jet partons plowing through QGP, referred to as "Molière scattering" after the person who first discussed the QED analogue [11][12][13], they can find its weakly coupled quark and gluon constituents [10,14] and begin to study how the strongly coupled liquid emerges from its microscopic structure.
One idea for how to look for large angle scattering is to look for deflections of an entire jet [10] by looking for an increase in the "acoplanarity" of dijets or gamma-jets (meaning Figure 1. Kinematics of the thought experiment that we analyze. An incident parton of "type" C (type meaning gluon or quark or antiquark) with energy p in impinges on a "brick" of QGP with thickness L. An outgoing parton of type A with energy p is detected at an angle θ relative to the direction of the incident parton. We shall calculate the probability distribution of p and θ for a given p in and for all possible choices A and C.
the angle by which the two jets or the photon and jet are not back-to-back) in heavy ion collisions relative to that in proton-proton collisions. The acoplanarity is already quite significant in proton-proton collisions because many dijets (or gamma-jets) are not back-toback because they are two jets (or a photon and a jet) in an event with more jets. This makes it challenging to detect a rare increase in acoplanarity due to rare large-angle scattering, but these measurements have been pursued by CMS [15,16], ATLAS [17] and ALICE [18] at the LHC and by STAR [19] at RHIC, and it will be very interesting to see their precision increase in future higher statistics measurements. The same study can be done using events with one (or more, unfortunately) jets produced (only approximately) back-to-back with a Z-boson, albeit with lower statistics [20]. It was realized in Ref. [14] that Molière scattering can also be found by looking for rare large-angle scattering of partons within a jet shower, rather than of the entire jet. We shall see that this is advantageous in that it allows one to consider energetic partons within a jet with only, say 20 or 40 GeV in energy, whose kinematics allow for larger angle scattering than is possible if one considers the deflection of (higher energy) entire jets. However, the jet substructure observables needed to detect rare large angle scattering of partons within a jet (via measuring their modification in jets produced in heavy ion collisions) are of necessity more complicated than acoplanarity. It is very important that such observables are now being measured [21][22][23][24][25] and analyzed in heavy ion collisions, as it remains to be determined which substructure observables, defined with which grooming prescription, will turn out to be most effective. Quantitative predictions for experimental observables, whether acoplanarities or substructure observables, require analysis of jet production and showering at the level of a jet Monte Carlo, first in proton-proton collisions and then embedded within a realistic hydrodynamic model for the expanding cooling droplet of matter produced in a heavy ion collision. We shall not do such a study here; our goal is to provide a key theoretical input for future phenomenological analyses, not to do phenomenology here. Nevertheless, we expect that at a qualitative level our results can provide some guidance for planning experimental measurements to come.
In this paper, we set up a thought experiment in which we "shoot" a single energetic parton (quark or antiquark or gluon) with initial energy p in through a static "brick" of QGP of thickness L in thermal equilibrium at a constant temperature T , c.f. Fig. 1. For simplicity, we shall model the medium within our brick as a cloud of massless quarks and gluons, with Fermi-Dirac and Bose-Einstein momentum distributions, respectively. This is surely only of value as a benchmark. Although treating the partons as massless is appropriate if the momentum transfer is high enough, as we shall quantify in Section 3.3, adding thermal masses would surely be a worthwhile extension of our study. Also, our calculations could be repeated in future using any proposed model for the momentum distributions of the quarks and gluons as seen by a high-momentum probe. Indeed, it is hard to imagine a better possible future than the prospect of making experimental measurements that reveal the presence of rare large-angle Molière scattering, seeing quantitative disagreements with predictions obtained via incorporating our calculation within a jet Monte Carlo analysis, and reaching the conclusion that the momentum distributions of the quarks and gluons seen by a high-momentum probe differ from the benchmark distributions that we have chosen.
We shall then compute F (p, θ), the probability distribution for finding an outgoing hard parton with energy p and angle θ relative to the direction of the incident hard parton. We choose to normalize the distribution F (p, θ) as π θ min dθ ∞ p min dp F (p, θ) = N hard (θ min ) , (1.1) where N hard (θ min ) denotes the number of outgoing hard partons in a specific region of the phase space θ ≥ θ min , p ≥ p min per single incident parton. We have introduced a somewhat arbitrary hard energy scale p min so that we can refer to a parton with p > p min as a hard parton. We will specify p min as needed in Sec. 3, and will always choose p min to be significantly greater than T . F (p, θ) will depend on the temperature of the plasma, T , on the energy of the incident parton, p in , on the time that the parton spends traversing the brick of QGP, ∆t ≡ L/c, as well as on whether the incident parton and the outgoing parton are each a quark, antiquark or gluon, but we shall keep all these dependences implicit in our notation in this Introduction. It should be evident that our thought experiment is only that. The droplet of QGP produced in a heavy ion collision expands and cools rapidly; its dynamics is certainly not that of a constant temperature static brick. And, a jet shower is made up from many partons and has a complex showering dynamics of its own. In order to do phenomenology, our results for F (p, θ) must be incorporated within a Monte Carlo calculation of jet production and showering, with the jets embedded within a realistic hydrodynamic description of a droplet of QGP. Such a future calculation, in which the dynamics of a jet (including the splitting and propagation) and of the droplet of plasma is described ∆t by ∆t by ∆t, for some small value of ∆t, after each ∆t our result for F (p, θ) could be applied to each parton in the shower. In this way, our results can be used to add large-angle Molière scattering to a jet Monte Carlo calculation which does not currently include it, like for example the Monte Carlo calculations done within the hybrid model in Refs. [26][27][28][29]. In the case of a Monte Carlo calculation in which hard two-to-two scattering is already included, for example those done within JEWEL [30][31][32][33][34][35], MARTINI [36] or LBT [37][38][39], our results can be used in a different way, namely as a benchmark against which to compare for the purpose of identifying observable consequences of large-angle scattering. The other way in which the results of our calculation will be of value is as a qualitative guide to experimentalists with which to assess how large the effects of interest may turn out to be, namely as a qualitative guide to what the probability is that a parton with a given energy in a jet could scatter by an angle θ. In Section 3.4 we shall illustrate our results by plotting what we obtain for partons with p in = 25 T = 10 GeV and p in = 100 T = 40 GeV and p in = 250 T = 100 GeV incident on a brick with T = 0.4 GeV and ∆t = 3 fm.
Although we believe that our results will be of value as a qualitative guide for planning and assessing future experiments, giving a sense of just how rare it should be for a parton in a jet to scatter at a large enough angle that the jet grows a new prong that can be discerned via high-statistics measurements of suitably defined jet substructure observables, there should be no illusion that this will be a straightforward program. We do not anticipate any smoking guns to be found. As an object lesson, it is worth considering the question of how to detect evidence, in experimental data, for the Gaussian distribution of transverse kicks q ⊥ that all the partons in a jet must pick up as they traverse the plasma. As we noted above, the probability distribution for small q ⊥ is Gaussian, with a width often denoted byqL, after passage through plasma over a distance L and this can be understood either via holographic calculations at strong coupling or as a consequence of multiple scattering in a weakly coupled picture. Constraints on the measured value ofq all come from comparing calculations of energy loss (not transverse kicks themselves) to experimental data on observables that are sensitive to energy loss within a weakly coupled formalism in whichq also controls parton energy loss [40]. There is at present no clear experimental detection of the Gaussian distribution of transverse kicks themselves. The natural way to look for them is to look for increases in the angular width of jets, jet broadening, due to propagation through plasma, as all the partons in a jet accumulate Gaussian-distributed transverse kicks. In fact, it is with this in mind that these kicks are typically referred to as transverse momentum broadening. There are many extant measurements of the modification of jet shape observables in heavy ion collisions [18,19,21,23,25,[41][42][43][44][45], and many theorists have made efforts to turn these measurements into constraints on transverse momentum broadening, for example see Refs. [28,35,37,[46][47][48][49][50], but there are two significant confounding effects that obscure transverse momentum broadening [28]. The first effect is that the energy and momentum "lost" by the jet becomes a wake in the plasma which then in turn becomes soft particles spread over a large range of angles around the jet direction, carrying momentum in the jet direction. Some of this momentum gets reconstructed as a part of the jet, meaning that this contributes to jet broadening unless soft particles are groomed away [28,35,38,49,[51][52][53][54][55]. The second effect arises from the interplay between the fact that higher energy jets are less numerous than lower energy jets and the tendency for narrow jets to lose less energy than wide jets. (This tendency is seen at weak coupling [56,57], in holographic models for jets at strong coupling [58], and in the hybrid model [28].) As a consequence, the jets that remain in any given energy bin after an ensemble of jets passes through a droplet of QGP tend to be narrower than the jets in that energy bin would have been absent the QGP: wider jets are pushed into lower energy bins, where they are much less numerous than the narrower jets found there [28,57,59,60]. So, even though individual jets may broaden, at the ensemble level there is a strong tendency for the jets with a given energy to be narrower after passage through the plasma than jets with that energy would have been. Before an experimental measurement of transverse momentum broadening can be made, careful work must be done to find ways to evade, or precisely measure, both of these confounding effects. Relative to our goals in this paper, this is a cautionary tale. Although what we are looking for (jets sprouting an extra prong due to a parton within the jet scattering at a large angle) sounds more distinctive, because such events will be rare the effort will require high statistics, judicious choice of observables, and a very considerable phenomenological modeling effort. Our results provide an initial input for such an effort.
The probability for picking up a given transverse momentum q ⊥ via a single hard scattering off a parton in the plasma was calculated previously [10,61], but only in the limit of infinite parton energy which means zero angle scattering. That is, these authors calculated the probability that an infinite energy parton picks up some significant transverse momentum q ⊥ in a Molière scattering, without changing its direction. Since what is most relevant to any experimental observable is the scattering angle, it is hard to use these results per se to gain guidance for what to expect in future experimental measurements. Here, we remedy this by providing a leading order perturbative QCD calculation of the Molière scattering probability for incident partons with finite energy, computing the probability distribution for both the scattering angle and the energy of the outgoing parton.
The computation of F (p, θ) in weakly coupled QGP, even a static brick of weakly coupled QGP, is a multiscale problem and, in addition, there are different phase space regions where F (p, θ) is governed by different processes, as discussed schematically in Ref. [14]. We specifically focus here on the kinematic regime in which the angle θ is sufficiently large that the dominant process is a single binary collision between the incident hard parton and a medium parton (a scatterer in the medium). For sufficiently large θ, the contribution from multiple scattering is not relevant since one single collision is more likely to give a large angle than multiple softer collisions in sum. At smaller values of θ, multiple softer collisions do add up and dominate, yielding a Gaussian distribution in the momentum transfer as discussed above. We shall focus on the large θ regime which is more likely to be populated via a single Molière scattering incident parton + target medium parton → outgoing parton + X . (1. 2) The second important way in which our calculation extends what has been done before is that we include all relevant channels. The parton that is scattered by a large angle need not be the incident parton, as in Rutherford scattering or deep inelastic scattering; it could be a parton from the medium that received a kick from the incident parton. We include this channel as well, and we shall see that in some kinematic regimes it is dominant. That is, in Eq. (1.2) the outgoing hard parton (the one that we imagine detecting via its contribution to some jet substructure observable or, if the incident parton represents an entire jet, via its contribution to an acoplanarity), as well as the X which goes undetected in our thought experiment, can each be either the deflected incident parton or the recoiling parton from the medium that received a kick. F (p, θ) describes the energy and momentum transfer of the incident parton to the medium and contains information about the nature of the scatterers in QGP.
In this work, we shall evaluate F (p, θ) for sufficiently large θ by following the standard methods of perturbative QCD. We then determine the probability distribution P (θ) for the angle of an outgoing hard parton by integration over p: Finally, we integrate P (θ) over θ to obtain N hard (θ min ), see Eq. (1.1). Our calculation allows us to estimate how rare large angle scatterings with some specified θ are and in this way can be used to provide qualitative guidance for the ongoing experimental search for evidence of point-like scatterers in QGP. This paper is organized as follows. In Section 2, we derive the expressions which relate F (p, θ) to a summation over all possible 2 ↔ 2 scattering process and obtain a compact expression involving the phase-space integration over the scattering amplitudes weighted by the appropriate thermal distribution function. We then describe how to sum over the individual processes as well as how to simplify the phase-space integration. The reader only interested in results, not in their derivation, can jump to Sec. 3, where we present our results and compare them to previous studies, including the computations done in the p in → ∞ limit in Refs. [9,10]. By considering incident partons with finite energy and including all relevant channels, our goal is to provide a quantitative tool for incorporation in future jet Monte Carlo calculations as well as qualitative guidance for how to use future precise, high statistics, suitably differential measurements of jet modification in heavy ion collisions to find the scatterers within the QGP liquid.

Kinetic Theory Set-up and Calculation Details
In this Section, we explain how we derive the probability distribution F (p, θ) for finding an outgoing parton with energy p at an angle θ relative to the direction of the incident parton.
Our key ingredient is the phase-space distribution f a (p, t) f a (p, t) ≡ Probability of finding an energetic parton of species a in a phase-space cell with momentum p at the time t, averaged over helicity and color states, where a can be u,ū, d,d, s,s or g. As emphasized in the definition, we neglect the dependence on helicity and color configurations. Although the phase-space distribution in principle can depend also on these variables, we assume that the medium is unpolarized and has no net color charge. Furthermore, if we average over the possible helicity and color configurations for the incoming hard probe, we are allowed to use the averaged distribution introduced in Eq. (2.1). We shall set our calculation up as a calculation of the time evolution of f a (p, t) in kinetic theory in which this distribution initially has delta-function support, describing the incident hard parton, and later describes the probability of finding an energetic parton of species a that has ended up with momentum p after a binary collision.

Initial conditions
We imagine a static brick of quark-gluon plasma, and we then imagine shooting an energetic parton with energy p in and momentum p in at it. The on-shell condition reads p 2 in = p 2 in , therefore p in denotes both the energy and the magnitude of the momentum for the incoming parton. (We shall assume that this parton does not radiate, split or shower during the time ∆t that it is traversing our brick of plasma, since our goal is to focus on large-angle scattering caused by a single binary collision. In future phenomenological studies in which our results are used within a jet Monte Carlo, results from our calculation would be used ∆t by ∆t by ∆t, with the value of ∆t chosen small enough that radiation or splitting is negligible during a single ∆t.) If the energetic parton of species a enters the medium at the initial time t I , the initial condition for the phase space distribution function reads where V is a unit volume that will not appear in any results. Here, we have fixed the initial energy and direction. Without any loss of generality we can take the z-axis to lie along the direction of the incident parton, which fixes cos θ in = 1. We normalize the expression in Eq. (2.2) in such a way that the incoming flux is one incoming parton per unit volume. The degeneracy factor ν a is defined as accounting for helicity and color configurations, with N c the number of colors. And, for later convenience we have introduced the definition of a function f I (p), where I refers to initial and is not an index, that describes the species-independent momentum-distribution in the initial condition.

Evolution of the phase-space distribution
We wish to answer the following question: if an incoming parton enters the medium at the time t I , what is the probability of finding an energetic parton of species a (not necessarily the same as that of the incident parton) exiting on the other side with a given energy and at a given scattering angle? In order to give a quantitative answer, we need to track the evolution of the function f a (p, t). At time t = t I , f a is zero for all p other than p = p in ; at later times, because the incident parton can scatter off partons in the medium f a can be nonzero at other values of p, and in particular at nonzero angles θ. Henceforth, we shall evaluate f a (p, t) at some nonzero angle θ, meaning that a labels the species of the energetic parton detected there. The calculation of the time evolution of f a (p, t) is performed in Appendix A, we report only the final result here. We assume that the probe scatters off a constituent of the medium at most once during its propagation through the medium over a time ∆t. We will later come back to this approximation and check when it is legitimate, namely when ∆t is sufficiently small and/or when θ is sufficiently large so that no summation over multiple scattering is needed. Within this approximation, the phase space distribution at the time t I + ∆t when the parton exits the medium takes the form f a (p, t I + ∆t) = ∆t ν a processes (2.4) The form of this expression can be readily understood for all scattering processes except qq ↔ gg or qq ↔ q q , where q and q are different flavors, as follows (although it applies to those processes too). Our convention is that the parton a detected in the final state comes from parton c in the initial state, and the undetected parton b comes from parton d. So, the n c f d term in the result (2.4) corresponds to the case where the outgoing hard parton a that is detected came from the medium, having been kicked out of the medium by the incident parton d, whereas the f c n d term corresponds to the case where the detected parton a came from the incident parton c, which scattered off parton d from the medium.
The [1 ± n B ] factor (where the sign is + if b is a boson and − if b is a fermion) describes Bose enhancement or Pauli blocking and depends on the occupation of the mode in which the undetected particle of species b in the final state is produced. The sum runs over all possible binary processes ab ↔ cd, with p , k (p, k) the momenta of c, d (a, b). The phase space integral is written in a compact form (2.5) The squared matrix elements are summed over initial and final helicity and color configurations, without any average. The term with the Kronecker delta function accounts for the cases when c and d are identical particles. Finally, we must specify the "soft" medium distribution functions n a (p). As we discussed in Section 1, we shall choose to use distributions as if the quarks and gluons seen in the QGP by a high-momentum probe were massless, noninteracting, and in thermal equilibrium, meaning that n a (p) depends only on the statistics and energy of the particle in the medium that is struck and is given by Note that we are considering a medium in which the chemical potential for baryon number vanishes, meaning that the equilibrium distributions for quarks and antiquarks are identical. For this locally isotropic medium, the equilibrium distributions depend on the parton energy p but not on the direction of its momentum. They are also time-independent, since we are considering a static brick of plasma with a constant T . By taking a noninteracting gas of massless quarks, antiquarks and gluons, in thermal equilibrium, as our medium we are n Process Table 1. List of the binary collision processes that can produce a hard parton in the final state with large transverse momentum with respect to the incoming probe. Here, q and q are quarks of distinct flavors,q andq the associated antiquarks, and g is a gauge boson (gluon). The third column lists explicit leading order expressions for the corresponding QCD squared matrix elements, in vacuum, summed over initial and final polarizations and colors, as a function of the standard  defining a benchmark, not an expectation. As we noted in Section 1, we look forward to the day when comparisons between experimental data and predictions made using our results incorporated within a jet Monte Carlo are being used to determine how n a (p) for QGP differs from the benchmark that we have employed here. A future program along these lines could be thought of as the analogue, for a thermal medium, of determining the parton distribution functions for a proton.
Initially, at time t I , f a takes on the form (2.2) and is zero for all p except for p = p in . The expression (2.4) encodes the fact that after the incident parton has propagated through the medium for a time ∆t, because there is some nonzero probability that a 2 → 2 scattering event occurred there is now some nonzero probability of finding a parton with any p.

QCD matrix elements
The formalism set up so far is valid for a generic theory with arbitrary degrees of freedom and arbitrary interactions giving rise to binary scattering processes, and relies principally just on the kinematics of the binary collisions. The specific dynamics becomes relevant only when we have to specify the matrix elements in Eq. (2.4). We do so here, in so doing specializing to QCD. We collect the results for the matrix elements for all processes relevant to our study in Table 1. We label each process with an integer index (n = 1, 2, . . . , 11), and we write the associated matrix element summing over initial and final colors and polarizations. We also assign to each process a degeneracy factor w (n) , different for each degree of freedom involved in the collision, which will be useful shortly. With these matrix elements in hand, we can evolve the initial phase-space distribution given in Eq. (2.2) by plugging it into Eq. (2.4). In this way, we obtain the phase-space probability after the incident parton has spent a time ∆t in the medium.
In addition to neglecting all medium-effects in the distribution functions (2.6) as we discussed in Section 2.2, we shall do the same in the QCD matrix elements for 2 → 2 collisions. This means that we are assuming weak coupling throughout and furthermore means that we can only trust our results in the kinematic regime in which the energy and momentum transferred between the incident parton and the parton from the medium off which it scatters is much larger than the Debye mass. We shall check this criterion quantitatively in Section 3.3.

Probability distribution after passage through the medium
Having derived the evolution of the phase-space distribution in Eq. (2.4), we can now define and compute the probability distribution, which is the main result of this paper. Thus far, we have denoted different parton species with lower case letters (i.e. a = u,ū, d,d, s,s, g). It is convenient to introduce uppercase indices denoting different types of partons: gluons, quark and antiquarks (i.e. A = G, Q,Q). We use this notation to define the probability distribution that we introduced in Fig. 1: F C→A (p, θ; p in ) ≡ Probability of finding a parton of type A with energy p at an angle θ with respect to the direction of an incoming parton of type C with energy p in .
This quantity is given by the sum over all possible processes with C and A in the initial and final state, respectively. Its explicit expression reads The prefactor in front of the sum is the Jacobian of the phase-space integration The sum runs over all the lowercase indices corresponding to parton species of the type A. For example, if A stands for a quark, the sum runs over the values a = u, d, s. The degeneracy factor ν a appears because our distribution functions are averaged over colors and polarizations; the detector cannot resolve these quantum numbers, we account for all of them by this multiplicative factor. Finally, we note that the distribution function f a (p, θ; t I + ∆t) appearing in Eq. (2.8) is the time-evolved quantity given in Eq. (2.4), evolved from an initial condition at time t I given by for all other values of a (2.10) where the function f I (p in ) was defined in Eq. (2.2). (For example, if C = Q meaning that the incident parton is a quark then f a is nonzero for either a = u or a = d or a = s, and the flavor of the incident quark makes no difference to our calculation.) We have defined the probability (2.7) such that it does not distinguish between quarks of different flavors, but it does distinguish between quarks, antiquarks and gluons. So, if our goal is to find the total probability of finding any energetic parton in the final state with energy p and angle θ, we have to sum over the different types of partons. As an example, if we consider an incoming quark, the probability of getting any energetic parton in the final state reads (2.11) In the last step in our derivation, we directly plug the expression for the time-evolved phase-space distribution given in Eq. (2.4) into our expression for the probability distribution (2.8). Before doing that, it is useful to introduce some notation to make our final expression more compact. We define the generalized Kronecker delta functionsδ a,G ≡ δ a,g , δ a,Q which equals 1 if a = u or d or s and which vanishes for other values of a, andδ a,Q which equals 1 if and only if a =ū ord ors. Moreover, we define the generalized medium "soft" distribution function where n B.E. (p) and n F.D. (p) are the Bose-Einstein and Fermi-Dirac distributions from Eq. (2.6), respectively. With this notation in hand, we can now write the complete leading order expression for the probability function defined in Eq. (2.7): (2.13) Here, we have defined a dimensionless parameter κ multiplying the overall expression via (2.14) κ becomes large either for a thick brick (large T ∆t) or for a large value of the QCD coupling constant g s that controls the magnitude of all the matrix elements for binary collision processes. Note that the V in the prefactor of Eq. (2.13) cancels the 1/V from Eq. (2.2), meaning that no V will appear in any of our results. Henceforth we shall not write the factors of V . Note also that neglecting multiple scattering as we do is only valid when N hard , the integral over F C→A (p, θ; p in ) defined in Eq. (1.1), is small. For any given choice of p and θ, if κ is too large multiple scattering cannot be neglected and our formalism breaks down. Equivalently, for any given κ our formalism will be valid in the regime of p and θ, in particular for large enough θ, where F C→A (p, θ; p in ) is small and multiple scattering can be neglected. The sum over n in Eq. (2.13) runs over all the 11 processes in Table 1. The deltaδ a,A ensures that only processes with a parton of type A present in the final state are accounted for. Crucially, each process is multiplied by the C-dependent weight factor w (n) C , given explicitly in the last three columns of Table 1. As an example, if we are considering the production of A = Q from an incident gluon, C = G, via gg →qq, the weight factor w (8) G is N f since we can produce this final state by pair-production of any flavor of light quark. Thus, this multiplicative factor accounts for the multiple ways a given process can produce the energetic parton A in the final state. When such an outgoing parton originates from an incident parton c, the matrix element has to be multiplied by the thermal weight δ c,C f I (p ) n d (k ), whereas when the incoming parton is d this factor isδ d,C f I (k ) n c (p ).
The expression Eq. (2.13) is the central result of this paper, albeit written in a compact and hence relatively formal fashion. We note again that this relation is valid only as long as ∆t is much shorter than the characteristic time between those binary collisions between the incident parton and constituents of the medium that produce scattered partons with a given p and θ. We will see in Section 3 that this is true as long as the scattering angle is larger than some θ min , where θ min will depend on p, p in and κ. Before turning to results in Section 3, in Section 2.5 we shall write the expression (2.13) more explicitly in specific cases and in Section 2.6 we shall describe some of details behind the computations via which we obtain our results.

How to sum over different processes
In order to write the expression (2.13) more explicitly and in particular in order to sum the various different phase space integrals over various different matrix elements that contribute to a given physical process of interest, it is convenient to define the following set of phase space integrals: where the index n spans the 11 different binary collision processes listed in Table 1. The ± sign in both equations correspond to the cases where B is a boson or a fermion, respectively. For processes with identical incoming partons (and also for process 8 in Table 1), we have If we look back at Eq. (2.13), we notice that we can always express F C→A (p, θ; p in ) as a weighted sum over (n) D,B and ( n) D,B . Obtaining such expressions is the goal of this Section. There are 3 × 3 = 9 different cases, corresponding to three options for both the incoming and outgoing parton: quark, antiquark or gluon. We shall first list 4 cases, corresponding to choosing either quark or gluon. Replacing quarks by antiquarks gives 3 more cases, with identical results. We shall end with the 2 cases where the incoming and outgoing partons are quark and antiquark or vice versa. The brick of quark-gluon plasma is assumed to not carry a net baryon number, therefore the results for these last 2 cases are also identical. In the remainder of this subsection, we give explicit expressions for these 5 independent results. For each case, we define the partial contributions as follows That is, we decompose the total probability that we are interested in into a sum of up to 11 different terms, one for each of the processes listed in Table 1. As we will see shortly, only a subset of them will actually contribute in each case. For example, in order to understand which ones are relevant to F Q→Q (p, θ; p in ) we need to look at Table 1 and identify those processes with at least one quark in the initial and in the final states. The final result for each case can then be expressed in terms of the functions defined in Eqs. (2.15a) and (2.15b). Individual processes in Table 1 can contribute in more than one case; for example, process 9, quark-gluon scattering, contributes to the probabilities for four cases: F Q→Q (p, θ; p in ) ("incident quark, outgoing quark"): We start from the case where both the incoming and the outgoing parton are quarks. The relevant processes are the ones labeled by n = 1, 3, 4, 6, 7, 9 in Table 1 with individual expressions given as follows. First, where the factor 1/2 is a symmetry factor (see Eq. (2.13)), and w Q is read from Table. 1. In the last step, we have used the fact that (1) Q,Q = 1 Q,Q according to the relation (2.16). Likewise, since the squared matrix elements for the processes 4 and 6 are identical. And, Upon summing the above, we find the final result F Q→G (p, θ; p in ) ("incident quark, outgoing gluon"): This case gets contributions from the processes labeled by n = 8, 9. We identify again the individual contributions to the total probability where we have used the relation (2.16). And, which add up to give the final result for this case F G→Q (p, θ; p in ) ("incident gluon, outgoing quark"): The calculation for this case is analogous to the previous one. The partial contributions read which, after summing, result in F G→G (p, θ; p in ) ("incident gluon, outgoing gluon"): When both the incoming and outgoing energetic partons are gluons, the processes contributing to the probability distribution are the ones labeled by n = 9, 10, 11. The individual terms are where we have take into account the fact that processes 9 and 10 have identical squared matrix elements. And, where once again we have used the relation (2.16). Consequently, we find The last case we consider is when a quark enters the medium and an energetic antiquark exits on the opposite side. The processes that contribute to this case are where we use the fact that processes 6 and 4 have identical squared matrix elements. In addition, where we have use the relation (2.16). The total probability for this case is (2.27)

Phase space integration
After performing the summation over different processes, our final task is to evaluate the phase space integrals in Eqs. (2.15a) and (2.15b). The expression in Eq. (2.15a) involves a 9-fold integration in the phase space (p , k , p). We first integrate over a 4-dimensional delta function in Eq. (2.5). The integration over the azimuthal angle is straightforward. Finally, we perform two more integrations by taking the advantage of the delta function in f I . (See Appendix. B.1 for details.) Upon following techniques widely used in the literature (see e.g. Refs [63][64][65]), we find (2.28) Here, k T denotes the energy of the thermal parton from the medium whose momentum we shall denote by k T and k X = k + ω denotes the energy of the undetected final state parton. The integration range starts from the value corresponding to the minimum energy allowed by kinematics for the thermal parton from the medium. Moreover, φ is the angle between the two planes identified by the pair of vectors (p, q) and (q, k T ), and we use ω ≡ p in − p and q = p − p in to denote energy and momentum difference between the incident parton and the outgoing parton that is detected. The matrix elements M (n) that appear in Eq. (2.28) are to be taken from Table 1, with the Mandelstam variables t and u occurring within them specified in terms of quantitiest andũ that can be expressed as functions of q, ω, k T and φ as follows where in the matrix elements in Eq. (2.28) we have simply t =t and u =ũ but where we will need to set t =ũ and u =t below in our result for ( n) D,B . Here, q, andt can be expressed as functions of p, p in and cos θ thus: Following a calculation that proceeds along similar lines, the quantity in Eq. (2.15b) can be expressed as where the role oft,ũ are interchanged in the squared matrix element with respect to Eq. (2.28). There are two integrations left in Eqs. (2.28) and (2.33), over φ and k T . Remarkably, the integration over φ can be performed analytically, as explained in Appendix B.2. The remaining integration over k T has to be performed numerically.

Results and discussion
The purpose of this work is to evaluate F C→A (p, θ), the probability distribution for finding an outgoing hard parton of type A with energy p and angle θ relative to the direction of an incident hard parton of type C with energy p in . (For simplicity, here as in the Introduction we shall write F C→A (p, θ; p in ) as just F C→A (p, θ).) Recall that by "type" we mean gluon or quark or antiquark. We consider a static brick of a weakly interacting QGP, and have included the contributions from a single binary collision between the incident hard parton and a medium parton. In Section 2, we have presented a careful derivation of the expression for F C→A (p, θ) in Eq. (2.13), and have provided further technical details on the summation over different processes in Section 2.5, as well as the simplification of the phase space integration in Section 2.6. By summing over different types, we obtain the probability distribution for finding final parton of any type, Integration of F C→all (p, θ) over p using Eq. (1.3) then yields P (θ), namely the probability distribution for the angle θ.

Comparison with previous work
Before we present our results, we shall briefly sketch how they agree with results obtained previously where they should. The details of this comparison are found in Appendix C. The probability distribution for an energetic parton that travels for a distance L through a weakly coupled QGP to pick up transverse momentum q ⊥ , which we shall denote P(q ⊥ ), was analyzed in Ref. [10]. These authors confirmed that for sufficiently small L or for sufficiently large q ⊥ , P(q ⊥ ) will approach P single (q ⊥ ) (denoted by P thin (q ⊥ ) in Ref. [10]), the probability distribution obtained upon including at most a single scattering between the incident parton and a scatterer from the thermal medium. This is expected on physical grounds since the most probable way of picking up a large q ⊥ is via a single scattering. Expressions for P single (q ⊥ ) were calculated previously under the condition q ⊥ T in Ref. [66] and under the condition q ⊥ T in Ref. [61]. The calculations of Ref. [10] do not assume any ordering between q ⊥ and T , and their results agree with the older results in the appropriate limits. In all of these previous studies, however, the calculations are performed by first taking a limit in which p in /T → ∞ while q ⊥ /T remains finite, meaning in a limit in which θ → 0. In this limit, Rutherford-like scattering in which an incident parton scatters off a parton from the thermal medium is dominant over all other 2 ↔ 2 processes, including those in which a parton from the medium is kicked to a large angle as well as processes such as qq ↔ gg. We shall not take the p in /T → ∞ limit, meaning that we must include all 2 ↔ 2 processes and that we can describe scattering processes that produce a parton at some nonzero angle θ and hence can compute P (θ), the probability distribution for the scattering angle θ.
To compare to the previous results referred to above, we take the limit θ 1 in our result for P (θ) and compare what we find there with P single (q ⊥ ) from Refs. [10,61,66].
Consequently, to compare to previous results we evaluate F C→all (p, θ) in the regime and then perform the necessary integrations to obtain P (θ) in this regime. In Eq. (C.5) in Appendix. C.1, we show that our results agree with those from the literature if P (θ) is given by . In subsequent parts of Appendix C, we confirm in detail that our results do indeed match those found in Refs. [10,61,66] in the kinematic regime where they should.
In this work, we have extended the previous studies by considering finite (but large) p in /T meaning that ω/p in and θ need not vanish. Consequently, there are new features in our computations. In particular, we have included all 2 ↔ 2 scattering processes, as given in Table 1, in our evaluation of F C→A (p, θ). Furthermore, when ω/p in is finite, either the deflected incident parton or the recoiling thermal parton or both can show up with energy p and angle θ. Indeed, we shall see in the subsequent sections that P (θ) at nonzero θ differs qualitatively from that obtained by extrapolating its behavior in the small θ limit. In particular, the large-angle tail of P (θ) is in reality fatter than one would guess from such an extrapolation. This makes the inclusion of all 2 ↔ 2 processes as we do important and interesting, not just necessary.
Next, we note that by working at finite p in /T we introduce a kinematic cutoff on the momentum transfer, meaning that when we increase θ the probability distribution P (θ) must eventually be suppressed since (because of energy/momentum conservation) the minimum energy of the thermal parton needed to yield a specified θ will become much larger than T . We shall illustrate this quantitatively later, see the blue curves in Fig. 5. The analogous kinematic cutoff on q ⊥ in P single (q ⊥ ) computed in the limit in which p in /T → ∞ and θ → 0 is less constraining [10].
Finally, we note that in Ref. [38] quantities analogous to F (p, θ) or integrals of F (p, θ) have been computed in the Linear Boltzmann Transport (LBT) model for energetic partons shooting through a brick of weakly coupled QGP as in our calculation, albeit largely with a focus on a kinematic regime in which p, and hence the momentum transfer, are only a few GeV. These authors also compute a quantity directly related to the transverse momentum distribution P (q ⊥ ) using the LBT model for q ⊥ out to around 10 GeV, and provide a very interesting study of how the distribution becomes more and more Gaussian as the thickness of the brick is increased. However, even for the thinnest brick that they consider the values of q ⊥ that they investigate are not large enough for single scattering to be dominant. It would be interesting to extend these LBT calculations to larger q ⊥ where the probability of multiple scattering is negligible and compare them to our results, upon taking into account the appropriate Jacobian.

Results for the probability distributions F (p, θ) and P (θ)
We shall now present results from our numerical calculation of F C→all (p, θ) /κ as well as for P (θ)/κ, both of which are independent of κ. Recall that κ ≡ g 4 s T ∆t. The probability for a single 2 → 2 scattering with any specified kinematics is proportional to g 4 s at tree-level, and is proportional to ∆t ≡ L/c, the time that the incident parton would spend traversing the brick if it did not scatter. Hence, increasing κ (either via increasing the coupling or via increasing T ∆t) must increase F C→all (p, θ) and P (θ). Upon increasing κ, though, at some point the assumption that single scattering dominates must break down, and along with it our calculation. The criterion here is that N hard (θ min ), defined in Eq. (1.1), must remain small and this defines an upper limit on the value of κ at which our calculation can be used for angles θ greater than any specified θ min , or a lower limit on the angle θ at which our calculation can be used for any given value of κ. We shall illustrate this quantitatively in Section 3.4. Note that in this Section we shall work in the weak coupling limit g s → 0 in which κ → 0 and our expression for F C→A (p, θ) in Eq. (2.13) is valid for any nonzero θ and any finite ∆t.
We shall consider N f = 3 throughout and we shall only consider QGP with no net baryon number, meaning zero baryon number chemical potential and meaning that the distribution of quarks in our thermal medium is the same as that of antiquarks.
We begin our discussion by considering an incident gluon with p in /T = 100. In the top row of Fig. 2 that the probability distribution is peaked at p ≈ p in , meaning that outgoing partons with a very small angle are likely to have a small value of ω/p in , where ω = p in − p. This implies that computing F G→all (p, θ) in the limit (3.3) is sufficient to obtain P (θ) for θ 1, as we mentioned earlier. However, the dependence of F G→all (p, θ) on p changes qualitatively as we increase θ. F (p, θ) at θ = 0.4 and θ = 0.8 are both largest at small values of p and decrease monotonically with increasing p. To understand this, let us recall that the difference between p and p in , i.e. ω, measures the energy transfer during a binary collision, with a smaller p corresponding to a larger energy transfer ω. Likewise, a larger θ means a larger transverse momentum transfer. Since the typical energy of a thermal parton is quite soft, of order T , a large momentum transfer in a single collision between an incident parton and the thermal scatterer is more likely to be accompanied by a large energy transfer. That is why we see F G→all (p, θ) telling us that when we ask about scattering at large θ we find that it most often corresponds to scattering with a large ω and hence a small p. Equivalently, although in different words, we note that in this regime the detected parton is most likely to be a parton from the medium that was kicked to a large angle θ by the incident parton, with the incident parton having lost only a small fraction of its energy to the parton that is detected. The energy transfer defined as ω is large because the detected parton is the parton from the medium, not the incident parton.
In Fig. 2, in addition to plotting F G→all (p, θ) we have also shown its separate components corresponding to detecting an outgoing gluon or an outgoing quark or antiquark, namely F G→G (p, θ) and F G→Q (p, θ)+F G→Q (p, θ). (Note that F G→Q (p, θ) = F G→Q (p, θ).) While F G→Q (p, θ) F G→G (p, θ) at small θ, meaning that at small θ the outgoing parton is most likely to be a gluon when the incident parton is a gluon, we see that F G→Q (p, θ) + F G→Q (p, θ) eventually becomes comparable to F G→G (p, θ) at larger values of θ. This confirms that what is being seen at large values of θ and small values of p is to a significant extent partons from the medium that have been struck by the incident parton. The quarks and antiquarks seen in this regime also include those coming from the process gg → qq. And, this observation convincingly demonstrates that Rutherford-like scattering is not dominant over other processes at larger values of θ We now consider an incident gluon with a lower initial energy, i.e. p in /T = 25, and plot F G→all (p, θ)/κ for this case in the second row of Fig. 2. As before, we have selected three representative values for θ, from left to right choosing θ = 0.1, 0.4 and 0.8. The behavior of F G→all (p, θ) as a function of p is qualitatively similar to that with p in /T = 100: F G→all (p, θ) features a peak at p ≈ p in at small θ, but it then becomes a decreasing function of p/T at the larger values of θ. At a quantitative level, we observe that for θ = 0.1, the peak value of F G→all (p, θ) with p in /T = 25 is much larger than that with p in /T = 100. This is due to the dominance of Rutherford-like scattering at small θ, since the probability of Rutherford scattering decreases with increasing q ⊥ ≈ p in θ and we are comparing two values of p in at the same small θ. As with p in /T = 100, we see that when we choose θ = 0.8 we find a probability that is peaked at small p and we see that the contribution of quarks and antiquarks is not much smaller than that of gluons. Hence, at this large value of θ we are seeing partons kicked out of the medium. We see that with p in /T = 25 the choice of θ = 0.4 represents an intermediate case.
For completeness, in the third and fourth rows of Fig. 2 we plot F Q→all (p, θ) for an incident quark with p in /T = 100 (third row) and 25 (fourth row) at three values of θ. We have multiplied our results for an incident quark by the ratio of Casimirs C A /C F , which is 9/4 for N c = 3, to simplify the comparison to our results for an incident gluon. After taking this Casimir scaling factor into account, the resulting F Q→all (p, θ) are very similar  to those for incident gluons with the same choice of p in /T . Similar to what we found for gluons, if we look at small θ and p close to p in , we see that the Rutherford-like Q → Q process makes the dominant contribution whereas if we look at larger θ and small p we see that Q → G is comparable to, and in fact slightly larger than, Q → Q. This demonstrates that Rutherford-like scattering is not dominant here and suggests that the detected parton is most often a parton that was kicked out of the medium. To complement Fig. 2, which illustrates the dependence of F G→all (p, θ) on p with fixed θ, in the top row of Fig. 3 we show the dependence of F G→all (p, θ) on θ at three fixed values of p/T . In another words, in Fig. 3 we are looking into the angular distribution of an outgoing parton with a fixed p/T , considering three different values of p/T , namely 80, 40 and 20. We have chosen an incoming gluon with p in /T = 100 in all three panels. In the second row of Fig. 3, we show results for an incoming quark with the same p in /T . As before, we see that after, after multiplying by the ratio of Casimirs 9/4, F Q→all (p, θ) is reasonably similar to F G→all (p, θ). From our results with p/T = 80, we see that when we look at outgoing partons whose energies are not much lower than those of the incident parton, smaller values of the scattering angle θ are favored and the scattered parton is dominantly the same type as the incident parton. In contrast, in our results at smaller p/T we see a much broader θ distribution and, in particular at larger values of θ, we see comparable contributions from quarks or antiquarks and gluons in the final state, confirming that the detected parton was a parton from the medium that was struck by the incident parton.
We now present our results for the probability distribution P (θ), which we obtain by integrating F C→all (p, θ) over p, following Eq.  In addition to plotting the probability distribution for finding any outgoing parton at a given θ as the red curves, we also present its breakdown into the cases of an outgoing gluon (blue curves) and an outgoing quark or antiquark (orange curves). In the right column, we plot P (θ) for an incident gluon (red) as well as for an incident quark times 9/4 (black dashed curves) as well as the θ 1 result P AD (θ) from Eq. (3.4) first obtained by Arnold and Dogan [61].
somewhat arbitrary choice of p min /T , we will consider two different choices, p min /T = 10 and p min /T = 20, and check the sensitivity of P (θ) to this variation in this choice. We observe that for sufficiently small θ, P (θ) is insensitive to the choice of p min /T . This is to be expected, given our discussion of F C→all (p, θ): recall that it is peaked at p ∼ p in p min for small θ, meaning that where we place p min does not matter much in this case. However, when we choose a larger value of θ the magnitude of P (θ) becomes much smaller if we increase p min /T from 10 to 20. This is also expected since at large θ we have seen that F (θ, p) is a rapidly decreasing function of p. In the bottom-left panel of the figure, we see similar behavior in the case in which the incident gluon has p in /T = 25. When θ is not small, P (θ) is highly suppressed when we choose p min /T = 20. This is no surprise since for this choice p min is close to p in , meaning that the phase space included in the integration (1.3) is quite restricted. In both the panels in the left column of Fig. 4, we have in addition plotted P (θ) for an outgoing gluon, G → G, and for an outgoing quark or antiquark, G → Q. At small angles Rutherford-like scattering dominates and since the incident parton is a gluon we see that the probability to find an outgoing gluon is much greater than that for an outgoing quark or antiquark. At larger angles Rutherford-like scattering is no longer dominant, the parton that is detected most likely comes from the medium, and we see that the probability to find an outgoing quark or antiquark becomes comparable to the probability to find an outgoing gluon.
In the right panels of Fig. 4, we compare P (θ) for an incident gluon with that for an incident quark with the same choice of p in /T and p min /T multiplied by C A /C F . We observe that, after taking into account the appropriate Casimir scaling factor, P (θ) is almost identical for both cases.
As we discussed in Section 3.1, the transverse momentum distribution due to a single binary scattering P single (q ⊥ ) has been obtained previously in the small θ limit (3.3) [10,61]. If in addition q ⊥ T , P single (q ⊥ ) reduces to the expression first derived by Arnold and Dogan (AD) in Ref. [61] which we shall denote P AD single (q ⊥ ) and which we provide explicitly in Eq. (C.7). (See also Ref. [10]).) In the small θ limit, we can convert P AD single (q ⊥ ) to a probability distribution for the angle θ that we shall denote by P AD (θ) using the Jacobian (C.3). We obtain where ζ(3) ≈ 1.202 is the Riemann zeta function. Here, the incident parton is a gluon; for the case of an incident quark, one has to replace C A with C F in Eq. (3.4). In the two panels in the right column of Fig. 4, we have compared P (θ) with P AD (θ) extrapolated to finite θ. We observe that, as expected, P AD (θ) agrees very well with P (θ) at small θ. However, the large-angle tail of P (θ) is much fatter than that of P AD (θ) when p in /T = 100 for all p min /T under consideration, as well as when p in /T = 25 for p min /T = 10. This implies that when p in p min , it is important to include all 2 → 2 scattering processes as we have done, not only the Rutherford-like scattering process that dominates at small θ.
The results that we have illustrated in this Section are the principal results of our calculation. We have presented them here upon dividing F (p, θ) and P (θ) by κ ≡ g 4 s T ∆t. This is the appropriate form in which to provide them to anyone incorporating them in a future jet Monte Carlo calculation, since the values of the coupling g s and the time-step ∆t will be provided by that calculation and in such a calculation the local value of T will come from the description of the expanding cooling droplet of QGP which the Monte Carlo jet is traversing. As described in the Introduction, we also wish to provide some qualitative guidance for the planning of future experiments and for how to use future precise, high statistics, suitably differential measurements of jet substructure modification in heavy ion collisions to find the scatterers within the QGP liquid. To this end, in Section 3.4 we shall illustrate our results for P (θ) and its integral N hard (θ min ) using phenomenologically motivated values for various input parameters including κ. First, though, in the next Section we shall discuss the regime of validity of our calculation.

Regime of validity of the calculation
In this Section we pause to discuss the domain of applicability of the calculations presented in the previous Section. We have assumed that single scattering dominates, neglecting multiple scattering. This assumption is valid when N hard (θ min ) is much smaller than one, a criterion that depends on the value chosen for κ. We therefore leave the assessment of this criterion to Section 3.4, in particular to Fig. 6. We shall focus here on a different limitation of our calculation. Since we are neglecting all medium-effects in the QCD matrix elements for 2 ↔ 2 collisions, our results are trustable only in the kinematic regime in which the energy and momentum transferred between the incident parton and the parton from the medium off which it scatters are both much larger than the Debye mass m D . That is, our results are trustable only in the regime where Here, we will denote the square of the four momentum difference between the incident parton and the detected outgoing parton and that between the incident parton and the undetected parton byt andũ respectively, as in Section 2.6. By using Eq. (2.32), in which t is expressed in terms of p in , p and θ, we can determine the region in the (θ, p/T ) plane where the condition −t m 2 D is satisfied for any given p in and m D . Furthermore,ũ can be written as where p X and θ X are determined from transverse momentum conservation and energy conservation, respectively, and are given by where k ⊥ denotes the transverse momentum of the thermal scatterer. While in generalũ also depends on the magnitude of the momentum of the parton from the thermal medium k = |k|, we can expressũ in terms of p in , p, and θ for any value of θ that is not too small because the characteristic values of k ⊥ and k are of the order of T . First, since p T , the transverse momentum of the outgoing parton, p sin θ, will be much larger than T when θ is not too small. To balance such a large transverse momentum, we need to have p X sin θ X ≈ p sin θ. Second, we have observed from our study of F C→all (p, θ) in Section 3.2 that when the momentum transfer is large, the energy transfer in a binary collision is also likely to be large, i.e. ω T . We therefore have from energy conservation (3.7) that p X ≈ p in −p = ω. Combining the above two approximations and substituting into Eq. (3.6), we obtainũ  Fig. 4. The solid blue curves in both panels are determined by the condition k min = 7T where k min is the minimum possible value of the energy of a parton in the medium that, when struck by a parton with incident energy p in , can yield an outgoing parton at a given point in the (θ, p/T ) plane. k min is given by the expression (2.29), and we have used p in /T = 100 (left panel) or p in /T = 25 (right panel) in our numerical evaluation of k min . All our results become smaller and smaller farther and farther above the blue curves. Hence, our calculations are valid and our results are not small in the region below the blue curves and above the red and orange curves.
In Fig. 5, we illustrate the regimes in the (θ, p/T ) plane where the conditions (3.5) are satisfied for incident partons with p in /T = 100 and p in /T = 25. We use the standard expression for Debye mass squared: 9) choosing N c = N f = 3 and, as described in the next Section, choosing g s = 1.5. The red dashed and orange dotted curves are determined by solving −t = 10m 2 D and −ũ = 10m 2 D , respectively. We observe that the conditions (3.5) are satisfied for sufficiently large θ, although how large θ needs to be depends on the values of p in /T and p/T . The blue curves in Fig. 5 do not represent limits on the validity of our calculation. However, above the blue curves the results that we obtain must be small in magnitude, for the following reason. For scattering processes to yield outgoing partons with values of (θ, p/T ) above the blue curves, the only partons from the medium that can contribute are those with energies k greater than 7T , whose n a (k) in (2.6) are smaller than 10 −3 . For this reason, the probability for scattering events that yield outgoing partons above the blue curves must be small. Hence, the regime in the (θ, p/T ) plane where medium effects can be neglected in the matrix elements for 2 ↔ 2 scattering as we do and where our calculations yield a significant scattering probability is the region above the red and orange curves and below the blue curves.

Estimating P (θ) and N hard (θ min ) for phenomenologically motivated inputs
In Fig. 4 in Section 3.2, we have evaluated P (θ)/κ. By dividing the probability distribution P (θ) by κ we obtained and plotted κ-independent results. And, as we noted in Section 3.2, this is the form of our results that we should provide for use in a future jet Monte Carlo analysis, which is the path to phenomenologically relevant predictions for experimental observables. It may also be interesting to study the importance of processes in which a photon is radiated [67] as well as 2 → 3 scattering processes in future phenomenological studies. This is for the future. In the present paper, we would like to get at least a qualitative sense of P (θ) for incident partons with several values of p in . This means that we need to input phenomenologically motivated values of g s , ∆t, and T -and hence κ.
Since we are interested in those binary collisions with characteristic momentum transfer which is of the order 10 GeV, following Ref. [68] we will use g s = 1.5 as our benchmark value in the following analysis. Of course in reality g s runs, meaning that in a future calculation that goes beyond tree-level one should allow g s to depend on the momentum transfer in a particular collision. Working at tree level as we do, it is consistent just to pick a value of g s , and we shall choose g s = 1.5. We shall pick T = 0.4 GeV as the temperature of our brick of QGP and ∆t = 3 fm as the time that a parton spends in our brick of QGP. With these choices of parameters, κ ≈ 30. (The actual value is 30.84, but this would be misplaced precision. We shall use κ = 30 in plotting results in this Section.) While we should only expect our calculation to be quantitatively reliable for g s 1, we hope our results with g s = 1.5 will be of qualitative value in estimating the magnitude of P (θ) as well as its θ-dependence. (We also note that g s = 1.5 corresponds to α QCD ≈ 0.18, in many contexts a weak coupling.) Of course, any reader who has their own preferred values of g s , T and ∆t that they like to use to make phenomenologically motivated estimates should feel free to do so. Our result for P (θ) is simply proportional to κ = g 4 s T ∆t. We will concentrate on the case where the incident parton is a gluon. We plot P (θ) in the left column of Fig. 6 for p in /T = 25 (upper left) and 100 (middle left), in each case for p min /T = 10 and 20. These curves correspond to results shown in Fig. 4, multiplied by κ = 30. Taking T = 0.4 GeV, they correspond to incident gluons with p in = 10 GeV and 40 GeV and scattered partons with p > 4 GeV and 8 GeV. In the lower left panel, we plot P (θ) for p in /T = 250, corresponding to p in = 100 GeV, for scattered partons with p > 10 GeV and p > 40 GeV. As we have demonstrated in Fig. 4, P (θ) for an incident quark can be well described by multiplying P (θ) for an incident gluon by the ratio of Casimirs C F /C A = 4/9.
In the right column of Fig. 6, we integrate P (θ) over θ and obtain N hard (θ min ), defined in Eq. (1.1). (Since P (θ) drops very quickly for large values of θ, when we evaluate N hard (θ min ) numerically we stop the integration in Eq. (1.1) at θ = 1.5.) Among the quantities that we can calculate, N hard (θ min ) is perhaps the most useful for the purpose of obtaining a qualitative sense of how large the effects of point-like scatterers in the QGP will be. For example, reading from the dashed red curve in the middle-right panel of Fig. 6, we see that if an incident gluon with p in = 100T = 40 GeV traverses 3 fm of QGP with a temperature  of 0.4 GeV, the probability that a parton with an energy p > p min = 20T = 8 GeV is detected at some angle θ > 0.8 is around 1/1000, while this probability rises to around 1/100 for detection at an angle θ > 0.5, and the probability that a parton with p > p min = 10T = 4 GeV is detected at an angle θ > 0.8 is around 1/20. This gives a sense of the probability of kicking partons to these angles and as such is helpful in making qualitative assessment of how small (meaning how improbable) the effects that will need to be looked for via detecting suitable modifications to jet substructure observables may be. We would be happy to provide curves depicting our results for N hard (θ min ) or P (θ) for different choices of p in , p min , T , ∆t and g s . In the middle row of Fig. 6, where we consider incident partons with p in = 100 T , we have also included results where we only count scattered partons with p > p min = 40 T (the red dotted curves). This allows us to look at the dependence of our results on p in in two ways. If we compare the red solid curves above (p in = 25 T and p min = 10 T ) to the red dotted curves in the middle panels (p in = 100 T and p min = 40 T ) we see that increasing p in while increasing p min proportionally rapidly reduces the probability for large angle scattering. This corresponds to increasing the momentum transfer in the binary collision, and is qualitatively as one would expect based upon intuition from Rutherford scattering. On the other hand, if we compare the solid red curves in the top and middle panels, or the dashed red curves in the top and middle panels, we see that increasing p in while keeping p min fixed results in a much smaller change in the probability for large angle scattering. This corresponds to the observation that the probability for kicking a parton with p p min for some fixed p min out of the medium at some fixed large angle θ increases slowly with increasing p in . This further highlights the importance in our results of processes other than Rutherford scattering where what is detected is a parton that was kicked out of the medium.
In Fig. 6, we have only plotted our results (the red solid, red dashed and red dotted curves) for P (θ) and N hard (θ min ) at large enough values of θ and θ min that the the condition (3.5) is satisfied. As we discussed in Section 3.3, our calculation breaks down at smaller values of θ. For example, for p in /T = 100 and p min /T = 20 we observe from Fig. 5 that the orange curve (determined by (−ũ) = 10m 2 D ) intersects with p/T = 20 at θ = 0.27, meaning the condition (3.5) will be satisfied for θ ≥ 0.27. We have therefore plotted P (θ) and N hard (θ min ) for θ ≥ 0.27 and θ min ≥ 0.27 respectively.
Our results can also only be trusted where N hard (θ min ) 1, since if N hard (θ min ) approaches 1 this tells us that we cannot neglect multiple scattering. Including only single scattering, as we have done, is only valid where N hard (θ min ) 1. We see in the right column of Fig. 6 that, for the values of parameters used, N hard (θ min ) < 0.1 wherever we have shown our results, e.g. wherever we have plotted the red solid or dashed curves. This means that, for κ = 30, everywhere that the condition (3.5) is satisfied we also have N hard (θ min ) < 0.1. If we had chosen a larger value of κ this would not have been the case, and we would have needed to enforce a separate constraint.
At values of θ and θ min that are smaller than those for which we have plotted our results for P (θ) and N hard (θ min ), multiple scattering will become important, making the calculation much more difficult. At small enough angles, where many scatterings contribute, the result will simplify as the probability distribution for the transverse momentum transfer P(q ⊥ ) becomes a Gaussian at small enough q ⊥ [10]. As we noted in the Introduction, this is also the result that must be obtained in the regime in which the momentum transfer is small enough that the hard parton sees the QGP only as a liquid, without resolving the partons within it. The transverse momentum picked up by an energetic parton traversing a strongly coupled liquid is Gaussian distributed. Hence, whether we think of this from the perspective of a hard parton traversing a strongly coupled liquid or from the perspective of multiple scattering in a weakly coupled plasma, at small q ⊥ we expect P(q ⊥ ) to take the form where we have written the width of the Gaussian asqL, denoting the mean transverse momentum squared picked up per distance travelled byq as is conventional. The physics of multiple soft scattering in a weakly coupled plasma or the physics of how an energetic probe "sees" a liquid then determine the value of the parameterq. Following Ref. [28], it is convenient to introduce a dimensionless parameter K to parametrize the magnitude ofq viaq We can then use Eq. (C.3)) from Appendix C to convert P GA (q ⊥ ) in Eq. 3.10 to a probability distribution P GA (θ) for the angle θ, obtaining where we have used the approximation q ⊥ ≈ p in sin θ, valid for small θ where p ≈ p in . Hence, the behavior that we expect for P (θ) is that it should take the form (3.12) at small θ, for some value of K, and should then have a tail at larger angles θ that is due to single scattering of partons in the QGP, a tail that we have calculated and that is illustrated by the red curves in Fig. 6. To get a sense of how this might look, in Fig. 6 in addition to plotting the results of our calculations, in red, we have plotted P GA (θ) from (3.12) for two benchmark values of K, namely K = 5 and K = 12. (K = 5 is the value obtained by the JET collaboration [40] upon comparing calculations of observables sensitive to parton energy loss in a weakly coupled framework in which K controls energy loss as well as transverse momentum broadening. K = 12 is half of the value found for an energetic parton traversing the strongly coupled plasma of N = 4 SYM theory [8][9][10]; since this theory has more degrees of freedom than QCD, its strongly coupled plasma would have a larger value of K than the strongly coupled QGP.) Plotting P GA (θ) in addition to our own results in Fig. 6 is useful for two reasons. First, it helps us to imagine how these quantities may behave in a more complete calculation, following one of the black curves at small angles and then behaving along the lines of our results in red at large angles where single Molière scattering off partons in the QGP dominates. Second, by comparing the red curves to the black curves we can get a sense of at how large values of θ single hard scattering off partons in the QGP is likely to dominate over multiple soft scattering or the physics of the strongly coupled liquid. From the middle panels of Fig. 6 we see that the situation is rather clean for incident partons with p in = 100 T = 40 GeV: as long as we look at partons that scatter into a direction that deviates from the direction of the incident parton by θ > 0.3, we will be seeing Molière scattering. And, the probability for scattering at these angles can be quite substantial. If it proves possible to look at the scattering of even higher energy jet partons, for example as in the bottom panels of Fig. 6 where we take p in = 250 T = 100 GeV, Molière scattering and multiple soft scattering or the physics of the strongly coupled liquid separate even further. And, the probabilities for seeing large angle scattering remain quite significant as long as one looks for scattered partons with p > p min for a small enough p min , for example p min = 25 T = 10 GeV as in the solid red curves in the bottom panels of Fig. 6. The situation is less clear when we look at incident partons with p in = 25T = 10 GeV, in the top panels of Fig. 6. We see there that in order to see a red curve above the black curves at a probability above 10 −3 we need to look at the solid red curves, meaning we need to look at scattered partons with energies down to p min = 10T = 4 GeV and we need to look at rather large angles. It will be hard to separate final state hadrons coming from scatterings with these parameters from final state hadrons coming from the wake that the jet leaves behind in the plasma.
To the extent that one can draw conclusions from a calculation of scattering off a brick of plasma with T = 0.4 GeV and ∆t = 3 fm, our results suggest that experimentalists should look for observables sensitive to phenomena along the following lines: 40 GeV partons within a jet scatter off a parton in the plasma, yielding partons with energies greater than 8 GeV at angles θ > 0.5 with probability 1/100 and at angles θ > 0.8 with probability 1/1000. We would be happy to work with anyone planning future experiments to provide them with results along these lines for other values of the various parameters. But, the real path to predictions for observables is to take our results, formulated as in Section 3.2, and to incorporate them into a jet Monte Carlo analysis that also includes a realistic description of the expanding cooling droplet of plasma produced in a heavy ion collision.

Summary and outlook
We have analyzed the thought experiment depicted in Fig. 1 in which an incident parton (quark, antiquark or gluon) with energy p in traverses a brick of QGP with some thickness L and some constant temperature T and computed the probability distribution F (p, θ) for finding a parton (quark, antiquark or gluon) subsequently with an energy p that has been scattered by an angle θ relative to the direction of the incident parton. By integrating over p we obtain P (θ), the probability for finding a parton with p > p min scattered by θ, and then by integrating over θ we obtain N hard (θ min ), the number of hard partons scattered by an angle θ > θ min . We only consider binary collision processes in which the incident parton strikes a single parton from the medium, once. Because we neglect multiple scattering, our results are relevant only in the kinematic regime in which N hard turns out to be small, which means at large momentum transfer, and in particular at large values of θ. Because we are focusing on binary collisions with a large momentum transfer, for our medium we choose a gas of massless quarks, antiquarks and gluons with Fermi-Dirac or Bose-Einstein momentum distributions. Although we have ensured that we work only in a regime in which the momentum transfer in the binary collisions that we analyze is large enough that it is reasonable to neglect the Debye masses of the partons in the plasma, choosing their momentum distributions as if they were a noninteracting gas is relevant only as a simple benchmark. Ultimately, we look forward to the day when experimental measurements that are sensitive to the Molière scattering that we have analyzed can be used, first of all, to provide tangible evidence that the liquid QGP that we see today really is made of point-like quarks and gluons when probed at high momentum transfer and, second of all, via deviations from predictions based upon our calculations, to learn about the actual momentum distributions of these quarks and gluons. This would realize the vision of using the scattering of jet partons to learn about the microscopic structure of liquid QGP and would be analogous to learning about the parton distribution functions for QGP.
Realizing this vision will require incorporating the results of our calculations within jet Monte Carlo analyses in which realistic jets are embedded within realistic hydrodynamic models for the expanding cooling droplets of QGP produced in heavy ion collisions. Our results as we have obtained them here are based upon a thought experiment and cannot be compared directly to experimental data. It would be interesting to use comparisons between our results and results from Monte Carlo analyses in which binary collisons are already included (set up with jets probing a static brick like ours) to identify observable consequences of large-angle scattering. With a view toward Monte Carlo calculations which do not currently include binary collisions, we have presented our results in Section 3.2 in a form in which they could be incorporated into such analyses.
We note that we have worked only to leading order in perturbative QCD. This can certainly be improved upon in future work. However, it is our sense that incorporating these results in more realistic (Monte Carlo) modeling of jets probing more realistic (hydrodynamic) droplets of QGP is a more immediate priority than pushing our "brick calculation" beyond leading order.
Although the road ahead toward quantitative comparison to experimental measurements is a long one, our present results can already be used to reach several interesting qualitative conclusions. Perhaps the most interesting aspect of our results from a theoretical perspective is the importance of channels that are not Rutherford-like. It is only at small angles θ (where high momentum transfer requires large p in , as in previous calculations done in the p in → ∞ limit) where the dominant binary collision process is the Rutherford-like process where the parton that is detected is the incident parton, scattered by an angle θ. We have checked that our results reproduce the results of previous calculations in this regime. At the larger values of θ that are of interest, though, processes in which the detected parton is either a parton from the medium that received a kick or a parton that was produced in the collision (cf gg ↔ qq) are much more important. Consequently, also, we realize that at the values of θ that are of interest it is important to look for scattered partons that are still hard but that have substantially smaller energy than the incident parton.
Even though quantitative predictions for experimental measurements await further steps down the road ahead as we have discussed, the second place where our results are of qualitative interest is in the context of gauging what sorts of observables experimentalists should aim to measure. To get a sense of this, in Section 3.4 we have considered a brick of plasma that is 3 fm thick and that has a temperature T = 0.4 GeV, and have set g s = 1.5, corresponding to α QCD ≈ 0.18. (This exercise can easily be redone with other values of these parameters.) With these values, we find that it would be quite a challenge to look for the Molière scattering of jet partons that have p in = 10 GeV before they scatter. Doing so would require looking for observables that are sensitive to scattered partons with energies down to 4 GeV, and even if that were possible it would be hard to differentiate between partons scattering off particulate structures within the liquid QGP and partons picking up a Gaussian distribution of transverse momentum just from soft interactions with the liquid QGP. The picture is much more promising if instead we look for the Molière scattering of jet partons that have p in = 40 GeV (or more). before they scatter. Molière scattering is the dominant contribution if we look for scattering with θ > 0.3. And, although these processes are rare (they have to be rare in the regime in which they are the dominant contribution), the relevant probabilities are not tiny, given the high statistics data sets for jets in heavy ion collisions anticipated in the 2020s. For an incident parton with p in = 40 GeV, the probability of seeing a scattered parton with p > 8 GeV deflected by θ > 0.5 (θ > 0.8) is around 1/100 (1/1000). Getting a sense of the kinds of values of p in , p and θ where one should look, and a sense of the scale of the probability for the Molière scattering that one is looking for, should be of value both to experimentalists planning future measurements and to theorists exploring which jet substructure observables may be the most promising to measure.

A Full Boltzmann Equation
In this Appendix, we present a full derivation of the Boltzmann equation describing the evolution of the phase space density. After presenting the general formalism, we show how we recover Eq. (2.4) in the limit of a single binary collision. The expression (2.4) is then the starting point for the derivation of all of our results.
Beginning with greater generality than in Eq. (2.4), we define the phase-space distribution as follows F a (p, λ a , χ a ) ≡ Phase-space probability of finding a parton of species a (u, d, s,ū,d,s or g) with momentum p, helicity λ a and color state χ a . (A.1) This function depends on the time t, but we leave this dependence out of our notation for the present. The Boltzmann equation describing the time-evolution of this phase-space distribution takes the schematic form On the left-hand side, we have the time derivative of the phase-space distribution. On the right-hand side, we have the reason why such a function evolves with time: (binary) collisions. The collision operator C a is a functional that depends on the phase-space distribution of the parton a under consideration.
The collision operator has two distinct contributions that we denote via because there are two different ways to alter the distribution: • a binary collision produces the parton a with momentum p in the final state, which is accounted for by C (+) a [F a (p, λ a , χ a )] that appears with a plus sign; • a parton a with momentum p in the initial state is involved in a binary collision, which is accounted for by C (−) a [F a (p, λ a , χ a )] that appears with a minus sign.
We are interested only in the phase space distribution for the momentum, meaning that later in our derivation we will average over the helicity and color states.

A.1 Collision Operator for a Specific Binary Process
The expression in Eq. (A.2) is very general. Once we have a specific theory for the interactions mediating the binary collisions (in our calculation, QCD), we can derive an explicit expression for the collision operator. In this Appendix we shall not specialize that far, considering here a specific binary process In our derivation, we account for this process going both from left to right and from right to left. In the former case, it contributes to C (−) a (it can destroy a parton a with the given momentum p), whereas in the latter case it can contribute C (+) a . The explicit expressions for both contributions are given by: Here, the sign of the ± in a term like [1 ± F c ] is positive for bosons and negative for fermions, and these factors describe the Bose enhancement or Pauli blocking for the particles produced in the final state. Note that we are using the short-handed notation The squared matrix elements |M ab→cd | 2 are for a given polarization and color configuration, and we explicitly sum over such configurations for the states b, c, d. The prefactor with the δ cd accounts for the case where c and d are identical particles, where we must not double count. Upon assuming CP invariance, valid in particular for strong interactions, we have the identity Thus we can combine the two contributions together, and write the collision operator as The total collision operator appearing in the Boltzmann equation for species a is then the sum of all the individual ones accounting for each binary collision process in which a is involved: where n is the index labeling the different processes (e.g. n = ab ↔ cd).

A.2 Average over helicity and color states
We are not interested in keeping track of helicities and colors, since they cannot be resolved by the detector. We will average over them by introducing a new distribution The degeneracy factor ν a is the sum of all helicity and color configurations. Upon applying this definition to the Boltzmann equation in Eq. (A.2) we find ∂ f a (p) ∂t = 1 ν a λ,χ ∂F a (p, λ a , χ a ) ∂t = 1 ν a n λaχa C a [F a (p, λ a , χ a )]| n . (A.12) Focusing on a specific binary process n = ab ↔ cd, we can then write the explicit expression Finally, we replace all the distributions occurring on the right-hand side with the those averaged over polarizations and colors as defined in Eq. (A.11). In doing so, we are assuming that the medium has no net polarization and no net color charge. We also average over the helicity and color state of the incoming parton probing the medium. We end up with the expression where we have defined the collision operator accounting for the process ab ↔ cd by Here, we have introduced the matrix elements in the form that we use them in Section 2, namely |M ab↔cd | 2 ≡ λ abcd χ abcd |M ab↔cd | 2 , (A. 16) summed over initial and final polarizations. For the QCD processes of interest to us, these matrix elements are given in Table 1 of Section 2.3. The full evolution of the averaged phase space distribution reads with the sum accounting for all possible processes affecting the phase space distribution of the parton a.

A.3 Single Scattering Approximation
The results found so far allow for the possibility of multiple binary collisions. Next, we make the further assumption that the incoming probe scatters off a constituent of the medium just once before escaping on the opposite side. In order to so do, we find it convenient to employ the decomposition f a (p) ≡ n a (p) + f a (p) , (A. 18) where the "soft" thermal part n a is constant in time, and the residual piece can be interpreted as the "hard" part of the distribution, describing energetic partons. The collision operator for a specific binary process ab ↔ cd, whose explicit expression is given in Eq. (A.15), can then be simplified as follows. First, we observe that once we employ the decomposition in Eq. (A.18) the contribution with only thermal distributions vanishes because of the detailed balance principle. Next, we observe that we are only interested in collisions in which an energetic parton collides with a soft parton from the medium. (If we included many collisions, somewhere downstream from the first collision an energetic parton might collide with another energetic parton. This is impossible in the first collision, which for us is the only collision.) We furthermore observe that in the "hard region" of phase space (i.e. p T ) where we shall focus, we have n a (p) 1 and f a 1 also. Looking at the second and third lines in Eq. (A.15), describing the process cd ↔ ab, we find that via these considerations they simplify: (A.20) Upon making this single scattering assumption, and upon noting that the medium thermal distribution functions for our brick of noninteracting QGP are known and independent of the time, the Boltzmann equation takes the form The sum still runs over all the different binary processes involving species a, and the collision operator takes the final form where we have now added explicit mention of the time dependence to our notation. Focusing on just a single binary process ab ↔ cd, the solution reads (A.24) The probability for the parton a to have momentum p at the time t I + ∆t, namely the left-hand side of the above equation, is the sum of two contributions, the two terms on the right-hand side. First, we could already have a parton a with momentum p at the initial time t I and then have no further momentum transfer. Or we could achieve a momentum p at the time t I +∆t by a binary scattering. In this paper, we only care about the latter, since we are studying binary collisions with large momentum transfer resulting in the presence of a parton with a large angle deflection with respect to the incoming direction. That is, we shall always choose p to point in a direction that differs from that of the incident parton by some large angle θ, meaning that there is no parton a with momentum p at t I . Thus, for our purposes we need only consider the contribution in the last line of Eq. (A.24), which then becomes our Eq. (2.4) in the main text after summing appropriately over different processes. This is the key result of this Appendix, and the starting point for our analysis in Section 2.

B Phase space integration
where we used the spatial delta function to perform the integration over d 3 k and then shifted variables p to q 1 ≡ p − p, where φ 1 is the angle between the (q 1 , p) plane and the (q 1 , k ) plane, and where cos θ k q 1 and cos θ pq 1 denote the angles between k and q 1 and between p and q 1 , respectively. The integration over the azimuthal angle of q 1 has been performed trivially. To further integrate over the remaining delta function in (B.1), we follow the integration technology of Ref. [69] (see also Refs. [63][64][65]) and consider the identity The two delta functions in (B.2) can be recast as where we have used kinematic relations p = p 2 + q 2 1 + 2p q 1 cos θ pq 1 , k = (k ) 2 + q 2 1 + 2k q 1 cos θ k q 1 , The integration over cos θ pq 1 and cos θ k q 1 in Eq. (B.9) can be performed trivially when the following kinematic constraints are satisfied: The constraints (B.7) and (B.8) imply that |ω 1 | ≤ q 1 ≤ 2p + ω 1 and k ≥ q 1 −ω 1 2 . We consequently have where ∆θ 1 = θ 1 − θ and θ 1 denotes the angle between the directions of p and p in . Here, we used the relation which follows from the fact that We now substitute Eq. (B.9) into Eq. (2.15a) to obtain: To proceed, we express f I (p ) in Eq. (2.2) as a function of ∆θ 1 and ω 1 : Therefore, the integration over ω 1 and ∆θ 1 in Eq. (B.12) can be performed directly after substituting Eq. (B.13) into Eq. (B.12). As a result, we replace ω 1 with ω, q 1 with q, ∆θ 1 with θ and identify t, u witht,ũ as defined in Eq. (2.30). After relabeling the dummy integration variables k with k T and φ 1 with φ, we eventually arrive at Eq. (2.28).
The derivation of Eq. (2.33) follows similar steps.

B.2 Integration over
We then have:

C Comparison with previous results
C. 1 The relation between P(q ⊥ ) and P (θ) To elucidate the connection with previous studies [10,61,66] in which the two-dimensional probability distribution for the transverse momentum of the outgoing parton, P(q ⊥ ), has been computed, we need to relate this quantity, normalized as to the probability distribution P (θ) for the angle θ that we compute. Since the previous studies all work in a limit in which p in is large and θ is small, energy loss is negligible in these studies, i.e. p ≈ p in , and hence q ⊥ = p in sin θ .

(C.3)
We shall use this expression, with p replaced by p in , in Eqs. (3.4) and (3.12).
It simplifies the explicit comparisons that we shall make in Section C.2 if there we work in the small-θ limit in which q ⊥ ≈ p in θ and J ⊥ reduces to In Section C.2 our goal will be to check whether the following relation holds: lim θ→0 J ⊥ P (θ) = P single (q ⊥ = p in θ) , (C.5) where P (θ) is the result of our calculation and P single (q ⊥ ) is one of the results from Refs. [10,61,66] for P(q ⊥ ) due to a single binary collision.

C.2 Previous results, compared to ours
The expression for P(q ⊥ ) due to a single binary collision, P single (q ⊥ ), has been obtained in the limit m D q ⊥ T , by Aurenche, Gelis and Zaraket (AGZ) [66], who showed that (in our notation) in this regime, and in the limit q ⊥ T by Arnold and Dogan (AD) [61], who showed that in this regime. Each of these expressions is a limiting case of the more general expression for P single (q ⊥ ) computed by D'Eramo, Lekaveckas, Liu and Rajagopal (DLLR) [10]. In the limit q ⊥ m D their result can be written as (see Eq. (5.2) and Eq. (5.15) of Ref. [10]): where Im Π T,L are the imaginary parts of the gluon longitudinal and transverse self energy in QGP. To obtain Eq. (C.8), we have used the relation q z ≈ ω which is valid in the limit where we have used the relation (C F d F /d A ) = 1/2 and C A = N c . One important consequence of Eq. (C.18), in particular the fact that c = 0, can be found by substituting these generic results together with Eq. (C.13) which is valid in the small-θ limit into Eqs. (2.19), (2.21), (2.23), (2.25), (2.27) and discovering that F G→G (p, θ) F G→Q (p, θ) , F G→Q (p, θ) and F Q→Q (p, θ) F Q→G (p, θ) , F Q→Q (p, θ). That is, This simply reflects the fact that in the small-θ limit, Rutherford-like scattering (in which the parton that is detected is the incident parton after scattering) is much more important than other channels. We will focus on F G→G (p, θ) and F Q→Q (p, θ) from now on and write explicit expressions for them by substituting Eqs. (C.13) and (C.18) into Eq. (2.25), obtaining where we have used ν q = 2d F , ν g = 2d A , C A = N c . Comparing Eq. (C.20b) with Eq. (C.20a), we obtain the relation We can now compute the left-hand-side of (C.5) for an incident gluon by substituting Eq. (C.20a) into Eq. (C.9). We find that where we have used Eq. (C.14). For an incident quark, the resulting P (θ) can be obtained by replacing C A with C F thanks to the relation (C.21). Eq. (C.22) is a central result of this Appendix, as it will allow us to compare our results to those obtained previously in the limits in which such comparisons can be made. In order to compare our result to the AGZ result (C.6) [66] we must evaluate our expression (C.22) in the limit ω, q ⊥ T . We see from Eq. (C.11) that in this limit q T . Since the characteristic k T is of the order of T , we can set ω = 0 in n B.E (k T + ω) and where the Debye mass m D is given by Eq. (3.9). We observe that, as advertised earlier, Eq. (C.26) is equivalent to the AGZ result (C.6) through the relation (C.5). It is worth noting that the dominant contribution to the integration in Eq. (C.25) comes from ω ∼ q ⊥ p in , which justifies taking the limit ω/p in 1 in F C→all (p, θ). We now turn to comparing our result (C.22) to the DLLR result (C.8) [10]. To simplify the discussion, we will only include the contribution coming from thermal scatterers which are gluons. This amounts to setting N f = 0 in Eq. Correspondingly, we will only include the contribution to the gluon self-energy Π L,T in Eq. (C.8) that comes from gluon loops, and show that the resulting P DLLR single (q ⊥ ) is equivalent to Eq. (C.27) through the relation (C.5). The comparison upon including the contribution coming from fermionic thermal scatterers (quark and antiquark) is quite similar.
To proceed, we write the explicit expressions for Im Π L and Im Π T coming from the gluon loop as given in Ref. [10]: Finally, we substitute Eq. (C.30) into the DLLR result (C.8). It is now transparent that our expression (C.27) is equivalent to Eq. (C.8) through the the relation (C.5).
Noting that it has been demonstrated in Ref. [10] that the AD result (C.7) is obtained from the DLLR result (C.8) in the q ⊥ T limit, this concludes our verification that our result, in particular in the form (C.22), reduces to the previously known AGZ, AD and DLLR results in the appropriate limits.