Low Energy Protons as Probes of Hadronization Dynamics

Energetic quarks liberated from hadrons in nuclear deep-inelastic scattering propagate through the nuclear medium, interacting with it via several processes. These include quark energy loss and nuclear interactions of forming hadrons. One manifestation of these interactions is the enhanced emission of low-energy charged particles, referred to as grey tracks. We use the theoretical components of the BeAGLE event generator to interpret grey track signatures of parton transport and hadron formation by comparing its predictions to E665 data. We extend the base version of BeAGLE by adding four different options for describing parton energy loss. The E665 data we used consists of multiplicity ratios for fixed-target scattering of 490 GeV muons on Xe normalized to deuterium as a function of the number of grey tracks. We compare multiplicity ratios for E665 grey tracks to the predictions of BeAGLE, varying the options and parameters to determine which physics phenomena can be identified by these data. We find that grey tracks are unaffected by modifications of the forward production. Thus their production must be dominated by interactions with hadrons in the backward region. This offers the advantage that selecting certain particles in the forward region is unlikely to bias a centrality selection. We see a strong correlation between the number of grey tracks and the in-medium path length. Our energy loss model does not reproduce the suppression observed in the projectile region. We see an underprediction of the proton production rate in backward kinematics, suggesting that a stronger source of interaction with the nuclear medium is needed for accurate modeling. These results lay an important foundation for future spectator tagging studies at both Jefferson Lab and at the Electron-Ion Collider, where neutron and proton grey track studies will be feasible down to very small momenta.

Energetic quarks liberated from hadrons in nuclear deep-inelastic scattering propagate through the nuclear medium, interacting with it via several processes. These include quark energy loss and nuclear interactions of forming hadrons. One manifestation of these interactions is the enhanced emission of low-energy charged particles, referred to as grey tracks. We use the theoretical components of the BeAGLE event generator to interpret grey track signatures of parton transport and hadron formation by comparing its predictions to E665 data. We extend the base version of BeA-GLE by adding four different options for describing parton energy loss. The E665 data we used consists of multiplicity ratios for fixed-target scattering of 490 GeV muons on xenon normalized to deuterium as a function of the number of grey tracks. We compare multiplicity ratios for E665 grey tracks to the predictions of BeAGLE, varying the options and parameters to determine which physics phenomena can be identified by these data. We find that grey tracks are unaffected by modifications of the forward production. Thus their production must be dominated by interactions with hadrons in the backward region. This offers the advantage that selecting certain particles in the forward region is unlikely to bias a centrality selection. We see a strong correlation between the number of grey tracks and the in-medium path length. Our energy loss model does not reproduce the suppression observed in the projectile region. We see an underprediction of the proton production rate in backward kinematics, suggesting that a stronger source of interaction with the nuclear medium is needed for accurate modeling. These results lay an important foundation for future spectator tagging studies at both Jefferson Lab and at the Electron-Ion Collider, where neutron and proton grey track studies will be feasible down to very small momenta.

