Transverse Momentum Broadening of a Jet in Quark-Gluon Plasma: An Open Quantum System EFT

We utilize the technology of open quantum systems in conjunction with the recently developed effective field theory for forward scattering to address the question of massless jet propagation through a weakly-coupled quark-gluon plasma in thermal equilibrium. We discuss various possible hierarchies of scales that may appear in this problem, by comparing thermal scales of the plasma with relevant scales in the effective field theory. Starting from the Lindblad equation, we derive and solve a master equation for the transverse momentum distribution of a massless quark jet, at leading orders both in the strong coupling and in the power counting of the effective field theory. Markovian approximation is justified in the weak coupling limit. Using the solution to the master equation, we study the transverse momentum broadening of a jet as a function of the plasma temperature and the time of propagation. We discuss the physical origin of infrared sensitivity that arises in the solution and a way to handle it in the effective field theory formulation. We suspect that the final measurement constraint can only cut-off leading infrared singularities and the solution to the Markovian master equation resums a logarithmic series. This work is a stepping stone towards understanding jet quenching and jet substructure observables on both light and heavy quark jets as probes of the quark-gluon plasma.


Introduction
Jets are sprays of collimated particles produced in high energy collisions of hadrons and/or electrons. Their formation starts with a highly virtual parton generated from an initial hard scattering, followed by subsequent parton cascade and fragmentation. Jet production can be studied via perturbative QCD due to the large scales involved, for example, the virtuality of the initial parton produced. Therefore the calculation of the initial production of jets can be well-controlled theoretically, which makes jets powerful tools to probe the properties of the quark-gluon plasma (QGP) in heavy ion collisions.
Jet production is modified in heavy ion collisions, compared with that in proton-proton collisions, due to the jet-medium interaction. Jet quenching, a phenomenon of suppression of particles with high transverse momenta, has been studied intensively theoretically long before [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] and recently observed in experiments at both Relativistic Heavy Ion Collider (RHIC) [21][22][23][24] and Large Hadron Collider (LHC) [25][26][27]. The suppression mechanism is mainly the energy loss when jets traverse the hot medium. Both collisional and mediuminduced radiative energy loss contribute, but the latter dominates at high energy. The key to understand jet quenching and jet substructure modifications in heavy ion collisions is to understand how the jet interacts with the expanding medium. There has been tremendous theoretical effort to study the jet energy loss mechanism (see Refs. [28][29][30][31] for recent reviews). But this is not a simple problem because it involves multiple scales such as the jet energy, the transverse momentum with respect to the jet axis and thermal scales of the QGP. Furthermore, in current heavy ion collision experiments, the temperature achieved fits roughly the range 150 − 500 MeV, and may not always be a perturbative scale. Thus, a fully weak coupling calculation may not be valid. A hybrid model has been developed to address this problem [32][33][34][35][36], in which the initial jet production and vacuum-like parton shower are calculated perturbatively, while the subsequent jet energy loss in the medium is calculated by mapping the field theory computation in the strong coupling limit to a weak coupling computation in the classical gravity theory [37][38][39][40][41][42], i.e., by using the AdS/CFT correspondence [43].
From the perspective of field theory, a powerful tool to deal with multi-scale problems is effective field theory (EFT). The EFT that is particularly useful for jet studies is Soft-Collinear Effective Theory (SCET). There are also formulations of SCET (known as SCET G ) treating the Glauber gluon, which is a type of mode appearing in forward scattering, as a background field induced by the medium interacting with an energetic jet. By making use of the collinear sector of the corresponding EFT, this formalism has been used to address the question of jet quenching in the medium [44][45][46][47][48]. In the same spirit, a new EFT for forward scattering has been developed recently [49] which also uses the Glauber mode to write down contact operators between the soft and collinear momentum degrees of freedom. The partons from the thermal QGP are generally soft, when compared with an energetic jet, which can be described by a collinear mode.
Another theoretical challenge in understanding the jet-medium interaction is the quantum interference effect. For example, in the process of the medium-induced single radiation from a high energy parton, multiple transverse momentum kicks from the medium can suppress the radiation spectrum by destructive interference, a phenomenon known as the Landau-Pomeranchuk-Migdal (LPM) effect. Progress in understanding the LPM effect in the single radiation has been achieved in recent years [50][51][52]. Extension to studying the LPM effect in multiple splittings has been explored in simplified cases [53,54]. Furthermore, when two collimated partons are close to each other spatially, the medium may not be able to resolve them completely. So they may lose energy coherently as a single parton. This interference effect caused by the finite resolution power of the QGP is also important and can change jet substructure observables dramatically [36].
To take into account the interference effect systematically, one can keep track of the time evolution of the system's density matrix. This can most easily be done by using the open quantum systems formalism (for introductory books, see [55,56]). For jets inside a QGP, if we only focus on jet observables, the jet can be treated as an open quantum system interacting with a QGP bath. The application of the open quantum system formalism in heavy ion collisions has been thriving in the study of color screening and regeneration of quarkonium [57][58][59][60][61][62][63][64][65][66][67][68]. There, the heavy quark-antiquark pair in the color singlet interacts with the medium destructively when they are close. Great progress in the understanding of quarkonium in-medium dynamics has been achieved by combining potential nonrelativistic QCD (pNRQCD [69][70][71], an EFT of QCD) and the open quantum system formalism [72][73][74][75]. For example, a semiclassical Boltzmann transport equation of quarkonium in the medium has been derived, under assumptions that are closely related with a hierarchy of scales [75][76][77].
We would like to combine the forward scattering EFT recently developed within the formalism of SCET with the open quantum system formalism and explore its physical implications on the jet-medium interaction. As a first step, we will study in this paper, the transverse momentum broadening of a high energy parton moving through the medium. For simplicity, we will assume the plasma temperature is high enough so the weak coupling calculation is valid. We will leave the inclusions of nonperturbative effects and radiation into the calculation to future work.
The use of an EFT formalism allows us to write a simple and hence easily calculable description of the system. On the other hand, as alluded to earlier, we hope that the use of an open quantum system approach will allow us to easily keep track of the quantum interference effects and construct a new class of calculable observables. The long term goal here is to develop a theoretically robust formalism for computing jet substructure observables for both light parton and heavy quark jets. For example, the bottom quark jets have been identified as an effective probe of the QGP medium and will be experimentally studied at LHC, as well as by the sPHENIX collaboration at RHIC. There has been recent work on computing jet substructure observable for heavy quark jets in the context of protonproton collisions [78,79]. The objective would then be to compute the same observables in heavy ion collisions and study modifications caused by the medium. This paper is organized as follows. In Section 2 we discuss the importance of the forward scattering regime in jet-medium interactions. Based on this, we then introduce the physical system that we wish to study and the relevant physical scales that play an important role in its description. The next Section 3 reviews the basics of SCET and Glauber EFTs as tools to model jet propagation in the QGP. We then introduce the concept of open quantum systems and derive a master equation for the jet density matrix from the Lindblad equation in Section 4. The master equation is solved analytically and the transverse momentum distribution of a jet is studied in Section 5 along with a comparison with previous results in literature. We also discuss possible infrared (IR) divergences that show up in the Markovian limit. Finally we conclude and discuss future directions in Section 6.

