Quark branching in QCD matter to any order in opacity beyond the soft gluon emission limit

Cold nuclear matter effects in reactions with nuclei at a future electron-ion collider (EIC) lead to a modification of semi-inclusive hadron production, jet cross sections, and jet substructure when compared to the vacuum. At leading order in the strong coupling, a jet produced at an EIC is initiated as an energetic quark, and the process of this quark splitting into a quark-gluon system underlies experimental observables. The spectrum of gluons associated with the branching of this quark jet is heavily modified by multiple scattering in a medium, allowing jet cross sections and jet substructure to be used as a probe of the medium's properties. We present a formalism that allows us to compute the gluon spectrum of a quark jet to an arbitrary order in opacity, the average number of scatterings in the medium. This calculation goes beyond the simplifying limit in which the gluon radiation is soft and can be interpreted as energy loss of the quark, and it significantly extends previous work which computes the full gluon spectrum only to first order in opacity. The theoretical framework demonstrated here applies equally well to light parton and heavy quark branching, and is easily generalizable to all in-medium splitting processes.


I. INTRODUCTION
The attenuation of the production cross section of energetic particles and jets in high energy reactions with nuclei is one of the primary signatures of inelastic parton scattering in dense nuclear matter [1,2]. The rapid development of heavy ion programs at fixed target and collider experiments fueled tremendous interest in medium-induced bremsstrahlung processes and radiative parton energy loss in Quantum Chromodynamics (QCD) [3], often discussed in analogy with the Landau-Pomeranchuck-Migdal effect for photon emission in Quantum Electrodynamics (QED) [4,5]. To this effect, initial efforts have focused on the energy loss of energetic quarks and gluons as they propagate in the quark-gluon plasma (QGP), a deconfined state of strongly-interacting matter that existed in the early universe and is recreated today in relativistic heavy ion collisions. Radiative energy loss in QCD is synonymous with soft gluon bremsstrahlung, a process in which hard quarks and gluons shed energy in small quanta during propagation through matter. As a result, the leading parton always remains the most energetic. This does not preclude the possibility that it may dissipate a sizable fraction of its energy, but this is achieved through multiple gluon emission. All radiative parton energy loss approaches rely on perturbative techniques and treat the interactions of the jet with the quasi-particles of the medium primarily through t-channel gluon exchanges [3].
Important considerations, where various approaches available in the literature deviate, are the kinematic regimes in which the parton system is produced and the size of the nuclear medium. To give an example, let us denote by l f the quantum-mechanical time for the splitting process or gluon emission to occur, by λ the scattering length of the parent parton, and by L the size of the medium. In practice, it is well understood that none of those quantities are fixed. Scattering lengths depend on the density and dynamics of nuclear matter, and the size of the realistic medium cannot be defined in a model-independent fashion. The time for the splitting process to take place depends sensitively on the kinematics and interactions in the medium. If one considers the simultaneous limit l f /λ ≫ 1 and L/l f ≫ 1, a continuous approximation to the process of parton system propagation can be made, where a Schroedinger-like equation can be solved based on a Gaussian approximation [6,7]. Similarly, this limit can be treated in a path integral approach [8,9] where a Gaussian approximation is again essential to obtain a tractable resummed result. A hard thermal loop approach to the problem of gluon emission has also been developed [10].
The situation where the number of scatterings in an energetic parton system is very large is not often encountered in reactions with heavy nuclei [11]. Even more importantly, the splitting time which is the inverse virtuality of the two parton system before and after the branching 0 < l f < ∞ is different for different diagrams and does not obey any particular hierarchy with respect to fixed length scales. This motivates the development of formalisms where the exact kinematics of parton propagation in matter are kept [12,13]. In this spirit, the result for gluon bremsstrahlung is written as a series in the correlations between multiple scattering centers, or the opacity. In the calculation of more inclusive jet quenching observables, the contribution of higher orders in opacity is power suppressed, which is most easily seen in the high twist approach [14,15]. Differences can be expected, however, for observables that are more sensitive to the double differential splitting kernels, such as the ones related to jet substructure [16,17]. The evaluation of higher orders in opacity can become numerically expensive, even in the simplifying limit of soft gluon emission. It was shown, however, that in certain simplifying geometries corrections up 9 th order in opacity can be evaluated, albeit with decreasing numerical accuracy [18]. It was pointed out that the path integral approach, without the approximations above, is formally equivalent to the opacity series and can be re-expanded into diagrams to recover the opacity expansion [19]. Finite size effects have been introduced in the hard thermal loop approach as well [20].
The popular energy loss approach has been extended to account for the mass of heavy partons [21][22][23][24]. Initial-state inelastic processes have also been discussed [25,26], first motivated by fixed target experiments [27]. Higher order corrections to the jet-medium interactions [28] and soft splitting interference have been studied [29,30]. Last but not least, soft photon emission has been calculated using techniques similar to the ones described above [10,[31][32][33].
Advances in the theoretical understanding and experimental measurements of reconstructed jets necessitate more precise control over in-medium branching processes. This requires the development of theories of hard-probe production in the presence of nuclear matter that transcend the thirty-year-old energy loss approach. An important step in this direction is to obtain in-medium splitting kernels beyond the soft gluon approximation. The corrections to vacuum branching in the higher twist approach were discussed in Ref. [34]. A full set of in-medium splitting kernels, q → qg, g → gg, q → gq g → qq, to first order in opacity was derived in Ref. [35] in an effective theory of jet propagation using Glauber gluon interactions sourced by an external potential [36,37]. The λ ≪ l f ≪ L limit has also been generalized to full longitudinal splitting kinematics for the g → gg process [38]. The energy loss limit has also been relaxed in Ref. [39]. Multiple parton branching beyond the soft gluon approximation has been rigorously calculated in Ref. [40]. To address heavy flavor observables in ultrarelativistic nuclear collisions, the heavy flavor splitting kernels Q → Qg, Q → gQ g → QQ were obtained to first order in opacity [41]. It is important to realize that to bridge the gap between high energy particle physics and high energy nuclear physics and to reduce the systematic uncertainties inherent in phenomenological models of jet quenching, resummation [42][43][44][45] and higher order calculations [41,46] are necessary. These are not only facilitated by, but, in fact, require the implementation of the full in-medium splitting functions.
In this paper, we focus on the calculation of the gluon emission spectrum beyond the soft gluon approximation to an arbitrary order in opacity, with semi-inclusive deep inelastic scattering (SIDIS) processes in mind. We will illustrate the technique on the example of a quark jet. This is the dominant channel at leading order, and in the kinematic regions of interest the next-to-leading order corrections to the flavor composition of jets are small. Our work sets the stage for the calculation of all in-medium splittings to high orders in opacity and for related phenomenology. We write down explicitly the result for the quark to quark + gluon splitting channel, including light and heavy quarks, to second order in opacity and expect that our results will find application in the description of semi-inclusive hadron suppression [42,[47][48][49][50] and further expand the program at a future EIC [51].
Although we approach the problem in the standard framework of Feynman perturbation theory, the Fourier transformation to work with fixed longitudinal positions of the scattering centers naturally brings in many of the intuitive elements of Light-Front Perturbation Theory (LFPT) [52,53]. As such, the parton branching will be formulated in terms of light-front wave functions (LFWF) and the associated energy denominators which describe a fluctuation in the virtuality of an intermediate state; these quantities possess clear and intuitive kinematic dependences because of the Galilean invariance of the light-front Hamiltonian [38,52]. Light-front dynamics is a natural description for many high-energy processes, and in that context, the phases due to multiple scattering which stimulate radiative energy loss are also crucial for the construction of single spin asymmetries in QCD [54][55][56]. Quark jet production at an EIC is also intimately related with the transverse-momentum-dependent parton distribution functions of a heavy nucleus, which have been previously explored [57,58].
Our work is organized as follows: in Sec. II we set up the calculation of the gluon-in-quark-jet distribution in SIDIS, defining our frame and coordinates and laying out the origin of medium modifications to this distribution. In Sec. III we formulate the coupling of the jet to the medium at the Lagrangian level, using this to compute both the vacuum distribution and the modification from the medium at first order in opacity. The results obtained here utilize the language of LFPT, reproducing and generalizing some previous results in the literature. In Sec. IV we generalize the first order in opacity results to construct the kernel of a recursion relation, which is used to derive expressions for the gluon-in-quark-jet distribution at any arbitrary order in opacity. To illustrate the use of this recursion relation, in Sec. V we explicitly compute for the first time the exact distribution to second order in opacity, verifying that it reproduces the known result in the literature in the soft-gluon approximation. The primary new results of this work are the recursion relation Eq. (73) with the accompanying kernel (the "reaction operator") Eq. (74) and the explicit second order in opacity results Eqs. (76) and (77). Finally, in Sec. VI we conclude by summarizing our key results and outlining the applications which we intend to pursue in future work.

II. COLD NUCLEAR MATTER EFFECTS IN DEEP INELASTIC SCATTERING
Consider the process of jet production in Semi-Inclusive Deep Inelastic Scattering (SIDIS) on a heavy nucleus, γ * + A → (jet) + X. We will work in the Breit frame, in which the virtual photon travels along the +z axis with a large momentum p + and the nucleus travels along the −z axis with a large momentum p − N per nucleon. Here and throughout this paper, we employ light-front coordinates v ± ≡ g +− 2 (v 0 ± v 3 ) and denote transverse vectors as v ≡ (v 1 ⊥ , v 2 ⊥ ) with magnitudes v T ≡ |v|. Different references adopt different conventions for the normalization of the light-front coordinates; here we will choose g +− = 1, although it is not uncommon to see the convention g +− = 2 in the literature. In this frame, the momentum of the virtual photon and struck nucleon are with the Bjorken variable The photon-nucleon center-of-mass energy √ s is given by s = and for power counting purposes we will consider p + ∼ p − N ∼ O (Q); that is, a situation close to the center-of-mass frame. 1 We will also work in the light-cone gauge A + = 0 for which only physical gluon polarizations exist; this gauge is also convenient because it is equivalent to Feynman gauge with eikonal accuracy. This setup is illustrated in Fig. 1. At leading order in the strong coupling, the virtual photon strikes a quark inside one of the nucleons, ejecting it with the large longitudinal momentum p + along the +z axis. This high-energy quark propagates through the rest of the medium, potentially undergoing multiple scattering in the process, before emerging in the final state as a quark jet with momentum p. For a simple process like SIDIS, factorization theorems have been proven which express the inclusive jet production cross-section p + dσ d 2 p dp + in terms of parton distribution functions (PDFs): TMD factorization [59,60] at low p T which is sensitive to the nonperturbative transverse momentum distribution within the nucleus, and collinear factorization (see e.g. [61]) at high p T which generates this momentum by a perturbative hard recoil. Note that the jet transverse momentum p T has a significantly different interpretation for jets produced by SIDIS in the Breit frame than for heavy-ion collisions. In the latter, the jet p T is a hard scale measuring the total energy of the jet, whereas in the former, the jet energy corresponds to p + and flows primarily in the longitudinal direction, with p T being a potentially soft scale.
In this work, we will focus on the effects of multiple scattering on the jet and on its substructure, so it is sufficient to regard the interaction of the virtual photon with the nucleus as a source current j which creates a high-energy quark jet at some point in the medium. For SIDIS, this current can be expressed directly in terms of the PDFs, but this abstraction also allows us to apply the general discussion of medium modification to other processes, such as heavy-ion collisions, in which the medium through which the jet propagates can be very different from the cold nuclear matter present in SIDIS.
Measurements of jet substructure [62][63][64][65][66] are more sensitive to the detailed QCD dynamics of in-medium parton branching than the overall inclusive or tagged jet production cross-sections, though the latter are larger and more accurately measured at present [67][68][69][70][71][72][73]. At leading order in the strong coupling, the first contribution to the substructure of the quark jet comes from its splitting into a quark-gluon system. Thus, we will consider in particular the distribution of gluons within the quark jet with k µ the momentum of the gluon, p µ the total momentum of the quark jet, and x ≡ k + /p + the longitudinal momentum fraction of the gluon within the jet. At leading twist in the 1/Q expansion, the gluon will be emitted as in the vacuum, without modification from the medium. Such splittings, if integrated out, contribute to the leading-logarithmic Dokshitzer-Gribov-Lipatov-Altarelli-Parisi [74] or Collins-Soper-Sterman evolution of the underlying quark PDF [57,59]. For this vacuum substructure of the quark jet to be modified by the medium, the gluon must be emitted from an interaction inside the medium itself. Because the medium is highly Lorentz-contracted in Bjorken kinematics, these modifications are necessarily suppressed by the Lorentz gamma factor and are intrinsically higher-twist effects. Let us, for the sake of discussion, take nuclear matter of fixed size and transport properties. If the length a parton propagates in the medium before escaping is L and its mean free path is λ ≡ 1 ρσ el , then the ratio of these quantities χ ≡ L/λ is referred to as the opacity: the average number of scatterings which that parton will undergo. (Equivalently, one can boost all these scales from the rest frame of the medium to the Breit frame, writing χ = L + /λ + .) The mechanisms by which scattering in the medium will modify the jet substructure are generally classified as collisional versus radiative energy loss. Collisional energy loss refers to the higher-twist corrections which generate a direct transfer of energy between the jet and the medium; these effects, and other higher-twist corrections, generically enter at O ⊥ 2 Q 2 with ⊥ 2 some associated transverse momentum scale. Radiative energy loss refers to parton branching and the stimulation of gluon emission, which is associated with the accumulation of phases as the jet propagates through the medium. As we will see, a typical phase ∆φ is given by the product of a wave number 1 l + f ∼ ⊥ 2 p + and the mean free path λ + between scatterings: with A 1/3 the number of nucleons (or scattering centers, more generally) which is proportional to the length of the medium. If the opacity χ is not extremely large (χ ∼ A 1/3 such that every available nucleon / scattering center is struck) then the phases associated with radiative energy loss are enhanced by the length of the medium and dominate over collisional energy loss and other higher-twist effects. Finally, the average number of scatterings which takes place in the medium is controlled by the opacity χ, with the inclusion of correlations between multiple scatterings in the medium contributing to higher powers of the opacity. A perturbative approach to the computation of medium modification thus generically translates the perturbation series in the strong coupling α s , here denoting the interaction between the jet and the medium, into a series in the opacity χ for each correlated rescattering, referred to as the opacity expansion. At relatively small values of χ, the opacity expansion can be truncated at finite order, while for very large values it must be resummed. The precise value of the opacity at which resummation is mandatory depends on the observable, but in general it need not be restricted to χ ≪ 1. Indeed, for radiative energy loss in the soft-gluon approximation, rapid convergence is seen up to O χ 3 even for opacities as large as χ ∼ 5 [13].
In this work, we will extend the opacity expansion approach to relax the soft-gluon approximation up to any finite order in opacity; this will yield the full solution for gluon radiation from both light and heavy quarks to any desired accuracy. For practical purposes, we implicitly assume that we have the ability to truncate the opacity expansion, such that χ = L + λ + ∼ O (few). Together with the assumption that the accumulated phases associated with gluon radiation may be large ∆φ ∼ λ + l + f ∼ O (1) but that collisional losses and other higher-twist effects are small ⊥ 2 this defines the parametric regime that motivates our calculation: The leftmost inequality enforces the eikonal approximation (twist expansion) ⊥ 2 Q 2 ≪ 1; the middle comparison permits the accumulated phases to be large ∆φ ∼ O (1); and the rightmost comparison reflects the ability to terminate the opacity expansion 1 < χ < few. For comparison, the continuous Schroedinger-like description is valid in the limit of small phases and high opacities, 1 [6,7]. Because we assume no particular strong hierarchy among the scales l + f , λ + , L + , our results should be general enough to reproduce the other known limits under the appropriate assumptions.

A. Scattering in an External Potential
The coupling of a jet to the external vector potential A ext sourced by the constituents of a nuclear medium can be introduced directly at the Lagrangian level [35]. To linear order in A ext and the coupling of the jet to the medium g ef f we write The added terms introduce ordinary quark/gluon Feynman rules for single scatterings in an external field, and we keep the effective coupling g ef f a free parameter (not necessarily fixed to be the same QCD coupling g in L QCD ). Jets arise from hard scattering processes, the Feynman rules for which are contained in L QCD ; still it is useful to consider for our case a quark source term j with dimensions of m 5/2 and L source =ψj +jψ.
The external field A µ a ext (x) is a superposition of the color fields a µ a i (x − x i ) of a large number N of scattering centers distributed throughout the medium with spacetime positions {x i }: The external potential in the A + = 0 gauge from a single scattering center is given in the eikonal approximation by with (t a ) i a matrix in the color space of the i th scattering center. The scattering potential v(q 2 T ) is proportional to the propagator of the exchanged gluon, which may acquire a thermal mass µ in a hot medium: and its square is proportional to the elastic scattering cross-section of a quark on the scattering center: Note that the color factors given here treat the scattering centers in the medium as being the fundamental representation. If the scattering centers are in the adjoint representation C F /2N c → 1/2. These external potentials will need to be averaged in the medium; for this, we will use Gaussian averaging as described below. This assumption corresponds to limiting the interaction with a given scattering center to two gluons each, and as such all target fields are averaged pairwise. For scattering centers which are locally color neutral, such as nucleons in cold nuclear matter, two-gluon exchange is the leading interaction, with higher-order field correlations suppressed by at least a factor of α s . Using Eqs. (7) and (8), the average of two color fields is The color averaging of the scattering centers yields zero unless i = j: giving Following Ref. [13], we convert the remaining sum over sites into an average, i (· · · ) = N · · · , which we evaluate as an ensemble average over the parameters of the medium. In this case, that corresponds to averaging over the transverse position x i and the longitudinal position x + i of the scattering center: where we take the medium to have transverse area A ⊥ and longitudinal extent 2R + centered around zero, and in the last step we have introduced the mean free path λ = 1 ρσ el with ρ the number density of scattering centers and σ el given by the integral of Eq. (10). Using Eq. (14) one readily obtains (c.f. Eq. (2.2) of Ref. [38]) with We have written Eq. (16) as if the medium had a uniform density, such that ρ and hence λ + were constants, resulting in a correlation γ which depends only on x − y. For more realistic phenomenology, one can generalize this straightforwardly to make these quantities dependent upon the longitudinal position x + or the impact parameter x+y 2 ; this latter generalization is more difficult with path-integral techniques than with the opacity formalism presented here [38].
With the Lagrangian Eq. (5) and the Gaussian averaging Eq. (15) it is straightforward to calculate scattering amplitudes in perturbation theory. For inclusive quark jet production from a source j the amplitude is As illustrated in Fig. 2, this denotes a source at point x 0 which is used to create a quark at point x ′ 0 , which then evolves until a quark q(p) is detected in the final state. The state |Ω denotes the interacting vacuum of the Lagrangian Eq. (5). The amplitude Eq. (17) is then straightforward to calculate using the standard techniques of Lehmann-Symanzik-Zimmermann (LSZ) reduction [75] and the Gell-Mann-Low theorem [76]. At 0 th order in g ef f we have the vacuum amplitude, unmodified by rescattering in the medium, with iS F the free fermion propagator and we note that, since j(∆x) has dimensions of m 5/2 , its Fourier transform j(p) has dimensions of m −3/2 . Squaring the amplitude Eq. (18), averaging over the quantum numbers, and integrating over the phase space of the final state gives Given that j(p) has dimensions of m −3/2 , we see that N 0 has dimensions of m 0 and can be interpreted as the number of produced particles (quark jets) in the final state. Equivalently, this gives the phase space distribution of quarks generated by the source current as 3. Diagrammatic representation of the radiation amplitude Eq. (21). A source at point x0 creates an energetic quark, which propagates through the external field Eq. (8), radiating a gluon in the process.
Similarly, we can use the Lagrangian Eq. (5) to compute the first contribution to jet substructure: the distribution of gluons within the quark jet. For this, we want a similar scattering amplitude to Eq. (17), but with a quark and a gluon in the final state (see Fig. 3): Again, at 0 th order in g ef f we have the unmodified vacuum spectrum of gluon radiation: We can perform the dℓ − integral by residues, closing the contour below for z + − x + 0 > 0 and setting ℓ − = p − − iǫ: where we have replaced the numerator of the quark propagator by a sum (implied) over spinors / ℓ → / p = U (p)U (p) after putting the quark momentum fully on shell. At this stage, all propagators have been put on mass shell, and we have an explicit dependence on the "time" coordinate z + ; this situation is equivalent to formulating the scattering amplitude directly in LFPT. As is natural in any time-ordered perturbation theory, the emission time z + is bounded causally from below by the source time x + 0 , but it is unbounded above. Note that the iǫ regulator from the denominator of the quark propagator automatically regulates the convergence of the integral at the upper limit z + → +∞. Integrating over the emission time z + yields and we recognize the combined factor in parentheses as the light-front wave function (LFWF) for a quark to split into a quark + gluon system, normalized in the conventions of Ref. [77] (although note that there, they use the metric g +− = 2 rather than g +− = 1 as used here) and depicted in Fig Here, x = k + p + is the longitudinal momentum fraction of the radiated gluon, λ is its spin, m is the mass of the quark, and it is convenient to express the dependence on the spins σ, σ ′ of the outgoing and incoming quarks, respectively, in terms of the unit and Pauli matrices [1], [ τ ] [57]. Note that the LFWF depends on transverse momentum only through the intrinsic transverse momentum of the splitting κ = k − xp. In a frame in which p = 0, the transverse momentum of the splitting is κ = k, but after a transverse boost to a frame in which p = 0, this intrinsic momentum becomes κ = k − xp. This characteristic form of the transverse momentum dependence in a LFWF is the origin of the "mysterious" structure in Ref. [7] and reflects the Galilean symmetry discussed in Ref. [38]. We note also that the LFWF's are boost-invariant (since they depend only on x and not on p + ) and explicitly conserve angular momentum Fock state by Fock state (since the orbital factor ǫ * λ · κ is an eigenfunction of longitudinal angular momentumL z with eigenvalue (−λ)). Thus our final expressions are manifestly boost-invariant, even though we started this calculation with the center-of-mass frame in mind.
In terms of the light-front wave function Eq. (25), the vacuum radiation amplitude Eq. (24) takes the simple form Before continuing to compute the gluon phase space distribution, we note the origin of the two phases in Eq. (26). There is a trivial "production phase" e ip·x0 associated with the source at point x 0 emitting a quark with momentum p, and there is an " associated with the interval z + ∈ (x + 0 , ∞) in which the gluon radiation can occur. Note that the energy scale entering the emission phase is exactly the same as the energy denominator entering the LFWF Eq. (25).
Because the upper limit of the emission time is unbounded, the first term in the emission factor is zero; when an in-medium scattering occurs after the splitting, this will introduce an upper bound on the z + interval, and the first term will be nonzero. Squaring Eq. (26) and averaging over the quantum numbers, we obtain Noting that the source current j(p) has dimensions of m −3/2 , the spinors U (p) has dimensions of m +1/2 , and the LFWF in Eq. (25) has dimensions of m −1 , we see that |R 0 | 2 has net dimensions of m −4 . This implies that, after integrating over the on-shell phase space for both the final-state quark and gluon, the result is dimensionless and can be interpreted as a number distribution: such that the phase space distribution of gluons within the quark jet is After summing over the quantum numbers of the outgoing particles and averaging over the spins of the initiating quark, the LFWF's are given by which are just the massive generalizations of the Altarelli-Parisi splitting functions. This gives the vacuum distribution as which, after some algebra, is seen to be equivalent to Eq. (2.29) of Ref. [41]. 2 In the massless limit the connection to the Altarelli-Parisi splitting functions becomes explicit: such that Eq. (30) can be rewritten as B. Leading Order Jet Substructure in the Opacity Expansion Although the distribution Eq. (34) of gluons within a quark jet represents the leading contribution to jet substructure, this contribution consists of gluon emission in the vacuum. The mechanism by which the modification of this substructure occurs is through scattering in the medium. In this process, phases associated with the scattering in the medium accumulate, and the endpoints of the emission phases in Eq. (26) are fixed to be located in the medium. Using the Gaussian averaging Eq. (15), the first contribution from scattering in the medium comes from the 17 diagrams shown in Figs. 5 and 6. If one gluon is exchanged with the jet in the amplitude (and another in the complex-conjugate amplitude), there are 3 2 = 9 "single-Born" or "direct" diagrams D 1−9 as illustrated in Fig. 5. At the same order, there are contributions in which the scattering center exchanges two gluons with the jet on the same side of the cut, generating the 2 × 4 = 8 "double-Born" or "virtual" diagrams V 1−8 shown in Fig. 6.
Single-Born "direct" scattering diagrams on a scattering center at position x + . To illustrate the calculation, consider the three building blocks of the 9 "direct" diagrams at the amplitude level First-order rescattering corrections in the medium. There are three distinct diagrams, in which the rescattering occurs before the splitting (R A 1 , left panel), after the splitting on the quark (R B 1 , center panel), and after the splitting on the gluon (R C 1 , right panel) shown in Fig. 7, starting with the case in which the rescattering happens before the splitting, as shown in the left panel of Fig. 7: whereÂ ≡ A a t a . The manipulations in going from coordinate to momentum space are exact and involve only conservation of momentum. The rest of the calculation proceeds along the lines of the amplitude R 0 for radiation in vacuum: we first perform the dq − integral by residues, enclosing the pole in the upper half-plane only if We again replaced the numerator / ℓ 2 − / q by a spinor sum after putting the propagator on shell. Next we perform the dℓ − 2 integral by residues, picking up the pole at The dz + integral provides the energy denominator to complete the LFWF Eq. (25) and leaves the boundary phases from z ∈ (x + i , +∞), giving where we used the eikonal form Eq. (8) of the external potential. As before, we note that the emission phase can be written in terms of the energy denominator Eq. (27): . We also note the appearance of a new type of phase arising from scattering in the medium: an "impulse phase" e −i[(p−q) − −p − ]x + i which reflects the instantaneous change in the wavelength (inverse formation time) of the quark due to a change in its transverse momentum. After multiplying by one of the complex-conjugate amplitudes, these impulse phases will also cancel or assemble into energy denominators, as we will show.
In the same way, we compute the amplitudes B and C from Fig. 7 to be: Now consider the combination of impulse phases that arises from the interference of any two such amplitudes. Noting that the averaging of the target fields as in Eq. (14) and Eq. (15) yields a delta function setting the momentum transfer q equal in the amplitude and complex conjugate amplitude, we readily see that the impulse phases cancel exactly for the square of any amplitude. For the nontrivial interferences the net effect of the impulse phases can always be expressed in terms of energy denominators in LFPT: This general feature can be understood as follows. The amplitude squared is related by the optical theorem to the imaginary part of a forward scattering amplitude, in which zero net momentum is transferred between the initial state in the amplitude and the "final state" -the initial state of the complex-conjugate amplitude. At this level, the momentum transfer from the medium flows into the jet, through various partons, and back out. Depending on that momentum routing in a given diagram, the momentum flow may pass through the parton branching, modifying the intrinsic transverse momentum κ = k − xp which enters the arguments of the LFWF and energy denominators. Depending on the diagram, the result may be a shift to the jet center of mass p → p − q leading to κ → κ + xq; the relative gluon momentum k → k − q leading to κ → κ − q; or to both, leading to κ → κ − (1 − x)q. This general feature indicates that the kinematic effect of a scattering in the medium can always be incorporated in a shift by one of −q, +xq, or −(1 − x)q in the ingredients of LFPT.
Combining these amplitudes, we form the 9 "direct" diagrams D 1−9 ; together with the two-particle phase space (29), they are: with Note that the terms D 4−9 are interferences of the various single-Born diagrams; as such, they contain nontrivial impulse phases in the bracketed prefactors. Note also that all of the phases here ∆E − z + are formally higher twist as seen in Eq. (27), but that, for relatively small opacities, the differences in positions can be parametrically large: z + − x + 0 ∼ L + , such that these higher-twist effects are length-enhanced as discussed in Sec. II.
Similarly, we can calculate the "double-Born" or "virtual" amplitudes shown in Fig. 8; the interference of these amplitudes with the vacuum amplitude Eq. (26) (along with the overall complex conjugates) generate the diagrams shown in Fig. 6. Note that, in addition to the amplitudes shown in Fig. 8, there are two others which could be drawn: amplitudes in which the scattering center interacts with the jet by exchanging one gluon before the splitting and another gluon after. These diagrams vanish in the eikonal approximation, because it is impossible to simultaneously put the intermediate propagators on shell, leading to poles which cannot be enclosed and hence to zero. Physically, the scattering is instantaneous (within the eikonal approximation) so that there is no time to radiate the gluon during that scattering.
The double-Born contributions are calculated in the same way as the single-Born "direct" contributions, with the only subtlety arising from the treatment of the pole of the intermediate propagator. For concreteness, consider the amplitude R F 2 in which the rescattering takes place on the radiated gluon: In the A + = 0 light-cone gauge, the numerator of the gluon propagator is given by with N µ+ = N +ν = 0 exactly, along with the corresponding component of the polarization vector ǫ + * = 0. Because the eikonal external potential only has nonzero component a − i , this eliminates 2/3 of the terms from the triple-gluon vertices, with only the terms contracting the polarization vector with the gluon numerators remaining: where we have performed the z 1 , z 2 , ℓ 2 , ℓ 3 integrals to enforce momentum conservation. The remaining effect of the numerators N µν is to shift the argument of the polarization vector; using the fact that q + 1 = q + 2 = 0 in the external potential and ǫ + * = 0 in this gauge, we find This evaluation of the numerator algebra also applies to the single-Born amplitude R C 1 Eq. (40) and in general to the rescattering of the radiated gluon. Using the averaging Eq. (15) of the target fields (note here that the second field is not complex-conjugated, leading to a relative minus sign in momentum space) and collecting the poles of ℓ − 1 and q − 1 as usual gives The subtlety associated with the double-Born amplitudes R D 2 , R E 2 , R F 2 resides in the remaining integral over q − 2 : while there is still the pole of the (k − q 2 ) propagator between the two exchanged gluons, there is no longer a Fourier factor regulating the convergence of the integral at infinity. As such, the integral is logarithmic at large |q − 2 | and must be carefully regulated. There are a variety of ways to do this, such as regulating the magnitude of the minus momentum to enforce the eikonal approximation |q − 2 | < p − N or by symmetrizing the integrand under q − 2 → −q − 2 ; either way leads to the result −i/2 for this "contact" limit of the integral: The calculation of the other double-Born amplitudes proceeds along the same lines with amplitudes R D 2 , R E 2 , R F 2 all possessing a "contact" integral, while the amplitude R G 2 has poles that can be collected normally: Again, these amplitudes can be written entirely in terms of the elements of LFPT: the wave functions and the energy denominators. Note the presence of a nontrivial impulse phase for amplitude R G 2 due to a rescattering on the quark-gluon system which redistributes the transverse momentum between the partons, thus altering the virtuality (formation time) of the state.
Combining these amplitudes with the vacuum amplitude Eq. (26) and including the complex conjugates, we form the 8 "virtual" diagrams; together with the two-particle phase space Eq. (29), they are: with As before, we note the presence of nontrivial impulse phases in the bracketed prefactors of V 7−8 , arising this time from double-Born scattering on the quark-gluon system which redistributes transverse momentum between them.
Finally, we combine the direct and virtual contributions, with the various phases adding to their complex conjugates to generate cosines. This gives the first order result in the opacity expansion with exact kinematics: where δz + ≡ z + − x + 0 and L + ≡ R + − x + 0 and we will suppress the explicit x dependence of the wave functions for brevity going forward. Products of wave functions are also implied to be summed / averaged over quantum numbers as in Eq. (31).
This general result at first order in opacity can be compared with a number of others in the literature. Many of these use the "broad source approximation," which assumes that the initial distribution p + dN0 d 2 p dp + is insensitive to shifts p → p − q in the momentum of the initial jet of order q T ∼ O (µ). This assumption corresponds to where we have introduced the constant scale 1 µ 2 for dimensional consistency. 3 Additionally, we introduce the short-hand notation and the gluon mean free path λ + g ≡ CF Nc λ + . (Note also that various references specify p T = 0, which we have not assumed here.) In terms of these quantities, the exact first order in opacity LO result Eq. (55) becomes in the broad source approximation with the explicit form in exact agreement with Eq. (2.51) of Ref. [41]. Having obtained agreement in the massive case, we also subsequently agree with the massless limit of Eq. (9.43) of Ref. [35] and even with the expressions for the single-and double-Born contributions separately in Eqs. (9.15) -(9.18), before the broad source approximation is employed. Finally, in the soft-gluon limit x ≪ 1 and also dropping the leading x 2 m 2 mass dependence we have A = B = k and C = D = (k − q) such that which can be compared with, e.g., Eq. (113) of Ref. [13], again in the broad source approximation. Thus, we see that the general form Eq. (55) contains all of these previous results, encoded into the ingredients of light-front perturbation theory: the light-front wave functions Eq. (25) and Eq. (31) and the energy denominators Eq. (27). The next step will be to construct a general recursion relation (the "reaction operator" [13]) to construct higher orders in opacity in terms of these ingredients.

IV. THE RECURSION RELATION AT FINITE x
Having constructed the first order in opacity contribution to the gluonic substructure of quark jets, we will now proceed to generalize this form to construct higher orders in opacity. The strategy will be to use the 17 diagrams shown in Figs. 5 and 6 as the kernel of a diagrammatic recursion relation (the "reaction operator"), generalizing the expressions in Eqs. (45) and (54). In principle, one may worry that, in going to higher orders in opacity, the color structures associated with multiple scattering may become increasingly complex and not reduce down to a form of the building blocks Eqs. (45) and (54). As a preliminary exercise, let us show that this does not occur, at least at the level of the Gaussian averaging Eq. (15): the color factors remain multiplicative at all orders in opacity.
The color structure for a given order in opacity can always be written in the form  for some color matrix M a , as illustrated in Fig. 9. At the next order in opacity under the approximation of Gaussian averaging Eq. (15), two additional gluons in a color-singlet configuration are inserted into the color trace in one of the 6 topologies shown in Fig. 9. For clarity, here we do not concern ourselves with the placement of the final-state cut or the kinematic factors associated with the scattering amplitude; our goal is to show that the color structure remains multiplicative at any order in opacity, without generating fundamentally new structures. In addition to the 1/N c present in C 0 , the insertion of an extra two-gluon rescattering generates a factor of 1/2N c (from the color averaging Eq. (12) in the medium) and a factor of 2N c /C F (from the conversion to the elastic scattering cross-section Eq. (10)), leading to an overall prefactor of 1/N c C F for the color structures C 1−6 . The evaluation of these color structures is straightforward using the Fierz identity, yielding The particular color factors obtained in Eq. (62) will be combined with various factors from the scattering amplitude, such as compensating factors of ±i and factors of 1 2 from the contact limit of double-Born interactions. But the important feature is that the color structure C 0 is only corrected multiplicatively by additional rescatterings. Since for the vacuum structure M a = t a and C 0 = C F , we can simply read off the color factors for each type of rescattering from the LO contributions in the opacity expansion, which were already computed in Eqs. (45) and (54). These color factors are therefore preserved as we iterate the scatterings to higher orders in opacity.
Given the above inductive proof regarding the color structure, we will now proceed to promote the 17 LO diagrams Eqs. (45) and (54) to the kernel of a recursion relation which expresses distribution at N th order in opacity in terms of the distribution at (N − 1) st order. We will do this by starting with the last scattering to take place on the jet and evolving the jet history backwards in (light-front) time in both the amplitude and complex-conjugate amplitude, which is possible because the procedure laid out in Sec. III B is manifestly ordered in the positions of the rescatterings. Because the rescatterings can take place differently in the amplitude and complex-conjugate amplitude, the extension of Eqs. (45) and (54) will in general be off-forward, possessing partons of different momenta, or even different partons entirely at a given time in the jet history. Keeping track of the upper limits of the gluon emission times separately in the amplitude and complex-conjugate amplitude is therefore very important in constructing the phases accumulated by the jet.
Let us define a set of four off-forward amplitudes at N th order in opacity, designating the partonic content of the amplitude and complex conjugate by F = "Final" for the quark-gluon system after splitting and I = "Initial" for the quark jet before splitting. The Final/Final , Initial/Final , Final/Initial , and Initial/Initial amplitudes are then written where k (k ′ ) is the relative transverse momentum of the gluon with respect to the quark, x + (y + ) is the upper limit on the gluon emission time in the amplitude (complex-conjugate amplitude), and p is the transverse center-of-mass momentum of the jet (quark before splitting, or quark + gluon after splitting). The final measurement of the gluon distribution within the quark jet is given by the Final/Final amplitude with forward kinematics, along with the two-particle phase space: with the understanding that e ±i∆E − ∞ + → 0 due to the iǫ regulator from the Feynman propagator, as in Eq. (24).
We have also factored out the vacuum color factor C F , so that any additional color factors are the standard ones written in Eqs. (45) and (54).
With this, we can write down the off-forward generalizations of Eqs. (45) and (54) as the kernel of a recursion relation. We emphasize that different diagrams can mix the different Initial/Final sectors with each other. For instance, the direct diagram D 2 takes place entirely in the Final/Final sector, such that preceding scatterings can also take place in the Final/Final sector. On the other hand, the direct diagram D 1 resolves the gluon emission in both the amplitude and complex-conjugate amplitude, such that preceding scatterings take place in the Initial/Initial sector. The recursion relation for the Final/Final sector, shown in Fig. 10 is given by and is the most general expression, receiving contributions from all 17 diagrams. The terms which resolve a splitting in the amplitude (complex-conjugate amplitude) include an explicit factor of the wave function ψ (ψ * ) during the transition from a Final-state to an Initial-state sector. If all of the medium scatterings in either amplitude occur in the final state, without resolving the splitting, then the corresponding wave function is contained in the initial conditions Eqs. (72) instead of the evolution kernel. We also note that the upper limit of the z + integration depends on the minimum of the times x + and y + , and since z + ≤ x + , y + , it is z + which will set the upper limit on the next iteration of this recursion relation.
The recursion relation for the Initial/Final sector is given by a subset of these contributions: and the Final/Initial sector is essentially the complex conjugate: The Initial/Initial sector has already resolved both gluon splittings, and thus depends on x + , y + only through the minimum which sets the upper limit of the z + integral: which is just the recursion relation ("reaction operator") for jet broadening, see e.g. Eq. (8.8) of Ref. [35].
The last ingredient to the recursion relations is the initial conditions, which simply correspond to the vacuum splitting (where applicable), along with the initial distribution of quark jets: It is interesting to note that the recursion relations between these four functions can be cast in the form of a matrix equation, with a particularly simple triangular form due to their causal structure: where the matrix of integral kernels K 1−9 is an explicit representation of the reaction operator. Note that, being complex conjugates, the Initial/Final and Final/Initial sectors are mutually exclusive and do not mix, leading to the vanishing of the {2, 3} and {3, 2} elements of the kernel in Eq. (73). We can write these integral kernels themselves explicitly using shift operators as K 9 = e −q·∇ p e +(z + −x + )∂ x + e +(z + −y + )∂ y + + − 1 2 e +(z + −x + )∂ x + + e +(z + −y + )∂ y + .
The recursion relations Eqs. (68), (69), (70), (71) (or equivalently, the matrix equation Eq. (73) with the reaction operator Eq. (74)), together with the initial conditions Eqs. (72) and the final observable Eq. (67), are one of the main results of this work. We emphasize that these results have expressed the jet substructure observable in terms of universal ingredients, the light-front wave function and the light-front energy denominator, which apply for any partonic splitting and employ exact kinematics. The only kinematic approximations necessary to obtain this result are the eikonal approximation Eq. (8) for scattering in the external potential and the Gaussian averaging Eq. (15). With this, the calculation of any desired order in opacity is cumbersome but straightforward. Although the output rapidly becomes so large that it is impractical to manage by hand, the algebraic computation of Eq. (73) is quick and straightforward on a computer.
We also note that the upper triangular structure of the reaction operator in Eq. (73) itself is particularly amenable to an analytic solution. Because the reaction operator in this form has essentially already been reduced under Gauss-Jordan elimination, the coupled set of four difference (differential) equations can be reduced down to four independent ordinary difference (differential) equations. For example, the recursion relation for the Initial/Initial sector Eq.
with σ(r) ≡ d 2 q (2π) 2 e iq·r dσ el d 2 q and p + dN0 d 2 r dp + ≡ d 2 p (2π) 2 e iq·r p + dN0 d 2 p dp + the Fourier transforms of the elastic scattering cross-section and initial source distributions, respectively. With the use of the explicit solution Eq. (75), the recursion relations for the Initial/Final and Final/Initial sectors Eqs. (69) and (70) both become decoupled equations for their respective functions and can be solved similarly. We leave the implementation of this strategy for future work.

V. RESULTS
Having fully constructed the recursion relations Eq. (73), one immediate application we can pursue is the construction of the second order in opacity gluon-in-quark-jet distribution. A straightforward but tedious application of Eq. (73) yields the general form with the functions N 1−4 containing the various color factors, wave functions, and cosines associated with the diagrams. We note that the expressions to follow represent the sum of the 221 Feynman diagrams which contribute at second order in opacity. In the resulting expressions, we denote δz 1 ≡ z + 1 − x + 0 and δz 2 ≡ z + 2 − z + 1 for brevity, leading to: The exact results Eqs. (76) and (77) at second order in the opacity expansion are the second major result of this work, presented here for the first time. Although there is no existing calculation in the literature against which this exact result can be compared, we can compare with Ref. [13] under the broad source and small-x approximations,

VI. OUTLOOK AND CONCLUSIONS
In the past several years, the QCD community has taken important steps toward the realization of an electron-ion collider in the near future. For this endeavor to be successful, the scientific potential of the new facility must be fully explored. One area that can benefit from new developments is the investigation of cold nuclear effects on hadron and jet production in semi-inclusive deep inelastic scattering. A dedicated program at an EIC can shed new light on the transport properties of large nuclei across across different energy regimes and complement the extensive studies of jet quenching and hadronization in relativistic heavy ion collisions. It can also unify communities and directions of research in high energy nuclear physics that have so far followed separate paths.
To this end, we set out to calculate the quark branching (gluon emission off of an energetic quark) beyond the soft gluon approximation and to an arbitrary order in the correlation between the multiple interactions of the parton system in the medium (opacity). To complete this task required the development of a new theoretical framework which is easily generalizable to all in-medium splitting processes. In this work, we established the new formalism by constructing a new recursion relation Eq. (73) and Eq. (74) in the opacity expansion approach to jet-medium interactions. This equation, which has been derived using only the eikonal approximation to the medium scattering, can be used to calculate the medium modification of jet substructure to any finite order in opacity. The fundamental ingredients, the light-front wave functions Eq. (25) and Eq. (31) and energy denominators Eq. (27), are universal and easily generalized to other partonic splitting processes or to include mass effects. With this recursion relation at hand, it is straightforward (although cumbersome) to calculate the medium-induced parton splitting kernel input to jet substructure observables to any finite order in opacity with exact kinematics.
We have also applied this recursion relation to calculate the exact gluon-in-quark-jet distribution at second order in opacity Eq. (76) and Eq. (77), presented here for the first time. The validation of this result against the known soft-gluon limit, together with multiple validations of the exact first order in opacity result Eq. (55), comprises a fairly robust test of the kernel Eq. (74). These new results are fully ready to be incorporated into existing phenomenology [16,41,[43][44][45][46]78] for calculations ranging from inclusive light and heavy hadron production to jet substructure. These second-order-in-opacity results beyond the soft-gluon approximation may improve these calculations in significant ways, especially for more differential observables such as the production of hadrons inside jets [17,79], by potentially leading to a smoother behavior of the splitting kernels as a function of the kinematic variables and reducing the regions of phase space where the distributions can be negative at any order in opacity [13]. While these methods rely on the ability to truncate the opacity expansion at finite order, the particularly simple triangular structure of the recursion kernel in Eq. (73) makes it possible to decouple the system of equations into 4 independent difference (differential) equations. In future work, we will attempt to solve this system sequentially to obtain a fully resummed expression for the jet substructure distributions. If successful, this approach will be valid for any value of the opacity.
It will also be interesting to study numerically the approach of these distributions to the known analytic limits. For instance, the soft-gluon limit Eq. (80) of the second-order in opacity distribution is tremendously more compact than the full result Eq. (77); a systematic study of the error involved in making such an approximation could be used to increase calculational efficiency. Similarly, by systematically extending the opacity expansion to higher orders, it will be possible to explore the matching onto the continuous path-integral descriptions of the medium interactions. Such a study could yield quantifiable, reliable comparisons of the different descriptions of the medium, along with a better understanding of the boundary between the region of applicability of finite-opacity versus resummed methods.
Finally, coming back to our original motivation, it will be important to incorporate these results into the detailed phenomenology of cold nuclear matter effects expected at a EIC. The use of jet tomography as a probe of QCD matter has been rightly studied in detail for use in heavy-ion collisions, but the jet component of the science program for a future EIC has not received nearly as much attention. In spite of this, the formalism for the how a jet couples to a QCD medium is universal, and the architecture developed for the heavy-ion program at RHIC and the LHC has the potential to add considerable strength to an EIC science program. Because of this, we believe the new results presented here constitute an important advance in the study of jets in both hot and cold nuclear environments.