I. INTRODUCTION
The process by which energetic quarks and gluons propagate through space and evolve into hadrons, referred to as hadronization or fragmentation, has not yet been understood in terms of the Lagrangian of Quantum Chromodynamics (QCD). The first ingredients for quantitative study of the process are the parton distribution functions (PDFs). Although PDFs for the proton have been studied for decades, new insights into their estimation continue to be unearthed [1]. For scattering from small objects such as the proton, the formalism of QCD factorization [2] in such reactions as semi-inclusive deep inelastic scattering (SIDIS) allows the observed final state to be parameterized as the convolution of PDFs with the conventional fragmentation functions [3] within many kinematic conditions [4], as successfully demonstrated by a recent global extraction of fragmentation functions from many different experiments [5]. However, the information provided on hadronization in that scenario is limited to quantities such as hadron multiplicities and particle production ratios, such as K/π ratios.
Access to completely new information on hadronization can be obtained by implanting the process inside atomic nuclei [6]. While the long distance nature of the interactions within the nucleus might be presumed to preclude naive application of factorization-based methods [7], there is nonetheless recent progress in extending the usual approaches for proton targets to describe data for parton distribution functions in nuclei (nCTEQ15 [8], EPPS16 [9], nNNPDF2.0 [10], and TUJU19 [11,12]) even at low energies and high-x Bj [13,14].
In this paper we press forward on a new and promising front to understanding hadronization in cold nuclear matter, namely the study of so-called grey tracks, protons in the range momentum 0.2-0.6 GeV/c, in lepto-nuclear scattering. This historical term comes from tracking technologies such as nuclear emulsion or streamer chamber measurements, where the visual appearance of the measured charged particle track differs from that of the highest energy particles. Grey tracks (n g ) correspond to lower momentum charged particles with an ionization energy loss much greater than that of minimum ionizing particles. In the context of this paper, we are studying them as a proxy for the interactions of the energetic parton and/or forming hadron with the nuclear medium. The factorization theorem states that the hadronization of a quark or gluon with high momentum is independent of the target, projectile and subprocess [2]. By performing detailed modeling of grey track production with the BeAGLE event generator enhanced by the upgrade of the PyQM module, we can frame the problem of getting microscopic information on such processes from lepto-nuclear scattering such as the data from the E665 fixed-target experiment at Fermilab. This work lays the foundation for future studies at fixed-target experiments at Jefferson Lab and other electron-beam facilities, but more importantly, at future lepton collider experiments such as the proposed electron-ion colliders [15][16][17][18][19][20]. In the context of collider experiments, it is feasible to study extremely low momentum particles such as grey tracks because they emerge from the nucleus having a velocity in the lab frame very close to that of the hadron beam. The charged component can be momentum analyzed using the magnetic fields already present for the recirculating beam transport, while the neutral component can be measured using the zero-degree calorimeter technique [21][22][23]. These techniques, sometimes referred to as spectator tagging or geometry tagging [24][25][26], often require special instrumentation to detect very lowenergy particles associated with highly energetic interactions [27,28]. Geometry tagging studies of grey tracks provide unprecedented new information on the microscopic interactions of the hadronization constituents with the cold nuclear medium.
The concept for studying grey tracks in experimental measurements originated in the 1960's [29] and was discussed through the following decades [30] generally in the context of hadron-nucleus collisions, although it was recognized that deep inelastic scattering by lepton beams would be simpler to interpret [31]. An experimental exploration of these ideas took place for proton-nucleus interactions in the E910 experiment at the Alternating Gradient Synchrotron at Brookhaven National Lab beginning in 1996, measuring slow protons and deuterons produced in the collison of 18 GeV/c protons with Be, Cu, and Au targets [32]. In comparing their data with several models, they concluded that there was a strong linear dependence of the number of projectile-nucleus interactions N int pA on the mean number of grey tracks measured, with a target-dependent constant of proportionality. Defining N int pA as the centrality of the collision, they concluded that measurement of the mean number of grey tracks can determine event centrality well for a given nucleus.
The E665 experimental collaboration published more than 20 papers in the 1990's based on measurements with a secondary muon beam of 490 GeV at Fermi National Accelerator Laboratory. One of these included measurements of grey tracks produced from a xenon target, making a comparison to those seen from a deuteron target [33]. These data stimulated renewed theoretical interest in using grey tracks to study the space-time development of hadronization via the nuclear medium. References [34] and [35] explored the dependence of the mean number of grey tracks on four-momentum transfer Q 2 and x Bj within a theoretical model. They found that the number of grey tracks rises with increasing Q 2 , and found fair agreement with the dependence measured in E665.