Relevant Hierarchy of Scales
In this section we examine the dominant interaction of a jet with a QGP medium and discuss the possible hierachy of scales that can appear as a function of the jet energy and QGP temperature.

Dominance of Forward Scattering Regime
We can gain intuition about the dominant interaction of a parton traversing a QGP medium by examining a 2 → 2 scattering. Consider the simple example of 2 → 2 scattering e − µ − → e − µ − . The lowest order Feynman diagram is just a t channel photon exchange. The differential cross section for this process has the form [80] where E cm is the center of mass energy and θ is the angle between the final and initial state electron/muon. α EM is the electromagnetic coupling. This cross section has a singular behavior as θ → 0 This singularity over the phase space is not integrable and is in fact a physical singularity that arises due to the infinite range of the Coulomb potential. In real experiments, however, this singularity gets cut off by some IR scale such as a dynamically induced photon mass or a finite interaction region introduced by localized beam wave packets at a finite impact parameter. It is worth noting that most of the contribution to the cross section comes from this small angle region of phase space. In the case of scattering inside the QGP, the corresponding 2 → 2 scattering would be forward scattering of quarks/gluons mediated by a gluon. For a QGP in thermal equilibrium with a temperature T , the interactions in the medium induce an effective gluon mass (Debye mass) m D which acts as an IR cut-off scale. As we will do in this paper, we may also impose some cuts on the final state measurements which will act to regulate some of the IR divergences. We wish to develop an EFT description in this region of forward scattering by expanding in the small scattering angle θ, which will also be the power counting parameter of our EFT. How exactly this parameter is related to the physical scales of the system will depend on the measurement that we impose on the final state. Roughly speaking, for a high energy parton (quark or gluon) with an initial energy Q, the angular parameter θ ∼ Q ⊥ /Q measures the transverse momentum of the final state with respect to its initial direction (considered as the longitudinal direction). We will discuss this in detail in the following subsection.

Jets in Quark-Gluon Plasma
The QGP at vanishing chemical potential 1 is characterized by its temperature T in thermal equilibrium. In this paper we will mostly be concerned with light quarks so all the partons are considered massless at the level of the Lagrangian. The finite temperature interactions between the partons induce a dynamical gluon mass m D , which is of the order of gT in perturbation. Here g is the strong coupling at the scale T . The Debye mass provides a screening effect and effectively shortens the strong interaction range to be of the order of 1/m D . The perturbative description of the QGP works well when the temperature T is far above the confinement scale Λ QCD , which is a dynamically generated scale of QCD.
The system we want to study is a highly energetic jet traversing a region of the QGP. The energy of the jet, Q, will be the hard scale in our process and is assumed to be much larger than all the other scales in the problem. In this paper, we will focus on the leading interaction of the jet with the QGP medium and leave the vacuum as well as medium-induced splitting to future studies. So effectively our jet is described by a leading parton. The aim is to calculate the final transverse momentum broadening of the jet when it comes out of the QGP. For our purpose in this paper, we are going to impose a measurement constraint on the final transverse momentum Q ⊥ of the jet, by concentrating on the forward scattering region, i.e., λ ∼ θ 1. Here λ is the power counting parameter for our EFT description. If the interaction of the jet with the medium is a single coherent scattering, then we can write λ ∼ θ ∼ Q ⊥ /Q. In this paper, however, we allow the jet to have a series of mutually incoherent interactions with the medium. In any of the intermediate incoherent interactions, a smaller transverse momentum p ⊥ Q ⊥ can be exchanged between the jet and the medium so along as the net value adds up to Q ⊥ . This means that our power counting parameter satisfies λ Q ⊥ /Q, with the smallest transverse momentum being effectively cut-off by m D . We will discuss this in detail in Section 5. Since the transverse momentum exchanged is always small compared to the jet energy, the jet (leading parton) can be treated as a collinear particle throughout its evolution. Medium partons will be treated as soft modes that carry energy and momentum of the order of the QGP temperature T . We are always going to work in a regime where T Q.
We can use the final measurement on Q ⊥ to probe the physics at different scales (while maintaining Q ⊥ /Q 1, otherwise our EFT framework described below does not apply). The natural question then is how the IR scale Q ⊥ compares to all the other scales such as T , m D and Λ QCD . Several possible hierarchies are possible here and we now discuss each of them.

• High temperature
The simplest case is when the temperature T is high enough that both T and m D are perturbative scales. In other words, we have Q T m D Λ QCD . We can also have the possibility of a somewhat lower temperature such that Q T ∼ m D Λ QCD , where the Debye mass scale is still perturbative. Under these hierarchies, we 1 We leave the case of a non-zero baryon density to future studies.
can do a perturbative computation of the final observable Q ⊥ to probe all the scales greater than Λ QCD and much smaller than Q. We will primarily focus on this regime in this paper.
• Intermediate temperature In this case, we still have T high enough to be a perturbative scale. However, the scale m D ∼ Λ QCD is now nonperturbative, i.e., Q T m D ∼ Λ QCD . But we can still set up a perturbative EFT for Q ⊥ ∼ T . In fact the EFT that we construct in the high T regime in the first case will be valid for the Q ⊥ ∼ T case here as well. On the other hand, if we want to probe the physics at lower scales, Q ⊥ ∼ m D , then we have to take into account nonperturbative effects.
• Low temperature Finally we have the low temperature regime in which case all our scales T ∼ m D ∼ Λ QCD are nonperturbative. In this case, we can still develop a perturbative EFT description for the case Q Q ⊥ T , and then appropriately make a transition to the nonperturbative regime of Q ⊥ .
As a side remark, we want to emphasize that jet observables with the same jet radius R can be very different at the RHIC and LHC energies. Since at LHC the collision energy is much higher, more energetic jets can be produced, whose energies Q can be much larger than those at RHIC with the same jet radius. Then the transverse momenta ∼ QR (with respect to the jet axis) of jets at RHIC and LHC can be very different, even though these jets are defined with the same jet radius R, and probably probe the physics at different scales. For example, it is likely that the transverse momentum at the LHC energy is in the perturbative regime while that at the RHIC energy sits in the nonperturbative regime. EFT approaches can help us to better understand the difference quantitatively.
Since we will concentrate on the region of forward scattering, we can use EFT tools already available in the literature to describe our system. One such a formalism that has been extensively used in collider physics is SCET. We now review the basics of SCET and discuss how it can be used in our problem.