II. THE E665 DATA
The E665 data being described by the simulations in this paper consisted of muon-xenon scattering compared to muon-deuterium scattering in Deep Inelastic Scattering (DIS) kinematics. The experiment achieved nearly 4π steradians of acceptance for charged particles by using a streamer chamber immersed in a 15 kG dipole magnetic field as a vertex detector. This enabled study of grey tracks in combination with the high-energy forward tracks, which were analyzed with a downstream spectrometer system. Scattered muons were identified with a downstream iron absorber followed by tracking devices. An electromagnetic calorimeter was used for photon detection to identify muon bremsstrahlung.
The experimental definition of grey tracks used by E665 starts with identification of the scattered muon in the downstream muon detector, with track matching to the trajectories measured upstream. Charged particles not identified as muons were treated as hadrons. Subsequently, a rigorous visual analysis of the photographs of hadron tracks in the streamer chamber was performed. As low momentum particles, they deposited much more energy per unit length than the minimum ionizing particles, thus leaving a much higher streamer density. Grey tracks were counted independently in two laboratories, and these analyses agreed to within 11%. In addition to this visual analysis, they were required to have momentum in the interval 0.2 -0.6 GeV/c. The majority of these particles were protons. As depicted in (a), for this particular event the quark q with impact parameter b1 has a lower energy which on average results in a shorter color lifetime d1/c, and the antiquarkq with impact parameter b2 has a higher energy which on average results in a longer color lifetime d2/c. The propagating quarks begin the evolution into jets or hadrons j1 and j2 either inside or outside the medium. In (b) is shown the simpler process that dominates at higher xBj. The energy and momentum of the virtual photon is absorbed by a single valence quark q with impact parameter b which evolves into jet or hadron j. The degree of misalignment between the virtual photon direction and the direction of leading hadron or jet j due to intrinsic k ⊥ and Fermi momentum is exaggerated in the right-hand diagram (b) to facilitate visualization. In (b) it is possible to speak of the "struck quark" γ * + q → q while (a) corresponds to γ * + g → qq, see text. Small arrows in both diagrams represent low-energy hadrons emitted from the nucleus due to various processes, e.g., potential grey tracks. Features not shown in the figure include: the falloff of the nuclear density at the boundary of the nucleus; other quarks and gluons produced in the process; possible fission of the nucleus; and neutron evaporation from the nuclear fission fragments.

III. GEOMETRY TAGGING FOR DIS
The term "Geometry Tagging" in this context refers to the measurement of particles emerging from the nucleus to determine geometric characteristics of the hard interaction and of the subsequent final-state interactions of the participants, including spectator particles such as evaporation neutrons, fission fragments, hadrons produced from gluonic bremsstrahlung from in-medium partons, and low energy knocked-out protons and neutrons, among others. BeAGLE is able to simulate such processes, and it has been employed in evaluations of the effectiveness of determining geometric features from various experimental signatures, such as those just mentioned, for collider experiments [24]. That work follows earlier studies with the DPMJET3 code for forward neutron tagging [25] in which it was found that the evaporation neutrons could give information on the event centrality.
The impact parameter b is a well-known feature of the scattering process. It is particularly relevant for hadronnucleus scattering since the probability of interaction of the incident hadron with the medium is large due to the magnitudes of hadronic cross sections. Thus, the projectile pathlength through the medium is strongly correlated to the impact parameter, and the number of interactions of the projectile within the medium is indicative of the centrality of the collision, as studied by the E910 collaboration discussed earlier. For lepton beams, which have electroweak cross sections, a more relevant quantity is the pathlength in the medium of the energetic quark and any subsequent hadron that contains it. This is because for higher x Bj interactions, leptons and photons easily penetrate into the nucleus so that the hard interaction can take place anywhere within its volume. For the energetic quark, the quark pathlength, labeled d, can be small for multiple reasons. It can be small because its color lifetime (the time interval between the first interaction of the quark with the medium and the neutralization of its color) is short, reflecting that its most probable value is zero as is the case with ordinary lifetimes, or it can be because the interaction is on the periphery of the nucleus. However, it can only be longest for central collisions for events in which the color lifetime has a large value. The color lifetime can be large either because it fluctuates to a large value, or because of kinematic conditions such as a relative energy z h ≡ E h /ν that is near the optimal value of ≈ 0.3, as indicated by the Lund String Model prediction [36]. As will be discussed later, in BeAGLE we find that on average a larger value of d is correlated with a larger value of the number of grey tracks. A maximal value of d is more suggestive of a colored energetic quark with a long color lifetime in a central collision that produces low energy particles through medium-induced energy loss, and less consistent with an inelastic hadronic collision, which would be more likely to produce higher energy hadrons.
In Fig. 1 is shown an illustration of b and d for low x Bj and for high x Bj , adapting the development of references [37,38]. In what follows we primarily consider the color propagation in the Target Rest Frame (TRF), while clarifying the connections to the Infinite Momentum Frame (IMF) for completeness. At lower x Bj < 0.1, which is the case for most of the E665 data, the quark pair production process becomes possible. The process of producing a virtual photon and qq pair in the target reference frame is distributed over a longer distance in space due to the Ioffe time ≈ 1/(2 · x Bj · M p ) or coherence length [39], with the pair formation naturally occurring as one of the DGLAP radiation loop processes that are proportional to log(Q 2 ). For higher x Bj > 0.1, the simpler process of absorption of the virtual photon by a valence quark, shown on the righthand side of the figure, is dominant.
In the TRF the DIS cross section σ DIS is determined by the scattering of the qq Fock state of the virtual photon on the target (and of higher Fock states at higher orders of α s ). At leading twist there are two qualitatively distinct configurations of the qq pair which give rise to σ DIS ∝ 1/Q 2 : 1. Quarks with similar momenta: the two quarks share the photon longitudinal momentum ∼ ν roughly equally, and have relative transverse momentum ∼ Q. Then the transverse size of the pair is ∼ 1/Q, and their scattering cross section in the target is ∝ 1/Q 2 (color transparency). The pair dominantly scatters through gluon exchange to the target. In the standard IMF frame this corresponds to the physical process γ * + g → qq, which occurs via higher orders of the photon's Fock space. The gluon needs to be highly virtual to resolve the small-sized quark pair, hence this process is suppressed by α s (Q 2 ). Nevertheless this process dominates at low x Bj due to the large gluon distribution, corresponding to the left-hand side of Fig. 1.

2.
Quarks with very different momenta: one of the quarks takes nearly all of the photon momentum ν, while the other's momentum remains fixed in the Bjorken limit. The transverse size of such pairs originates in the intrinsic k ⊥ distribution of the parton, ∼ 1 fm, and so they scatter with a large (Q-independent) cross section.
The probability of such a Fock state in the photon is ∝ 1/Q 2 , accounting for σ DIS ∝ 1/Q 2 . The quark with fixed momentum may be considered to be part of the target wave function. Thus this process corresponds, in the IMF, to the lowest order parton model γ * + q → q shown on the right-hand side of Figure 1.
There is a smooth connection between cases 1 and 2. At low x Bj the quark with Q-independent momentum nevertheless has a large longitudinal momentum ∝ 1/x Bj , and starts to interact via gluon exchange. Such quarks correspond to sea quarks in the IMF, generated by gluon splitting at higher orders of α s .
In the following section we will make use of the PYTHIA6 module within BeAGLE to quantify the relative contributions of these processes in the E665 kinematics.

IV. MONTE CARLO SIMULATION: BEAGLE
The Monte Carlo code "Benchmark eA Generator for LEptoproduction" (BeAGLE) [40] version 1.02 was used to simulate the E665 grey tracks in this work. BeAGLE is a general purpose FORTRAN program for simulating electron-nucleus (eA) interactions. We are testing and using the PyQM module which we updated in the current version of BeAGLE to explore the sensitivity of the grey track observables to partonic energy loss in the cold nuclear medium.
BeAGLE is a hybrid model that uses the DPMJet [41], PYTHIA6 [42], PyQM, FLUKA [43,44], and LHAPDF5 [45] codes to describe high-energy leptonuclear scattering. The geometric density distribution is provided primarily by PyQM while the quark distributions within that geometry are provided by nPDF EPS09 [46]. The parton-level interactions and subsequent fragmentation is carried out by PYTHIA6. Hadronic formation and interaction with the nucleus is described by DMPJet, the impact of the scattering on the nucleus is described by FLUKA, including nucleon and light fragment evaporation, de-excitation by photon emission, nuclear fission, and Fermi breakup of the decay fragments. The PyQM module implements the Salgado-Wiedemann quenching weights to describe partonic energy loss [47].
BeAGLE includes a variety of options to control the phenomena included in the simulation, some of which are mentioned here. Nuclear shadowing is described via two different approaches. The hadron formation time is accounted for in the DPMJet intra-nuclear cascade. Fermi motion of the nucleons in the nucleus can be described with several different mechanisms, or turned off completely. The nuclear geometry parameters from PyQM can be overridden, and deformed nuclei can be described. The PyQM actions available include specification of theq transport coefficient to adjust the degree of interaction between energetic partons with the nuclear environment. Some details of the partonic energy loss process in PyQM can also be selected, such as the fraction of the recoil in PyQM transferred to the nucleus, the modeling options for quark/hadron transverse momentum generated by the medium, the multiplicity and energy of the generated gluons, and other options. The main program is DPMJet, which uses PYTHIA6 to handle elementary interactions and fragmentation. PyQM handles this directly after elementary interactions in PYTHIA6, while DPMJet handles nuclear geometry and, after fragmentation by PYTHIA6, DPMJet takes care of the nuclear evaporation by FLUKA. I. Percentage of the processes in which a "point photon" is used and in which interaction mechanisms are used in PYTHIA6 for the BeAGLE parameter settings used in this work. The BeAGLE hadron formation time parameter tau0 was set to 7 fm. For the "point photon" processes, the hard-scattering final state is composed of a single q orq (LO), a qqg triplet (QCD Compton), or a qq pair (PGF). In all cases an appropriate soft nucleon remnant or remnant cluster is produced by PYTHIA6. For "vector meson production" the final state is V + N, V * +N, V + N * , V * +N * , where * means a diffractively broken up state. In the "resolved process," the γ * is treated as a hadron and interacts through a hard QCD process with the nucleon. Thus, a quark, antiquark or gluon from the γ * interacts with a quark, antiquark or gluon from the nucleon, leaving two remnant systems. The "low-pT" processes collect together remaining categories that contribute to low pT interactions. The numbers in this table are calculated for xenon, but they are very similar for deuterium. As described earlier, at high x Bj the process γ * + q → q dominates, and at low x Bj the process γ * + g → qq dominates. Here we quantify the relative probability of those two kinds of processes according to the modeling of PYTHIA6 in the E665 kinematics. These are quantified in three categories, following Reference [33]: x Bj < 0.02 (shadowing region), x Bj > 0.02 (no shadowing region), and the full range accessed 0.002 < x Bj < 0.3.
As is evident in Table 1, the PYTHIA6 assessment of the various processes indicates that the high-x Bj process described in the previous section, and shown in Fig. 1 (right), is the dominant process by far. Even though the E665 data go as low as 0.002 in x Bj , according to PYTHIA6 it is not low enough to have a significant production of the low-x Bj process described in the previous section and shown in Fig. 1 (left). This will presumably be much more prevalent in high energy colliders such as the future high energy Electron-Ion Collider.