Review of Soft-Collinear Effective Theory
SCET is a theory of both soft and collinear particles. Collinear particles have a large momentum along a particular light-like direction, while soft particles have a small momentum, and no preferred direction. For each relevant light-like direction, we define two reference vectors n µ andn µ such that n 2 =n 2 = 0 and n ·n = 2. The typical choice of n µ = (1, 0, 0, 1) andn µ = (1, 0, 0, −1) will be used below. The freedom in the choice of n, as in the case of the label velocity in Heavy Quark Effective Theory, is represented in the EFT by a reparametrization invariance [81,82]. Any four-momentum p can be decomposed with respect to n µ as The SCET is defined by a systematic expansion in terms of a formal power counting parameter λ 1, which is determined by the measurements or kinematic restrictions imposed on the QCD radiation. The momenta for different modes in the SCET scale as Collinear : n·p, n·p, p ⊥ ∼n·p 1, λ 2 , λ , Ultrasoft : n·p, n·p, p ⊥ ∼n·p λ 2 , λ 2 , λ 2 .
A theory with only collinear and ultrasoft modes is typically referred to as SCET I, while that with only collinear and soft modes is referred to as SCET II [83] 2 . In this paper, we will only be concerned with SCET II along with the Glauber mode.
In order to expand the fields in the full theory QCD around a particular direction, the momenta are decomposed into the labelp µ and residual k µ components Then for a collinear particle,n ·p ∼ Q andp ⊥ ∼ λQ, where Q is a typical scale of the hard interaction, while k µ ∼ λ 2 Q describes small fluctuations around the label momentum. Field modes with momenta of definite scaling in the SCET are obtained by performing a multipole expansion of the fields in the full theory QCD. SCET involves independent gauge bosons and fermions for each collinear direction A n,p (x), ξ n,p (x), which are labeled by their collinear direction n and their large label momentump, as well as (ultra)soft gauge boson fields A (u)s (x). Independent gauge symmetries are enforced for each set of fields, which have support for the corresponding momentum carried by that field [84]. Overlap between different regions is removed by the zero-bin subtraction procedure [85]. This ensures no double counting of momentum regions. The leading power SCET II Lagrangian that we shall be concerned with in the following, takes the form Here L G is the leading power Glauber Lagrangian [49], which describes the leading power coupling between the soft and collinear degrees of freedom through potential operators. We will discuss the Glauber interaction in more detail in the next section.
Hard scattering operators involving collinear fields are constructed out of products of Wilson line dressed fields that are invariant under collinear gauge transformations [86,87]. For example, the gauge invariant gauge boson operator is given by where D n⊥ is the collinear gauge covariant derivative, and W n (x) is a collinear Wilson line where P µ is an operator that returns the label momentum of the fields on its right. The collinear Wilson line, W n (x), is localized with respect to the residual position x so that B µ n (x) can be treated as local gauge boson fields from the perspective of the ultrasoft degrees of freedom. For the leading power calculation presented here, ultrasoft and soft fields will not appear explicitly in our hard scattering operators, other than through the Wilson lines in the field redefinition which is performed on each collinear sector. For a general representation, r, the ultrasoft Wilson line is defined by 3 where P denotes path ordering. This so-called BPS field redefinition has the effect of decoupling ultrasoft and collinear degrees of freedom at leading power [88], and it accounts for the full physical path of ultrasoft Wilson lines [89,90]. In the following, we will also need soft Wilson lines,

Glauber Mode for Forward Scattering
The process of near forward scattering is referred to as the Glauber exchange which involves the exchange of an off-shell gluon. The transverse momentum of this gluon (with respect to the forward direction) is parametrically larger than its longitudinal components so that |k ⊥ | 2 n · kn · k, which is different from the regime of a Coulomb exchange which satisfies |k| 2 (k 0 ) 2 . As a result, the Glauber mode is not a propagating mode and acts instantaneously along the light-cone time. 3 Here we give the explicit result for an incoming Wilson line. Depending on whether particles are incoming or outgoing, different Wilson lines must be used. When done correctly, the BPS field redefinition accounts for the full physical path of the particles [89,90].
A systematic study of the Glauber mode was carried out in [49] within the formalism of SCET. Depending on the process that we are interested in, the asymptotic (propagating) states that we deal with can be classified as collinear and soft modes as described in the previous section. Scattering processes at colliders can usually be factorized in terms of these modes that separate the momentum fluctuations at different scales. The Glauber is an off-shell mode that mediates between either two collinear, two soft or one collinear and one soft modes, thus violating factorization. These factorization violating interactions can be captured via effective operators in the Glauber Lagrangian L (0) G . At leading power, these operators have been derived in [49].
Our interest lies in utilizing this EFT formalism to study the near forward scattering of a jet inside a QGP in thermal equilibrium. For the 2 → 2 scattering process, the region of near forward scattering dominates the total cross section and various approximations can be made to simplify the result at leading power in the expansion parameter. In the current case the expansion parameter is the small scattering angle θ.
We can now discuss the effective interaction between the jet leading parton and the medium in the forward scattering region. The thermal QGP is mainly composed of soft particles whose energies and momenta are on the order of T . Their momenta p s scale uniformly in our expansion parameter λ ∼ θ 1, so we can write where Q is the hard scale in our process, i.e., the jet energy. The jet is composed of highly energetic collinear particles that are moving along the light-like direction n. The scaling of their momenta, in the light-cone coordinate, can be written as As the collinear particles move through the medium, they interact with the soft medium particles mainly via forward scattering where both the collinear and soft particles maintain their momentum scaling after the scattering. This interaction is therefore mediated by the Glauber mode with the scaling Adding or subtracting a momentum of the Glauber scaling from the collinear or soft mode does not alter their momentum scaling.
The Glauber can be integrated out of the Lagrangian, which leads to effective operators coupling the collinear and soft degrees of freedom. The effective gauge invariant operators for quark-quark (qq), quark-gluon (qg or gq) and gluon-gluon (gg) interactions have been worked out in the Feynman gauge in Ref. [49] O qq where B is the color index and the subscripts n and s denote the collinear and soft operators respectively. The soft operators O s are constructed from the gauge invariant soft quark and gluon building blocks that are built out of the soft fields dressed with soft Wilson lines: where the soft Wilson lines ensure that the operators are invariant under soft gauge transformations.
The collinear operators are built out of the collinear building blocks. In this paper, we will only work with collinear quarks which are constructed from bare collinear quark fields dressed with collinear Wilson lines: where ψ is the standard four-component Dirac spinor. Since the soft momentum puts the collinear particle off-shell and off-shell modes have been integrated our in the construction of the EFT, the collinear fields do not transform under the soft gauge transformations. The effective Lagrangian density for the Glauber exchange then looks like in which the operators O i are listed in Eq. (3.13). One point to note here is that the collinear and soft operators O n and O s are separately gauge invariant. However, a simple calculation shows that a change in the gauge choice for the Glauber propagator in the construction would lead to a different form of the gauge invariant operators. The derivation of Eq. (3.13) in Ref. [49] chooses the Feynman gauge. In general, the operators will be of the form where ∆ µν would be the Glauber gluon propagator in the chosen gauge. The operators O nµ , O sµ would still be separately gauge invariant 4 but the form of these operators would change according to the Glauber gauge choice. The final result for any scattering amplitude would remain gauge independent. Thus it suffices to work with any specific gauge. From now on we will choose to work in the Feynman gauge since all the effective operators have been constructed.

Lindblad Equation for Open Quantum System
In this section we review the basic concepts of the open quantum system formalism, which can be used to describe the quantum dynamical evolution of an open subsystem in contact with an environment. Later in this section, we will apply this formalism to the case of a jet (subsystem) interacting with the thermal QGP (environment) via the effective operators described in the previous section. We begin with the microscopic derivation of the Lindblad equation which closely follows the discussion in Ref. [56]. We assume that the Hamiltonian of the total system (subsystem and environment) is given by where H S is the subsystem Hamiltonian, H E is the environment Hamiltonian, and H I contains the interactions between the subsystem and the environment. The interaction Hamiltonian is assumed to be factorized as follows: denote the subsystem and environment operators respectively. Here α denotes all relevant quantum numbers. (For local quantum field theory, the factorized form is generally true and α includes the spatial coordinates for which the summation means an integration.) In our case, α would be the particle type (q or g), the color and spatial coordinates of the effective operators. We can assume O Here ρ E is the density matrix of the environment. Each part of the Hamiltonian is assumed to be Hermitian.
The von Neumann equation for the time evolution of the total density matrix in the interaction picture is given by where The density matrix and Hamiltonians without any superscript are in the Schrödinger picture. In the above formula, we have used the fact [H S , H E ] = 0. We will omit the superscript "(int)" in the following discussion. The symbolic solution is given by where the evolution operator is and T is the time-ordering operator. We will assume the subsystem and the environment are weakly interacting. We further assume the initial total density matrix factorizes which is generally true for weakly-coupled systems (factorization breaking terms come at higher orders in the coupling). The environment density matrix is assumed to be in thermal equilibrium: where T = 1/β is the temperature of the thermal environment. If we expand the interaction to second order in perturbation and take the partial trace over the environment degrees of freedom, we obtain the Lindblad equation: Each term in the Lindblad equation is defined as where {|a } forms a complete set of states in the Hilbert space of the subsystem. The Lindblad equation has two non-trivial terms: First, the term proportional to σ ab corresponds to a unitary evolution induced by the interaction with the environment, in addition to the usual time evolution driven by the subsystem Hamiltonian. Second, the term proportional to γ ab,cd generates a non-unitary evolution through which the subsystem dissipates and loses coherence. In our case, the subsystem is an energetic jet propagating through the environment which is a QGP in thermal equilibrium with a temperature T . For current heavy ion collision experiments, the measured jet energy Q is much larger than the highest temperature achieved in the collision: Q T . The jet starts out as a single, highly virtual parton which then showers even in vacuum. Now inside the QGP, the parton shower is modified by its interaction with the medium. In the simple case of forward scattering, the energetic parton exchanges a small amount (small compared with the jet energy) of transverse momentum with the medium. Other physical processes are also possible: A virtual parton can radiate off soft on-shell partons with energy ∼ T which then become part of the medium (the wake of a jet). The virtuality of the parton may come from the initial production as in vacuum parton shower, or be developed by a sequence of soft kicks from the medium. In the latter case, the radiation is medium-induced. In this paper, we will take a first step towards understanding the full consequences of the subsystem-environment interaction. To that end, we will ignore the subsystem evolution in vacuum (which is already well understood). We will concentrate on the jet evolution induced by the forward scattering, i.e., the transverse momentum broadening of a jet. We will leave the inclusion of medium-induced radiation into our framework to future studies.
In order to describe the subsystem evolution in terms of the Lindblad equation, we need to know all the three Hamiltonians that are involved: • The subsystem Hamiltonian is the same as the SCET Hamiltonian for collinear particles in vacuum. Usually in vacuum the modes that appear in the SCET Lagrangian are determined by the measurement performed on the jet. Here both the measurements and the medium scales will determine the modes in SCET.
• The environment Hamiltonian describes the thermal QGP with the temperature T . The sources of the Glauber exchange from the medium are soft modes that scale as Q(λ, λ, λ).
• The interaction Hamiltonian describes the effective interaction between the collinear particles (subsystem) and the soft partons in the QGP (environment). For the near forward scattering region, the interaction happens via Glauber exchanges. The Glauber mode scales as Q(λ, λ 2 , λ) with the power counting parameter λ ∼ θ. The interaction between collinear and soft modes mediated via Glaubers is described in Section 3.
Finally, to apply the Lindblad equation in our study, we need the one-to-one correspondence between the SCET operators and the general operators used in the Lindblad equation. Here we list them: where q/g denotes the quark/gluon operator, A is the color index and x is the spatial coordinate of the fields. Now we will apply the Lindblad equation to study the interaction between a collinear quark and a soft quark via the Glauber exchange.

Unitary Evolution Part
We first compute the piece for the unitary evolution of the subsystem density matrix induced by its interaction with the thermal environment where non-unitary Lindblad terms are omitted for the moment. To evaluate the expression, we first write At the same time, we note that where we have swapped a and b, α and β, t 1 and t 2 in the last line. Then we can simply write where we have replaced the environment correlator β (t 2 ) T , since the environment density matrix is assumed to be in thermal equilibrium. The subscript T indicates the finite temperature. As we discussed earlier, the indices α and β include the spatial coordinates of the field operators. To show the spatial coordinates more explicitly, we rewrite this term as in which x i = (t i , x i ) and now the indexes A and B are just color indexes and no longer include the spatial coordinates. If we assume the thermal bath is homogeneous in space and time, we have where O (E) is given by the soft operators (4.15) that are dressed with soft Wilson lines. Here we introduced the finite temperature Wightman function in momentum space D AB > (k). We will evaluate this finite temperature correlator for soft quark operators in the imaginary time formalism. Details are provided in Appendix A.

Subsystem Transition
Now we discuss the computation of the other piece in Eq. (4.21): a|O Since the purpose of the current paper is to study the transverse momentum kicks, in which splittings are not taken into account, all the relevant states in the subsystem are one particle states. They can be specified by their momentum (we neglect the quark mass), color and spin. Spin does not flip in the Glauber exchange process, which is explained in Appendix B. We will average (sum over) the colors of the incoming (outgoing) states. So we will focus on the momentum and use |p to label the subsystem state |a . In the following, all the abstract state labels a, b, · · · will be replaced with the momenta p 1 , p 2 , · · · and the summation over a will be replaced with an integration over the momentum.
We can then evaluate the correlator of the subsystem (collinear quarks) operators as follows. For two collinear momenta p 1 and p 2 , we define the relevant transition between them as where the indices α and β include both color and the spatial coordinates. Inserting a complete set of one particle states leads to where q 0 = E q = |q|. This now represents a t channel process which we are interested in describing, since we focus on a collinear quark. Since we are concerned with a jet initiated by an energetic quark, we use the appropriate SCET interaction operators in Eq. (3.13) This can be rewritten as where now q 0 is not constrained.