V. MEDIUM-INDUCED GLUON RADIATION: PYQM
In nuclear DIS processes, the partons emerging from hard scattering will, for a time τ c , travel as a colored particle [48], and experience multiple rescatterings on the nucleus remnant. The parton energy degradation -or, parton energy loss -caused by medium-stimulated gluon radiation is an important process to understand the suppression of large momentum hadron production observed in semi-inclusive nuclear DIS processes at various facilities [49].
In BeAGLE, the PyQM module is used to simulate medium-induced gluon radiation and quark energy loss. For each quark and gluon generated by PYTHIA6, the amount ∆E of energy radiated by a parton traveling through a nucleus is stochastically generated using "quenching weights" calculated by Salgado and Wiedemann [47], that provide the probability distribution P (∆E; ω c , R) in the energy ∆E . The model considers a static and uniform medium of length L with scattering power quantified by the transport coefficientq, that measures the average parton transverse momentum broadening per unit path length. The radiation process is regulated by the gluon characteristic energy ω c = 1 2q L 2 . In a finite size medium, large angle radiation suppression is controlled by the cutoff parameter R = ω c L. The geometry of the nuclear medium is taken into consideration by using a Wood-Saxon density distribution [50] to randomly generate the struck quark production point, and calculate its average in-medium path length d. The energy loss ∆E is then calculated using quenching weights with L = d. The partonic final state generated by PYTHIA6 is then modified by removing the appropriate amount of energy from each of the partons generated in the hard scattering, and furthermore adding gluons to the final state, or letting these be absorbed by the nucleus remnant as discussed below. Finally, hadronization of this modified final state is handled as usual by PYTHIA6's implementation of the Lund string fragmentation model. A more detailed description of PyQM can be found in [51]. The in-medium interactions of non partonic states generated by PYTHIA6 -especially soft hadrons -are handled by BeAGLE's DPMJET module.
We have implemented four different energy loss scenarios when interfacing PyQM to the BeAGLE simulation. These differ in the way the radiated gluons are added to the final state generated by JETSET and how much of their energy is absorbed by the nuclear medium, or escapes it.
In the first option, called no gluons, there is no compensation for the energy lost by the hard partons, namely, we neither add gluons to the final system to carry away the lost energy, nor do we redistribute this in the rest of the nuclear remnant. Energy conservation is therefore broken, and we use this option for cross-check purposes.
In the second option, called 1-hard gluon, we attribute the energy loss to the radiation of a gluon generated in a single hard re-scattering on the nuclear medium. This gluon is added to the final state with energy ∆E = E initial parton − E final parton , where E initial parton is the initial energy of the parton that will lose energy, and E final parton is the result of its energy loss, isotropically distributed as transverse momentum of magnitude q 2 T =qL, and a longitudinal momentum appropriate for an on-shell massless gluon.
In the third case, called 1-hard + soft gluons, we consider the energy of a radiated hard gluon based on the multiple soft scattering approximation [47], where the gluon energy spectrum per unit path length is: for ω < ω C . After integrating over the path length and ω, along with considering that the multiplicity of gluons emitted with energies larger than ω is N (w ) = ∞ w dw dI(w ) dw we obtain for one gluon: where L is the average path length,q is the transport coefficient and α s is the strong coupling at soft scales. The energy loss is calculated as ∆E = E initial parton − E final parton , just like the option above, which gives the total energy available to us as ∆E, where ω hard can be either less than or equal to ∆E. To convert energy in cases where ω hard is less than ∆E, we create soft gluons with energy: For technical reasons, FLUKA has difficulties to handle very high energy deposited in the nucleus. To circumvent this issue, we add a constraint on the soft gluons that they must not have energy greater than 5 GeV. In the scenario where they can have this energy or more, we produce a qgq triplet with energy: In summary, in this option we add a hard gluon with energy ω hard while the remaining radiated energy is transferred to the remnant nucleus and a qgq triplet can be created in case of energy excess. In the last option, called soft-gluons, we do not add a gluon to the PYTHIA6 list, but we conserve energy by redistributing the energy loss to the remaining nuclei. We model the process assuming that only soft gluons are radiated with a total energy ∆E = E initial parton − E final parton which is redistributed to the nuclear remnant. Similarly to the previous option, we have to limit the total energy sent to the nucleus and will create a triplet to absorb any extra energy beyond 5 GeV. This option is similar to option one, with the only difference that we conserve the energy.
As an additional remark, for the purposes of calculation of energy loss, we assume that the photon energy, proportional to ν, is large enough to boost the struck parton's lifetime outside of the nuclear target. Consequently, τ c d, and induced gluon radiation occurs all along the struck parton's path through the whole nucleus. This should be a fair approximation for forward tracks at the energy of the E665 experiment considered in this work, as well as at the future Electron-Ion Collider. It is not necessarily true for the target projectile region, even in E665 kinematics, since there is strong experimental evidence of intranuclear cascades associated with backward kinematics. BeAGLE uses an estimate for τ c to decide when each hadron is formed, and if hadrons are produced inside the nucleus, a detailed simulation of the resulting intra-nuclear cascade is initiated. At lower energies, such as at Jefferson Lab, hadronization may occur inside the medium even for forward tracks, and one would also need to include prehadron formation and nuclear absorption in the simulation.