Final Expression
We can now put all the pieces together for the driving term of the unitary evolution a,b σ ab L ab = − i 2 dp 1 dp 2 where the 1/N c factor comes from the average of the color of the incoming states. We have introduced the shorthand notation Integrating over x 1 and x 2 gives two delta functions for momentum conservation: δ 3 (p 1 − q − k) and δ 3 (p 2 − q − k). From this we can conclude p 1 = p 2 ≡ p. Since the external one particle states |p i (i = 1, 2) are on-shell, we further conclude E p 1 = E p 2 ≡ E p . Then the relevant time integrals in the limit t → ∞ become where ω = E p − k 0 − q 0 is the common energy conservation condition that the integral over t i yields. The t → ∞ limit is called the Markovian approximation, which is valid when the subsystem relaxation time is much bigger than the environment correlation time.
Physically it means during a typical evolution time of the subsystem, the environment loses all memory of the subsystem. Mathematically, it occurs when the environment correlator C αβ (t 1 , t 2 ) dies off quickly enough at large t i so the contribution to the time integral from the large time region is negligible. Markovian approximation is generally true when the subsystem is weakly coupled to a thermal bath. The argument is as follows: The typical energy scale of the thermal bath is the temperature. Its inverse gives the typical correlation time of the environment. The subsystem relaxation rate is roughly α s T , which is based on the leading order estimate. The inverse of the relaxation rate gives the relaxation time, which is much bigger than 1/T when α s is small. We will give an explicit estimate of the subsystem relaxation rate in Section 5.1.
Finally we obtain Since the various momenta in the formula scale in a specific manner, we can simplify the formula further by doing an expansion and keeping terms only at leading power in λ.
We notice that the unitary evolution driving term is diagonal in the subsystem state space. If the initial state density matrix is a pure state of a single parton ρ (4.31) So at least at leading order, no correction on the subsystem unitary evolution is generated from its interaction with the medium. In other words, at leading order, the single parton state energy is not corrected by the medium interaction.

Non-unitary (Dissipative) Evolution Part
The other terms in the Lindblad equation lead to a non-unitary or dissipative evolution for the subsystem density matrix: where higher order terms are neglected. Even though the evolution is non-unitary, it preserves the properties of a valid density matrix, i.e., hermiticity, positivity and unity trace (the total density of all states is conserved). In this case the correlator in the environment does not have any time ordering. We can simply relate it to one of the Wightman functions using translational invariance of the thermal bath, which can be related to the other Wightman function as well as the spectral function. For our bosonic environment operators (4.15), where n B is the Bose-Einstein distribution and ρ AB (k) is the spectral function. Details for the weak coupling computation of the Wightman functions are given in Appendix A. Now we compute the matrix element of the subsystem operator: Since we focus on one particle quark state with a collinear momentum, we can use the collinear momentum to label the states. As discussed in the unitary evolution, spin does not flip in the transition and we will average over the initial state colors and sum over the final state colors. We choose |a = |p 1 , |b = |p 2 , |c = |p 3 and |d = |p 4 . Then at leading order we have, We can now separately evaluate each term in our expression for non-unitary evolution. We first compute we can use this to eliminate the integral over one of the momenta. The integrals over x 1 and x 2 then set p 2 = p 4 . We further obtain E p 2 = E p 4 and E p 1 = E p 3 because of the on-shell particles. Then we can apply the same trick to the time integrals as we did for the unitary evolution in the previous subsection. The two time integrals will lead to one delta function for energy conservation, multiplied by the time length t in the Markovian approximation. When the dust settles, we are left with where E p = |p| and E q = |q|. The dummy color indexes A and B are summed over implicitly. This term appears as in the Lindblad equation. When the initial subsystem density matrix is |Q 0 Q 0 |, it gives back a projection onto |Q 0 Q 0 |. Then the anti-commutator gives the same result, which cancels the factor of one half. Due to the negative sign, this term in the Lindblad equation represents the loss in the probability of staying in the state |Q 0 Q 0 |.
We have one more term in the Lindblad equation: where the Lindblad operator is defined by L ij = |p i p j |. In the Markovian approximation Further simplification depends on the initial subsystem density matrix ρ S (0) and the final measurement operator applied.

Transverse Momentum Broadening
We now apply the Lindblad equation to compute specific observables on the subsystem density matrix. For an observable associated with a subsystem operator M , the measurement result at time t is defined by the expectation value In our case, the Fock state of the subsystem consists of only one particle states, labelled by a collinear momentum. So we can write In the future, when we include splitting in the study, we will need to include all possible states in the Fock space such as a two particle state. An interesting observable is the final transverse momentum distribution of particles within a jet (with respect to the jet axis). If we focus on the leading parton in the jet, the final transverse momentum vanishes in vacuum. But in the medium, due to the interactions with the medium, the leading parton can develop a non-zero final transverse momentum. In other words, the transverse momentum distribution is broadened by the interactions with the medium. To obtain the final momentum distribution, we set the measurement operator to be a projection operator P Q = |Q Q| for some momentum of Q. Then we need to compute We will derive an evolution equation for the observable P Q (t) that will correspond to a linear differential equation in time. By solving it, we can resum various effects. The resummation can be done given that the Markovian approximation holds. To obtain the measurement results, our job now is to compute the time evolution of the diagonal piece of the subsystem density matrix. 5 Rigorously speaking, one must first show the two frequencies associated with time in the exponents are equal and then apply Eq. (4.29) to obtain the energy conservation delta function multiplied by the time length. Here we write down two delta functions for energy conservation without the factor of time length t. Later we will show the two energy conservation functions are the same for our case here, which allows us to write one of them as the time length.
We first put all non-vanishing leading order pieces derived in the previous section together where we neglect the medium-induced unitary evolution since later we will assume the initial density matrix is a single parton with a given momentum. Sandwiching between Q| and |Q leads to which schematically can be written as where we defined the dissipation rate R(Q) and the fluctuation kernel K(Q, q) (5.8)

Markovian Approximation
To convert Eq. (5.6) into a different equation in time, we will use the Markovian approximation again. We will move the first term on the right hand side to the left hand side, divide the equation by t and take the limit t → 0. Then we infer the master equation for the probability of being in a specific momentum state P (Q, t) ≡ Q|ρ S (t)|Q One should be cautious here because previously we have taken t → ∞ when computing the time integral to obtain delta functions for energy conservation. These two seemingly contradictory limits are compatible with each other in the Markovian limit. The Markovian approximation is valid if the environment correlation time is much smaller than the subsystem relaxation time. When we take t → 0, we are thinking of the time length as a typical subsystem relaxation time. This time length is still much larger than the environment correlation time. What seems to be a short time for the subsystem is actually very long for the environment. The master equation (5.9) is coarse-grained. Physically, the environment has lost any information about the subsystem before it interacts again with it. In this way, at each interaction point, the environment has no memory about the past history of the subsystem. Using Eq. (5.9) derived above, we can check the validity of the Markovian approximation. The dissipation rate R(Q) is associated with a typical time scale of the subsystem relaxation 1/R(Q), which is originated from the Glauber exchange between a collinear parton and a soft parton from the medium. On the other hand, the typical time scale for the environment decoherence is of the order of 1/T , with T being the temperature of the thermal bath. So for the validity of the Markovian approximation, we require .
From Eq. (5.41), which will be explained in Section 5.4, we can estimate the relaxation rate as so the Markovian approximation is valid in the weak coupling limit, where we have used the fact that m D ∼ gT . The structure of the master equation (5.9) is simple: The first term on the right hand side is a loss term for the state |Q . The probability of being in the state |Q decreases with time because it may transition to other momentum states due to the Glauber exchange with the medium. The last term is a gain term for the state |Q . It originates again from the Glauber exchange. States with other momenta, say q, can turn into the state |Q by exchanging momentum with the medium.

Solution to Master Equation
Before we show the solution to the master equation, we want to elucidate our notations. We will use r ⊥ and k ⊥ to label the Minkowski transverse vectors while r ⊥ and k ⊥ to label the Euclidean transverse vectors. For the magnitude, we will use notations such as |r ⊥ | and |k ⊥ |.
We can rewrite the master equation (5.9) in a suggestive form (see Appendix B for the explanation) where we abused the notation: K(Q, k ⊥ ) here includes both K(Q, q) in Eq. (5.8) and some integration from the last term of Eq. (5.6). Due to the expansion based on our power counting, we see that the second term is a convolution only in the ⊥ direction. Hence the obvious way to solve this equation is to move to the space of impact parameter. Defining we can simplify Eq. (5.12) as We now obtain a simpler equation which suggests a solution of the form Then we can write out the final solution of our measurement results at time t as We will examine the distribution in transverse momentum, starting out with a collinear quark that has zero transverse momentum. Then our initial condition for the subsystem density matrix is 20) in which f (Q − ) is the overall normalization of the initial density. If we just focus on the transverse momentum distribution, we can set f (Q − ) = (2π) 2 . Our solution then becomes For the distribution at any time t, we only need to evaluate the T and Q ⊥ dependent Sudakov factor S(Q, r ⊥ ) ≡ −R(Q) + K(Q, −r ⊥ ). Using the results from Appendix B, we find for a collinear quark scattered off soft quarks of the medium 2Nc , T F = 1 2 and N f is the number of active quark flavors in the medium. The angular integration over φ k can be done by using the Bessel function of the first kind: Then the Sudakov factor can be written as The integrand over p − is from 0 to ∞. The integrand is regular as p − → ∞ or |p ⊥ | → ∞ due to the Fermi-Dirac distribution. On the IR side, the seemingly singular point p − = 0 is actually regular because the integrand scales as (p − ) −2 e −|p ⊥ | 2 /(2p − ) when p − → 0 for non-vanishing |p ⊥ |. If |p ⊥ | = 0 the integrand is vanishing. IR singularity exists when |k ⊥ | → 0. In this limit, the Bessel function behaves as So it partially cancels out the singularity of 1/|k ⊥ | 3 at k ⊥ = 0. But the Sudakov factor still has a logarithmic singularity. We will discuss this singularity in detail in Section 5.3. In our numerical studies shown in Section 5.4, we will cut this IR divergence by introducing the Debye screening. Another singular behavior can appear in the final solution because our initial density is a delta function in the transverse momentum. To make the initial delta function more explicit in the final solution, we will reorganize our result as follows: The physical meaning of the separation is as follows: The first term describes that the density of the initial state (Q ⊥ = 0) decays over time with the rate R(Q). The second term is the growing of the density of other states (Q ⊥ = 0). If we integrate over Q ⊥ , we will find the total probability is conserved which is to say that the time evolution preserves the trace of our density matrix.

IR Safety
We are interested in the physical regime Q ⊥ ∼ T to which only the second term in Eq. (5.26) contributes. We define If we expand out the exponent of the right hand side, we have where the rate and the kernel have the form of The function W(k ⊥ ) can be read from Eq. (5.24) and does not have any singularity at k ⊥ = 0. In fact, the function W(k ⊥ ) goes to a constant value at both low and high k ⊥ (see Fig. 1). For n = 0 and 1, we can work out the result easily For n = 0 and 1, the singularity at k ⊥ = 0 does not influence the computation of the integral. The IR singularity is cut-off by the constraint Q ⊥ in the final measurements. We can work out the n = 2 term in the series to see a non-trivial cancellation of some IR divergences: We have two regions of manifest IR divergence here: k ⊥ → 0 and k ⊥ → −Q ⊥ . The first term in Eq. (5.33) is symmetric under the interchange k ⊥ ↔ −Q ⊥ − k ⊥ . So the two regions have the same singular behavior. Expanding about these two singular regions (we only need to expand around one of them and then multiply by two), we obtain the leading singularity (LS): Given the form of R(Q) in Eq. (5.29), we see that the leading IR singularity cancels out. However, we can still have subleading IR singularities that do not cancel. To see this more explicitly, we can set W(k ⊥ ) to be a constant since it has a very mild dependence on k ⊥ . We can then write We can then expand in terms of small |k ⊥ |. As before, we have exploited the symmetry of k ⊥ ↔ −Q ⊥ − k ⊥ to account for both singular regions. After expanding, we can do the angular integral to obtain As we can see, the leading singularity cancels out, but a subleading logarithmic singularity remains. The final measurement constrain Q ⊥ can only cut-off the leading IR singularity but not the subleading one. This is connected with the Markovian approximation, as will be explained in the following. The W term is proportional to α 2 s . So the singularity first appears at O(α 4 s ). One might think that perhaps we are missing some pieces and if we include terms at higher order in the coupling constant, when deriving our master equation, this subleading singularity would be cured. But we can immediately see a higher order term would only contribute to O(α 6 s ) at O(t 2 ) and hence cannot cancel out our IR divergence at O(α 4 s ). This singularity must therefore have a physical origin.
The key point here is that we are working in the Markovian approximation in which all coherence is lost between successive interactions with the medium. Therefore we can treat successive interactions as products of independent scatterings. The only constraint that we are imposing is that the total transverse momentum accumulated by the collinear parton should be of the order of Q ⊥ . The value of Q ⊥ is determined by the final measurement we are conducting, and is assumed to be ∼ T . This implies that if the jet interacts more than once with the medium, it is possible for the jet to accumulate almost Q ⊥ transverse momentum in one interaction and very little in all the others. All the other interactions are therefore independent scatterings in which little or no transverse momentum is exchanged. For these scatterings, our power counting does not apply. We must instead use the EFT in the regime with λ ∼ m D /Q, making it sensitive to the scale m D . Therefore we have to put in a mass regulator such as the Debye mass m D to obtain a finite sensible result. The Debye mass defines the interaction range of these scatterings with little transverse momentum transfer.
This also suggests that in the Markovian approximation, we need smooth transition through all the EFTs from the scale m D up to the scale Q ⊥ ∼ T . The leading singularities discussed above cancel out at all orders in the expansion of Eq. (5.28). We sketch a proof for this in Appendix C. If we work to higher orders in the expansion, we observe that some of the subleading singularities also cancel out. So we conjecture that at the n-th order, we are only left with a ln n−1 singularity. If this is correct, then our solution is resumming a logarithmic series.