VI. RESULTS: COMPARISON TO E665 EXPERIMENTAL DATA
The analysis [33] made by the E665 experiment compares µXe and µD data and pXe scattering from the NA5 collaboration to study nuclear effects in the hadronic system, with the grey tracks in an event being the main tool in the comparisons. They reported that the average multiplicity n g of the grey tracks is significantly lower with µXe scattering than with pXe scattering, consistent with expectations due to the absence of initial state interactions between the incident particle and the nucleus in the muon scattering data.
The E665 measurements were performed with a muon beam with an average energy of 490 GeV. Since particle identification detector information was not available, a partial identification of protons was used based on the physics expectation, considering that in µD scattering with E µ = 490 GeV and W > 8 GeV about 50% of the positive hadrons with x F (m π ) < −0.2 were protons [33]. To define the DIS region, to avoid kinematic regions where the radiative corrections are large, and to exclude the region where the experimental resolution is poor, they applied the following kinematic cuts: Q 2 > 1 GeV 2 , 8 < W < 30 GeV and x Bj > 0.002. After all cuts, the total number of events in BeAGLE is approximately 4.6 million events for deuterium target and 4.5 million events for xenon target, and 6309 events for deuterium target and 2064 events for the xenon target in Adams et al. [33]. Although in the data the grey tracks were assumed to be predominantly protons, there was a large contamination by pions and kaons, estimated at about 40% and (15 ± 9)% for the µD and µXe data, respectively [33].
In the above reference, a comparison was made between the results for multiplicity ratios with and without grey tracks. It was concluded that grey tracks are very efficient for tagging events where cascade interactions occur because the nuclear effects on hadrons are strongly enhanced in the sample of µXe events containing grey tracks. The sample of µXe events without grey tracks appeared very similar to the sample of µD events, reinforcing the idea that diffractive scattering has a very similar pattern on deuterium and xenon.
In BeAGLE, the same kinematics cuts were applied, with the only exception that the grey tracks (n g ) are directly identified as protons with momentum 0.2 < p < 0.6 GeV.
The multiplicity ratio [33] was defined as R = n(n g ) µXe n µD (5) where n(n g ) µXe is the average number of hadrons in the final state including charged pions and kaons, protons, antiprotons in the xenon target, and n µD is the average number of hadrons and leptons in the final state for the deuterium target. Moreover, the same charge cuts were applied to the denominator as to the numerator, and grey tracks were also included in n(n g ) µXe and n µD . E665 concluded that there were significant diffractive scattering contributions to these data. They estimated that a lower bound of 18 ±3% of all events were diffractive for x Bj less than 0.02 on xenon. It can be anticipated that diffractive events would modify the measured multiplicity ratios. It is also commented in section 4.8 of [33] that diffractive scattering events would have a similar topology for xenon as for deuterium. It goes on to say that such scattering events, which involve particle combinations with the quantum numbers of the vacuum, would naturally produce fewer hadrons than the DIS processes, which involve color octet particles in propagation, on average. As a result, the multiplicity ratio, defined in eq. 5, changes its value if diffractive scattering occurs. The numerator and denominator of the multiplicity ratio are both affected by the presence or absence of diffractive scattering. The numerator, an average over the multiplicity for xenon as a function of the number of grey tracks, would be enhanced by diffractive scattering for low n g . For example, it contributes more events to the n g =0 bin than the pure DIS process would, on average. In the denominator, which is the average deuterium multiplicity, the diffractive scattering would also produce cleaner events than DIS, biasing the average to smaller numbers. If this trend were to be corrected, the denominator should be made larger.
Therefore, the expected effect on the multiplicity ratio, relative to pure DIS events, would be that the numerator is too large for low-n g events, and the denominator is too small for all measured events. Both of these effects imply that for pure DIS events, the multiplicity ratio should be smaller, particularly for low-n g events, with the n g =0 bin being the most strongly affected. As can be seen in Fig. 2, a somewhat smaller multiplicity ratio due to these effects could improve the overall agreement of BeAGLE with the E665 data. We do not attempt to make a correction to account for the effects due to diffractive scattering; the trigger used for acquiring these data, which required multiple charged particles, tended to partly suppress coherent diffractive scattering to an extent that is difficult to estimate.
The regions are divided into different values of rapidity calculated in the photon-lepton center of mass frame as where the first row is the target fragmentation region with rapidity values y < 1, the middle row is the central fragmentation region with rapidity values between −0.5 < y < 0.5 and the last row is the projectile fragmentation region with rapidity values y > 2. Each column indicates the charge of the hadrons and leptons in their final state: the first column stands for charged hadrons and leptons, the second for positive hadrons and the third for negative hadrons and leptons. The E665 collaboration observed that the multiplicity ratio R, eq. 5, with n g in different rapidity regions, reveals diverse characteristics of particle production on a nucleus. In Fig. 2, the multiplicity ratio as a function of the number of grey tracks n g is shown. The grey circles represent the E665 measurements. They conclude that the abundant production of hadrons in the target fragmentation region is due to cascade interactions. In the central rapidity region, additional production of hadrons due to multiple projectile collisions is seen in the pXe scattering, while little additional hadron production is seen in the µXe scattering [33]s. At large n g , depletion of hadrons in the projectile fragmentation region is seen, presumably due to energy loss in projectile or cascade interactions. In the µXe scattering, the events with large n g account for only a tiny fraction of all events. Compared to the pXe data, the µXe data showed a stronger depletion of fast hadrons than would be expected from the number of projectile collisions.
In Fig. 2 we can also see the multiplicity ratio simulated as a function of the number of grey tracks n g forq = 0.5 GeV 2 /fm. The data consider up to seven grey tracks, but in BeAGLE we predicted up to 12.
The different colors indicate the different gluon radiation options in the PyQM module. Red represents no gluon radiation option, pink represents a hard gluon whose energy is ∆E = ω hard , green represents the radiation of a hard gluon with at least energy ω hard (eq. 2) plus soft gluons with energy ω soft eq. (3). The last option in blue indicates that only soft gluons are calculated, which means that all the energy is given to the remaining nuclei, taking into account that the energy of the soft gluons never exceeds 5 GeV.
To compare the different options and the effects of induced gluon radiation in a lepton-nucleus collision, we simulated different values ofq, but we only showq = 0.5 GeV 2 /fm withq = 0.0 GeV 2 /fm as reference.This is because even at this relatively large number, there is essentially no large effect due to energy loss. In Fig. 2 we see the multiplicity ratios as a function of the grey tracks for each case and in Fig. 3 we have rescaled the projectile region to see the small effects between the options and subsequently the energy loss.
In the target fragmentation region, BeAGLE underestimates the ratio for positive particles, where we also have a large number of grey tracks in this region, see Fig. 2. In the case of negative particles, BeAGLE predicts with high accuracy up to 5 grey tracks, but does not show much difference from any of the PyQM options. It shows very large ratios, but independent of energy loss, even with a high number of grey tracks, at least in this region.
In the projectile fragmentation region, BeAGLE accurately predicts up to four grey tracks for positive particles, but overestimates a high number of grey tracks without showing a large suppression effect, even for large values ofq. Therefore, our calculations of energy loss are not considerably different. For negative particles, BeAGLE continues to overestimate the ratios, with no sign of a large suppression for a small number of grey tracks.
The central fragmentation region provides a good description for negative particles, but not for positive particles, underestimating the ratios. In any case, this region is a combination of processes that also come from the other two regions.
In Fig. 4 we present the average number of grey tracks, and all grey tracks per event, compared to the in-medium path length of the struck quark. This length is calculated by BeAGLE from the interaction point to the edge of the nucleus. We observe a strong correlation between the grey tracks and the in-medium path length, despite the fact that the distribution of grey tracks is very large, indicating a general trend that more grey tracks are observed as the distance of the struck quark increases.

VII. CONCLUSIONS
We find that the grey track production is dominated by interactions with in-medium hadrons in the backward region, where grey tracks are quite numerous. We have improved the PyQM module in BeAGLE to offer four options implementing partonic energy loss. Using a comparison of BeAGLE simulation to E665 grey track data we find that grey tracks are unaffected by such modifications for the forward production. We see a strong correlation between the number of grey tracks and the in-medium pathlength for lower values of n g , an important quantity needed for precise modeling and interpretation of geometry-tagged data. This offers the advantage that a selection of certain particles in the forward region is unlikely to bias a centrality selection using the backward region grey tracks.
Our energy loss model does not reproduce the suppression observed in the projectile region data, even with rather large values ofq. This raises questions on what could be at the origin of such a large suppression, unexpected at such high energies. We also see an unambiguous underprediction of the rate of positively charged grey track production in backward kinematics, suggesting that a much stronger interaction with the nuclear medium is needed. At the same time, there is very good agreement with the rate of negatively charged grey track production, indicating that the fragmentation process producing the negatively charged particles is well-described.
These results lay an important foundation for future spectator tagging studies both with CLAS12 at Jefferson Lab, and at the Electron-Ion Collider, where both neutron and proton grey track studies will be feasible down to very small momenta.