Numerical Results
We will calculate the density of states with Q ⊥ = 0. More specifically, we will compute Eq. (5.27) numerically. We will cut the IR divergence by introducing a gluon mass at finite temperature, which is the Debye screening mass in which C A = N c = 3 and N f = 3 is the number of active quark flavors in the QGP (we will assume the strange quark is massless for simplicity). In the following calculations, we will replace |k ⊥ | 4 in the denominator in Eq. (5.24) with (|k ⊥ | 2 + m 2 D ) 2 . The results shown in (5.27) depend on time t, the temperature T of the QGP and the transverse momentum Q ⊥ of interest. We will do the following scalinĝ so that the results at different temperatures fall onto a "universal" curve: TheĜ function is given bŷ in which the rate and kernel are given bŷ Heren F (x) = (e x + 1) −1 . TheĜ function seems to be independent of the temperature T . But it still depends on T through the coupling constant α s . The determination of α s relies on a scale. If we scale everything by T , then α s will depend on T .
The numerical results ofŴ andĜ with constant α s = 0.3 are plotted in Fig. 1. Running coupling effect will be studied in the future. The functionŴ does not vary significantly as |k ⊥ | changes. For |k ⊥ | > 10,Ŵ is almost a constant. In the plot ofĜ, we choose two different Debye screening masses to demonstrate the sensitivity to the IR scale: One is Eq. (5.45) The value of the Debye mass has a significant effect on the broadening rate of the jet. For the smaller Debye mass, the distribution is gradually broadened from the initial peak at the origin. For the larger Debye mass, the shape ofĜ saturates fast and the only change is the normalization, which is given by Eq. (5.45). This is intuitively clear from the fact that for a value ofm D much greater than |Q ⊥ | , the |Q ⊥ | scale becomes irrelevant for the IR physics. Then the shape of the curve is fixed and its amplitude merely scales with time.
We also notice that when |Q ⊥ | m D , the results ofĜ are less sensitive to the value of the Debye mass, where the Debye mass becomes irrelevant. Our simple calculations confirm the physical picture that the transverse momentum distribution of a collinear parton broadens as the parton traverses the QGP. This is still only a partial picture since we have not included other effects such as vacuum parton shower evolution, medium-induced splitting and non-zero chemical potential. So a comparison with data will only be meaningful when these corrections are included in our evolution equation.

Comparison with Previous Work
The solution to the Markovian master equation (5.21) agrees formally with Eq. (2.15) of Ref. [91]. To match the normalization used in Ref. [91], we need to choose f (Q − ) = (2π) 2 . The exponent structure in the solution comes from resumming length enhanced diagrams in Ref. [91] while in our case it shows up after solving the master equation. The solution  resums multiple scattering, as also explained in Ref. [91]. More specifically, the multiple scattering is incoherent. What Ref. [91] calls W R in Ref. [91] exactly matches the structure of the Lindblad equation.
Next we compare the transverse momentum distribution after single scattering in our study and those calculated in Refs. [91,92]. This quantity is defined as P single in Ref. [92]. In our notation, this quantity is given by Our EFT has been designed with an expansion in the small parameter Q ⊥ /Q along with the condition Q ⊥ m D . It is therefore valid for a wide range in the temperature T . It is possible to obtain simple analytic expressions in certain limiting hierarchies. To compare with earlier literature, we will focus on two kinematic regions: |Q ⊥ | T and T |Q ⊥ | m D . Our EFT framework is still valid in these two limiting cases (the only difference is that the power counting parameter λ is now given by Q ⊥ /Q − rather than T /Q − ). For Q ⊥ = 0, neglecting the Debye mass, we find where 4πα s = g 2 (note the difference in the definition of W in Eq. (5.30) andŴ in Eq. (5.42)). Since here we only focus on the Glauber scattering between a collinear quark and a soft quark, we will only compare the results relevant for quarks. The results of Ref. [92] for quark-quark scattering in these two regions can be written as Our result (5.47) agrees with Ref. [92] in these two limits. This can be easily checked using the numerical result ofŴ shown in Fig. 1 or referring to the limiting formulas (A. 19, A.20) in the Appendix A.

Conclusions
In this paper, we take the first step towards developing an effective field theory description for energetic jets traversing a region of the QGP medium. We treat the jet as an open quantum system interacting with a QGP environment in thermal equilibrium at constant temperature T . For now we restrict ourselves to the case when the scale T is perturbative, which is valid at high temperature. We focus on the transverse momentum broadening in this paper and neglect parton splitting. The interaction between the jet and the medium is encoded in effective operators mediated by a Glauber momentum exchange, which captures the dominant forward scattering regime. Tracing out the degrees of freedom of the environment yields a Lindblad evolution equation for the density matrix of the jet. The Lindblad equation turns into a master equation in the Markovian limit if the initial jet density matrix is diagonal in the transverse momentum space and a final projective measurement onto a specific transverse momentum Q ⊥ is imposed. Given that the Glauber mode connects the subsystem and the thermal QGP only via its transverse momentum, we observe that we can analytically solve the time evolution of the master equation by going to the impact parameter space. We work in the scale hierarchy Q ⊥ m D Λ QCD . The EFT works for a wide range of high temperatures and the expansion parameter for our EFT description is set by λ ∼ Q ⊥ /Q 1. We demonstrate that the leading singularities cancel out but subleading singularites still exist, which we conjecture are powers of logarithms to all orders. We surmise that the residual IR singularity is a consequence of the Markovian approximation used to describe the density matrix evolution. We calculate numerically the transverse momentum distribution by cutting off the IR singularity with a gluon Debye mass. The distribution becomes broadened as time increases. Our results agree with those previously derived in literature in various limiting regimes of our EFT.
Our ultimate goal is to understand the jet quenching inside a dynamically evolving QGP in a theoretically controlled way. Looking forward, there are several questions that we would like to address using our current approach. The most urgent one is to incorporate the effects of vacuum shower and other medium-induced effects (such as the medium-induced splitting) systematically. Furthermore, non-Markovian effect of the jet dynamics inside the QGP (such as the LPM effect 6 ) should be investigated in our framework. We can also apply our formalism to study a heavy quark jet traveling through the medium. For the heavy quark jet, we need to first construct effective field theory for the forward scattering of a boosted heavy quark and then use those operators in our open quantum system formalism. 6 At least three time scales are involved in the LPM effect: the environment correlation time τE ∼ 1/T , the Glauber exchange time scale τG ∼ 1/(αsT ) and the formation time of the radiated gluon τF ∼ where xE denotes the energy of the radiated gluon. The LPM effect is important when τF τG, when multiple Glauber exchanges can happen during the formation of the radiated gluon. These Glauber exchanges have to be resummed coherently in the amplitude level and their interference will lead to a suppression of the radiation. Thus, these Glauber exchanges cannot be treated as independent scatterings, and cannot be resummed in the Markovian master equation. One may construct effective operators in the Hamiltonian for the LPM modified radiation. If τF τE is valid, then each LPM modified gluon radiation can still be treated as independent processes, i.e., Markovian.
Finally, we should explore how to define and compute jet substructure observables in our formalism and how to incorporate the wake of a jet into the formalism. Finally, we want to understand whether universal nonperturbative effects can be suitably parametrized in our formalism and then extracted from experiment. Understanding these questions will help us make better use of jets as probes of the QGP.

Acknowledgments
We thank Krishna Rajagopal, Iain Stewart and Yi Yin for inspiring discussions. This work is supported by the Office of Nuclear Physics of the U.S. Department of Energy under Contract DE-SC0011090 and Department of Physics, Massachusetts Institute of Technology.

A Wightman Functions in Thermal Bath
We will use the imaginary time formalism to compute the finite temperature correlation function of the soft operators. The operator for the soft quark current is given by (3.14): where the soft quark operator ψ n s is dressed with a soft Wilson line ψ n s = S † n ψ s . In our weak coupling calculation, the soft Wilson line can be dropped at leading order. In the imaginary time formalism, we first compute the correlator in the Euclidean space and then analytically continue to the Minkowski space. The Euclidean correlator in momentum space is defined by where β = 1/T . The Euclidean coordinate is X = (τ = it, x), the Euclidean momentum is K = (k , k) and K · X = k τ − k · x. The Matsubara frequency for bosonic operator (our soft quark current operator is bosonic even though it comprises quark fields) is k = 2 πT where is an integer. Plugging the soft quark current into the correlator leads to where P = (p n , p) and Q = (q m , q). Here p n = (2n + 1)πT and q m = (2m + 1)πT are the Matsubara frequencies for fermionic operators. The notations here are Euclidean: The overall minus sign comes from the fermion loop. After the integral over the Euclidean coordinates, we obtain We now apply the trick known as the Saclay method to sum over the Matsubara frequencies.
To that end, we introduce a Kronecker delta function by writing , (A.5) where the function f contains all the trace factors in the numerator, with E 1 = |p|, E 2 = |p + k|. Moreover, p n and q m are fermionic Matsubara frequencies. Then we can write (A.6) We can further simplify the above expression by using the following relations: where n F (E) = (e βE + 1) −1 is the Fermi-Dirac distribution. After applying the summation formulas, the integral over τ can be easily done.
Once we obtain the correlator in the imaginary time formalism, we can obtain all the real time Green's functions and the spectral function via analytic continuation. First, the spectral function can be obtained by In our case, Eq. (A.8) leads to The four terms here present four different scattering processes respectively. The ones of our interest are the second and third term, which correspond to scattering processes with one incoming and one outgoing soft quark. The other two terms correspond to processes with either two incoming soft quarks or two outgoing soft quarks. In fact, these two terms (the term proportional to δ(k 0 + E 1 + E 2 ) or δ(k 0 − E 1 − E 2 )) do not contribute at the order we are working here, since the gluon exchanged in these processes scales as a soft mode rather than a Glauber mode. In our power counting, the soft mode scales as p s ∼ Q(λ, λ, λ) and the sum of any two soft modes scales similarly. If the exchanged gluon is soft, the collinear particle will become off-shell. So this process is suppressed and we can drop them. At the same time, we can use our power counting to drop the term n · k (which is Glauber and scales as ∼ λ 2 ) when compared withn · p or n · p (which are soft and scale as ∼ λ) If we define ρ AB (k) = ρ(k)δ AB , we can write the Wightman function D AB Plugging the spectral function gives where n F appears for the initial state while 1 − n F shows up in the final state. This is the standard Pauli blocking in quantum statistics. In order to apply the power counting in k, we need to express everything in light-cone coordinates. To that end, we make a change of variables p → −p in the second term and introduce a dummy four-momentum variable q to write The integral over q gives (A.14) where we define the integral as I. Now our task is to simplify the integral I(k − , k ⊥ ) ≡ d 4 p δ + (p 2 )δ + (k + p) 2 2(n · p) 2 n F (p 0 ) 1 − n F (k 0 + p 0 ) = |p ⊥ | d|p ⊥ | dp − dp + dφδ(p − p + − |p ⊥ | 2 )δ (p − + k − )p + − |p ⊥ + k ⊥ | 2 (p + ) 2 where we have dropped k + according to our power counting. For later convenience, we now calculate where p − > 0 and the values of p + and k − are fixed by integrating over the two delta functions: In the limit |k ⊥ | → 0, we find k − → 0, so . (A. 19) In the limit |k ⊥ | → ∞, we find k − → ∞, so we can approximate the (1 − n F ) term by one to obtain So far, we only considered one flavor of massless soft quark. If we assume the medium consists of N f = 3 massless soft quarks (we neglect the strange quark mass), we find after putting everything together

B Rate and Kernel
In this appendix, we will explain the computation of the rate R and the kernel K for a collinear quark scattering with soft quarks of the medium.

B.1 Rate
From Eq. (5.7), we have the expression for the dissipation rate that appear in the final master equation where we restore the spin indexes s, r here and we will show they are equal. To simplify this, we use the fact that D AB > ∝ δ AB and write From now on, we will omit the spin indexes s, r. We are left with which we have used our power counting to expand away any power corrections. Here Q is collinear and k is Glauber. As we have shown in Appendix A, the k + dependence in D > (k) is dropped out since it is subleading in our power counting. Hence, we can easily do the integral over k + which leads to Plugging Eq. (A.21) into R(Q) leads to where the integrand is independent of φ k and the integral over φ k can be done trivially.

B.2 Kernel
For the solution to the master equation, we also need to compute the kernel, which is defined by Eq. (5.14) in the transverse plane. Its expression can be obtained from Eqs. (5.6, 5.8 and 5.14) The evaluation of the kernel is almost the same as the rate, except for the extra phase e −ik ⊥ ·r ⊥ . The result can be written as where φ k is the relative angle between k ⊥ and r ⊥ .

C Proof for Cancellation of Leading IR Singularities
Here we sketch a proof of the statement that the leading singularities cancel out at all orders in the expansion of t, which is the expansion in Eq. (5.28). We will use the method of mathematical induction. The term at n + 1-th order in the expansion can be written as (we will drop the Q dependence in relevant functions) We have now written the result in terms of the n-th order term. We now only need to consider new singularities that arise as we go from the n-th to (n + 1)-th order term. We first check the first term in the expression above: Eq. (C.2) may have two new singularities at leading order: The first is when k ⊥ → 0 and the second is when k ⊥ → −Q ⊥ (see the explicit second order calculation in the main text).
In the first case the exponential is approximately unity and the expression reduces to which cancels out with the third term in the last line of Eq. (C.1). No other new leading order singularities can show up in the k ⊥ → 0 region because all leading order singularities cancel out in the n-th term In the second case when k ⊥ → −Q ⊥ , the term in Eq. (C.2) has no leading singularity since K(−r ⊥ ) − R has only logarithmic singularity, which is subleading. So for the consideration of leading IR singularity, Eq. (C.2) in the second region k ⊥ → −Q ⊥ reduces to − d 2 r ⊥ e −ir ⊥ ·Q ⊥ K(−r ⊥ ) − R n (C.6) which cancels out with the second term in the last line of Eq. (C.1).