Measurement of inclusive and leading subjet fragmentation in pp and Pb-Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV

This article presents new measurements of the fragmentation properties of jets in both proton-proton (pp) and heavy-ion collisions with the ALICE experiment at the LHC. We report distributions of the fraction $z_r$ of transverse momentum carried by subjets of radius $r$ within jets of radius $R$. Charged-particle jets are reconstructed at midrapidity using the anti-$k_{\rm{T}}$ algorithm with jet radius $R=0.4$, and subjets are reconstructed by reclustering the jet constituents using the anti-$k_{\rm{T}}$ algorithm with radii $r=0.1$ and $r=0.2$. In pp collisions, we measure both the inclusive and leading subjet distributions. We compare these measurements to perturbative calculations at next-to-leading logarithmic accuracy, which suggest a large impact of threshold resummation and hadronization effects on the $z_r$ distribution. In heavy-ion collisions, we measure the leading subjet distributions, which allow access to a region of harder jet fragmentation than has been probed by previous measurements of jet quenching via hadron fragmentation distributions. The $z_r$ distributions enable extraction of the parton-to-subjet fragmentation function and allow for tests of the universality of jet fragmentation functions in the quark-gluon plasma (QGP). We find no significant modification of $z_r$ distributions in Pb-Pb compared to pp collisions. However, the distributions are also consistent with a hardening trend for $z_r<0.95$, as predicted by several jet quenching models. As $z_r \rightarrow 1$ our results indicate that any such hardening effects cease, exposing qualitatively new possibilities to disentangle competing jet quenching mechanisms. By comparing our results to theoretical calculations based on an independent extraction of the parton-to-jet fragmentation function, we find consistency with the universality of jet fragmentation and no indication of factorization breaking in the QGP.

Properties of the QGP can be inferred by comparing jets in heavy-ion collisions, which traverse the QGP, to their counterparts in proton-proton collisions. Significant experimental and theoretical effort has been made to measure and calculate the modification of jet observables in heavy-ion collisions, known as jet quenching, in an ongoing attempt to achieve a unified description of the jet-QGP interaction [27][28][29]. It has been established that jets traversing the QGP emit soft medium-induced radiation outside of the jet cone, causing their observed yields to be significantly suppressed as compared to pp collisions, and that the fragmentation pattern of the resulting observed jets can be both narrowed and hardened in certain regions of phase space [30][31][32][33][34][35][36][37][38][39][40][41][42][43]. However, the precise role of several theoretical mechanisms of jet quenching, such as quark vs. gluon energy loss, factorization breaking, and color coherence remains unclear. It is essential to understand and disentangle these mechanisms in order to use jet quenching to reveal physical properties of the QGP.
In this article, we consider measurements of charged-particle subjets, defined by first clustering inclusive charged-particle jets with the anti-k T algorithm [44] with radius R, and then reclustering the jet constituents with the anti-k T algorithm with subjet radius r < R, as illustrated in Fig. 1 [45]. We focus on the fraction of charged-particle transverse momentum (p T ) carried by the subjet: where p ch (sub)jet T is the transverse momentum of the charged-particle (sub)jet.
In pp collisions, both the inclusive and leading subjet z r distributions have been calculated perturbatively for a variety of r and R values [46,47]. These calculations involve several interesting aspects that can be studied experimentally such as the role of threshold resummation, and, in the leading subjet case, nonlinear evolution of the jet fragmentation function. Similar effects have recently been examined in e + e − collisions [8,47].
In heavy-ion collisions, subjet observables have been proposed as sensitive probes of jet quenching [46][47][48][49]. The subjet z r observable presents several unique opportunities to study jet quenching: 1. Probe high-z fragmentation. The subjet fragmentation distribution z r is complementary to the longitudinal momentum fraction z of hadrons in jets [36][37][38][39]. The subjet fragmentation distribution can be understood as a generalization of the hadron fragmentation distribution, where in the limit r → 0 the two become equal. By generalizing to r > 0, subjet measurements offer the benefit of probing higher z values than hadron measurements; even for small subjet aperture the z r distribution exhibits a peak at large z r , whereas the hadron z distribution populates much lower z values. This enables a more precise handle on the overall hardness of jet fragmentation, and thereby makes it possible to access a quark-dominated sample of jets [47] at high z r . This in turn provides an opportunity to expose the interplay of soft medium-induced radiation with the relative suppression of gluon vs. quark jets, and introduces a method to do so based on inclusive jet samples alone, complementary to existing methods comparing inclusive and photon-tagged jet observables [37,39]. 2. Test the universality of jet fragmentation in the QGP. In vacuum, it is expected that the parton-to-R r Figure 1: Cartoon of a subjet of radius r inside a jet of radius R. We consider charged-particle subjets clustered using the anti-k T algorithm from the constituents of inclusive charged-particle jets. jet fragmentation function, J(z), is equal to the parton-to-subjet fragmentation function J r (z) [46]. However, it is unknown whether this universality of jet fragmentation functions holds in the QGP independently of the observable considered [50]. Measurements of z r distributions are directly sensitive to the medium-modified parton-to-subjet fragmentation function, J r,med (z), and can be used to extract it. The extracted J r,med (z) can then be compared to an independently extracted mediummodified parton-to-jet fragmentation function, J med (z), to test the universality of in-medium jet fragmentation across different jet observables in heavy-ion collisions and look for signs of factorization breaking. 3. Measure energy loss at the cross section level. A well-defined method of measuring out-of-cone energy loss at the cross section level has been proposed by computing moments of the leading subjet z r distribution [47]. The first moment or "subjet energy loss" describes the fraction of jet p T not carried by the leading subjet, and higher moments describe fluctuations in this energy loss. These quantities can be computed in both pp and Pb-Pb collisions for a variety of r and R values, and contrasted with other measures of jet modification.

Experimental setup and data sets
A description of the ALICE detector and its performance can be found in Refs. [51,52]. The data sample of pp collisions used in this analysis was collected in 2017 during the LHC Run 2 at √ s = 5.02 TeV using a minimum-bias (MB) trigger defined by the coincidence of the signals from the two V0 scintillator arrays in the forward region [53]. The Pb-Pb data set was collected in 2018 at √ s NN = 5.02 TeV. A central collision trigger was used which selects events in the 0-10% centrality interval based on the multiplicity of produced particles in the V0 detector acceptance [54,55]. In the event, the primary vertex was required to be within 10 cm along the beam axis from the center of the detector. Beaminduced background events were removed using timing information from the V0 detectors and, in Pb-Pb collisions, from two neutron Zero Degree Calorimeters located ±112.5 m along the beam axis from the center of the detector. Pileup events were rejected based on multiple reconstructed vertices and tracking selections [30]. After these selections, the pp data sample contains 870 million events and corresponds to an integrated luminosity of 18.0 ± 0.4 nb −1 [56]. The Pb-Pb data sample contains 92 million events in 0-10% central collisions, corresponding to an integrated luminosity of approximately 0.12 nb −1 .
This analysis uses charged-particle tracks reconstructed based on the information from the Time Projection Chamber (TPC) [57] and the Inner Tracking System (ITS) [58]. While track-based observables are collinear-unsafe [59][60][61], they can be measured with greater precision than calorimeter-based observables, and recent measurements have demonstrated that for many substructure observables track-based distributions are compatible with the corresponding collinear-safe distributions [4]. We define two types of tracks: global tracks and complementary tracks [52]. Global tracks are required to include at least one hit in the Silicon Pixel Detector (SPD) comprising the two innermost layers of the ITS and to satisfy multiple criteria on the quality of track reconstruction in the TPC and its pointing to the collision vertex. Complementary tracks are all those satisfying all the selection criteria of global tracks except for the request of a point in the SPD. They are refitted using the primary vertex to constrain their trajectory in order to preserve good momentum resolution, especially at high transverse momentum (p T ). Including this second class of tracks ensures approximately uniform azimuthal acceptance, while providing similar p T resolution to tracks with SPD hits. Tracks with p T > 0.15 GeV/c were accepted over the pseudorapidity range |η| < 0.9 and azimuthal angle 0 < ϕ < 2π. The momentum resolution σ (p T )/p T of the accepted tracks was estimated from the covariance matrix of the track fit parameters [52], and is approximately 1% at track p T = 1 GeV/c and 4% at p T = 50 GeV/c.
The instrumental performance of the detector was estimated with a simulation performed using PYTHIA8 Monash 2013 [62,63] for the event generation and the GEANT3 transport code [64] to propagate particles through the simulated ALICE apparatus. The tracking efficiency in pp collisions is approximately 67% at track p T = 0.15 GeV/c, and rises to approximately 84% at p T = 1 GeV/c, and remains above 75% at higher p T . Studies of the centrality dependence of the tracking efficiency in a HIJING [65] simulation demonstrated that the tracking efficiency is approximately 2% lower in 0-10% central Pb-Pb collisions compared to pp collisions, independent of track p T .

Analysis method
Jets were reconstructed from charged-particle tracks with FastJet 3.2.1 [66] using the anti-k T algorithm with E-scheme recombination and radius (or resolution parameter) R = 0.4 [44,67]. Subjets were reconstructed by reclustering the jet constituents using the anti-k T algorithm with E-scheme recombination and radii r = 0.1 and r = 0.2. The pion mass was assumed for all jet constituents. Jets containing tracks with p T > 100 GeV/c, corresponding to < 1% of the jet sample in the considered kinematic range, were discarded in order to ensure good momentum determination. Jets in heavy-ion collisions have a large uncorrelated background contribution due to the underlying event (UE) [68]. The event-by-event constituent subtraction method was used before jet finding was performed [69,70]. This corrects the overall jet p T and its substructure simultaneously by subtracting UE energy constituent by constituent. A maximum recombination distance R max = 0.25 was used [70]. In pp collisions, the UE was not subtracted and must be included as a model component in all theory comparisons. According to PYTHIA8, the impact of the UE in pp collisions on the leading z r distribution is <3% for z r < 0.95, and grows up to 13% as z r → 1. The jet axis is required to be within the fiducial volume of the TPC, η jet < 0.9 − R, where η jet is the jet axis pseudorapidity.
The jet reconstruction performance is studied by simulating pp events and particle transport through the ALICE detector material as described in Sec. 2. We compare PYTHIA8 generated jets at "truth level" (before the particles undergo interactions with the detector) to those at "detector level" (after detector simulation). The truth-level jet was constructed from the charged primary particles of the PYTHIA8 event, defined as all particles with a mean proper lifetime larger than 1 cm/c, and excluding the decay products of these particles [71]. For the Pb-Pb data analysis, we embedded the simulated pp events after track reconstruction into 0-10% centrality Pb-Pb measured events to account for background effects, and applied the constituent subtraction procedure described above. A jet matching procedure is used in order to associate jets at the truth level to jets at the detector level. In pp collisions, this matching procedure is based on geometrically matching jets within ∆R < 0.6 R, where ∆R = ∆y 2 + ∆ϕ 2 is their rapidity-azimuth separation. In Pb-Pb collisions, the matching procedure includes both this geometrical requirement and further requires that the jet contains at least 50% of the total track p T of the associated reconstructed jet from the embedded pp event at detector level. To study the subjet z r reconstruction performance, a similar matching procedure was adopted. In order to match inclusive subjets at the truth level with their counterparts at the detector level, we apply the same matching procedure as described above for jets, except with the subjet radius r replacing the jet radius R in the geometrical matching criteria. In the case of leading subjets, however, we match the leading subjet at the truth level to the leading subjet at the detector level, without requiring geometrical criteria. The jet energy scale shift  [72]. The mean jet energy scale shift for R = 0.4 charged-particle jets in pp collisions is approximately −13% at p ch jet T,truth = 20 GeV/c and decreases to −21% at p ch jet T,truth = 100 GeV/c. The jet energy resolution is approximately constant at 21% in pp collisions; in Pb-Pb collisions an additional contribution from the underlying event fluctuations further broadens the distribution of the jet energy scale shift, with the standard deviation of the background fluctuation contribution to the shift σ δ p T [68] being equal to approximately 11 GeV/c, independent of p ch jet T,truth . The z r reconstruction performance behaves qualitatively similarly as the jet p T described above, with a peak near zero in the relative residual distribution corresponding to z r,det = z r,truth . The z r reconstruction resolution is O(10%), with asymmetric tails that are broader for small z r than for large z r , and broader in Pb-Pb compared to pp collisions.
Local fluctuations in the underlying event of a heavy-ion collision can result in an incorrect subjet (unrelated to the hard scattering) being identified by the reclustering algorithm. This "mistagging" effect is in exact analogy to the case of identifying groomed jet splittings in the presence of a large underlying event [73], although with improved robustness to mistagging effects [74]. In order to address this issue, the measurement was performed by restricting to z r,det > 0.5, which mitigates these effects. The subjet purity due to these background effects was evaluated by embedding jets simulated with the PYTHIA8 event generator [62] into measured Pb-Pb collisions and following the procedure in Ref. [73] . The residual background contribution remains below 5% for z r > 0.6 and increases up to approximately 20% at z r = 0.5, for the p ch jet T,det range considered in this measurement. This level of background contamination is small enough to allow the results to be unfolded for detector effects and background fluctuations. The corresponding contamination in pp collisions is negligible.
The reconstructed p ch jet T and z r differ from their true values due to tracking inefficiency, particle-material interactions, and track p T resolution. Moreover, in Pb-Pb collisions, background fluctuations significantly smear the reconstructed distributions of z r . To account for these effects, we simulated events at the truth level and detector level and matched both jets and subjets at the truth level and detector level as described above. We constructed a 4D response matrix (RM) that describes the detector and background response in p ch jet T and z r : R p ch jet T,det , p ch jet T,truth , z r,det , z r,truth . For both the pp and Pb-Pb results, we use PYTHIA8 to generate the RM since previous studies indicate that the effect of changing fragmentation patterns on the response is small. This is due to the fact that there is only an indirect dependence via the response matrix [34,38,75]; an uncertainty due to this model dependence is assessed in Sec. 4. A simultaneous unfolding was then performed in p ch jet T and z r using the iterative Bayesian unfolding algorithm [76,77] implemented in the RooUnfold package [78]. The prior distribution for the unfolding procedure is taken to be the truth-level distribution from PYTHIA8. In Pb-Pb collisions, lower limits of p ch jet T,det > 60 GeV/c and z r,det > 0.5 are imposed on the data that is input to the unfolding, in order to reject combinatorial jets and mistagged subjets. No such limitation is imposed on p ch jet T,truth or z r,truth during the unfolding process. The distributions were corrected for "misses", in which a jet was generated inside the considered truth-level range but not inside the detector-level range. In pp collisions the rate of misses is < 3%, whereas in Pb-Pb collisions the rate of misses ranges from 24 − 40% due to the aforementioned fluctuations in the UE. The rate of "fakes", in which a jet exists inside the considered detector-level range but not inside the truth-level range, is < 1% and therefore negligible. The number of iterations, which sets the strength of regularization, was chosen by minimizing the quadratic sum of the statistical and systematic unfolding uncertainties described in Sec. 4. This results in the optimal number of iterations equal to 3 in all cases.
To validate the performance of the unfolding procedure, we performed refolding tests, in which the RM is multiplied by the unfolded solution and compared to the original detector-level spectrum. We also did closure tests, in which the shape of the input MC spectrum is modified to account for the possibility that the true distribution may be different from the MC input spectrum (using the same scalings as for the systematic variations in the unfolding prior described below in Sec. 4). In all cases, successful closure was obtained within statistical uncertainties. Additionally, we performed a closure test to quantify the sensitivity of the final result to combinatorial jets and background subjets. This consisted of redoing the entire analysis on "combined" events containing a PYTHIA8 event and a thermal background, in which "combined" jets were clustered from the combination of PYTHIA8 detector-level particles and thermal background particles. The background was modeled by generating N particles with p T taken from a Gamma distribution, f Γ (p T ; β ) ∝ p T e −p T /β , where N and β were fixed to roughly fit the R = 0.4 δ p T distribution in Pb-Pb data [68]. This background model was verified to describe the subjet purity to percent-level accuracy. The test consisted of constructing the combined detector-level jet spectrum, building the RM, unfolding the combined jets, and comparing the spectrum to the truth-level PYTHIA8 spectrum. Because the background does not contain hard jets, this test is able to quantify the extent to which the analysis procedure recovers the signal distribution and is free of background contamination. The unfolded combined jet spectrum was found to be consistent with the truth-level spectrum within statistical uncertainties, thus confirming a successful closure.

Systematic uncertainties
The systematic uncertainties in this measurement are due to the tracking efficiency, the unfolding procedure, the model dependence of the event generator, and in the case of Pb-Pb collisions, the background subtraction procedure. Table 1 summarizes the systematic uncertainty contributions for pp and Pb-Pb collisions for the two considered values of r. The total systematic uncertainty is calculated as the sum in quadrature of all of the individual sources described below.
The systematic uncertainty due to the uncertainty in tracking efficiency is evaluated by randomly discarding charged tracks before jet finding. The tracking efficiency uncertainty, estimated from the variation of the track selection criteria and a detailed study of the ITS-TPC track-matching efficiency uncertainty, is 4%. In order to assign a systematic uncertainty to the final result, we constructed an alternative RM with this random track rejection and repeated the unfolding procedure. The result was compared to the default result, with the differences in each bin taken as the systematic uncertainty. The uncertainty on the track momentum resolution is a sub-leading effect to the tracking efficiency and is negligible.
Several variations of the unfolding procedure are performed in order to estimate the systematic uncertainty arising from the unfolding correction: -The number of iterations in the unfolding was varied by ±2 units and the average difference with respect to the nominal result is taken as the systematic uncertainty.
-The prior distribution was simultaneously scaled by (p ch jet T ) ±0.5 and a linear scaling in z r by ±50% over its reported range. The average difference between the result unfolded with this prior and the original is taken as the systematic uncertainty.
-The detector-level binnings in z r were varied to be finer and coarser than the nominal binning.
-The lower bound in the detector level charged-particle jet transverse momentum p ch jet T,det range was varied by ±5 GeV/c.
The total unfolding systematic uncertainty is then the standard deviation of the results from the variations, and σ i is the systematic uncertainty due to a single variation, since they each comprise independent measurements of the same underlying systematic uncertainty in the unfolding correction.
The constituent subtraction introduces a bias in the observed distributions, since it implicitly makes a choice of how much p T to subtract from soft particles compared to the hard particles, and similarly for their angular distributions. To estimate the size of the systematic uncertainty related to the background subtraction, we varied R max from "under-subtraction" (R max = 0.05) to "over-subtraction" (R max = 0.7), around the nominal value of R max = 0.25. The maximum deviation of these two variations was assigned as the systematic uncertainty.
The systematic uncertainty due to the model-dependence of the generator used to construct the RM is estimated by comparing results obtained with PYTHIA8 [62] to those obtained with HERWIG7 [79] (in the pp case) or JEWEL 2.2.0 [80,81] (in the Pb-Pb case). For HERWIG7, the default tune was used, and for JEWEL, we adopted the settings described in Ref. [82], with an initial temperature T i = 590 MeV and no recoils. These RMs are then used to unfold the measured data, and the differences between PYTHIA8 and HERWIG7 (in the pp case) or PYTHIA8 and JEWEL (in the Pb-Pb case) are taken as a symmetric uncertainty. We note that in the Pb-Pb case, the results continue to exhibit little-to-no-modification when unfolded with the JEWEL-based response matrix despite that JEWEL itself exhibits large modifications as shown in Figs. 5,6 confirming that the experimental results remain under good control despite the large variations in the jet quenching models themselves.

Results
We report the z r distributions for r = 0.1 and r = 0.2 in both pp and Pb-Pb collisions. All presented results use R = 0.4 jets reconstructed from charged particles at midrapidity, and are corrected for detector effects and (in Pb-Pb collisions) underlying-event fluctuations. We report results for p ch jet T between 80 and 120 GeV/c in both pp and Pb-Pb collisions, as well as a result with finer binning in z r for p ch jet T between 100 and 150 GeV/c in Pb-Pb collisions. The distributions are reported as normalized differential cross sections, 1 where N jet (σ jet ) is the number (cross section) of inclusive charged-particle jets within the given p ch jet T interval, and N (σ ) is the number (cross section) of subjets. With this normalization, the integral of Eq.
(1) is equal to the average number of subjets per jet: Using the leading subjet distributions, we also compute the average "subjet energy loss": which describes the fraction of p T inside the jet that is not contained within the leading subjet [47].  Figures 2 and 3 show the measured z r distributions in pp collisions at √ s = 5.02 TeV for inclusive and leading subjets, respectively. For z r > 0.5 the leading and inclusive subjet distributions are identical, as expected. In this region, the amplitude of the z r distribution increases with z r , with a more pronounced peak at large z r for r = 0.2 than for r = 0.1 since larger subjets are more likely to capture a larger fraction of the jet energy. It is expected that as z r → 1, the distributions will eventually decrease due to the increased splitting probability of soft emissions [47]. This is, however, not visible in the data due to the coarseness of the bin sizes. As z r becomes small, the inclusive subjet distribution grows due to soft radiations emitted from the leading subjet, whereas the leading subjet distribution falls to zero. The fraction of p T inside the jet that is not contained within the leading subjet is z loss = 0.21 for r = 0.1 and decreases to z loss = 0.10 for r = 0.2.  The distributions are generally well described by PYTHIA8 [62,63], however some tension is observed in the largest z r bin. This may be due to threshold logarithms of 1 − z, which may contribute significantly at all orders in the strong coupling constant α s and are not directly included in PYTHIA8 [47]. In addition, hadronization effects are expected to be significant at large z r , since hadronization causes a smearing of the fragmentation across the boundary of the subjet and away from z r = 1. However, due to ill-defined perturbative accuracy in general-purpose MC generators such as PYTHIA and the fact that they are tuned to reproduce data, it is difficult to draw detailed physics conclusions from their comparison to data. In order to study these effects in greater detail, we turn our attention to comparisons with analytical calculations based on perturbative QCD (pQCD).

Subjet fragmentation in proton-proton collisions
Theoretical calculations of the z r distribution in pp collisions have been carried out within the Soft-Collinear Effective Theory (SCET) framework [83] for both inclusive and leading subjets for a variety of r and R values [46,47]. These calculations include all-order resummations of large logarithms of the jet radius and threshold logarithms to next-to-leading logarithmic (NLL ) accuracy. In order to compare these predictions to our measurement using charged-particle jets, a "forward folding" procedure based on MC event generators is applied to account for the fact that we measure only the charged component of jets [6]. Although the calculations are provided with hadronization corrections included [47], we additionally applied a bin-by-bin correction to account for multi-parton interactions using the procedure outlined in Ref. [6]. Figure 4 compares the measured z r distributions to NLL calculations for inclusive and leading subjets [46,47]. Two sets of NLL results are reported in the figure, which are obtained using either PYTHIA8 [62] or HERWIG7 [79] to account for charged-particle corrections, which show generally similar behavior. Uncertainties on the analytical predictions were estimated in Refs. [46,47] by varying the combinations of scales that emerge in the calculation. The softest of these scales determines a transition between the perturbative and non-perturbative regimes: where Λ is the energy scale at which α s becomes non-perturbative. To denote this transition, we draw a dashed vertical blue line at Λ = 1 GeV/c, taking p T to be the weighted average p ch jet T in the considered interval scaled by 120% to approximately translate the p T scale from charged-particle jets to full jets. While we draw the line at a discrete value in order to provide guidance, we remind the reader that the transition from values of z r that are dominated by perturbative versus non-perturbative physics is actually smooth. Note that we display the non-perturbative transition only at large z r , although a similar transition occurs at small z r which is not addressed in this study [84,85]. In the results shown in Fig 4, the cross section is scaled according to the integral of the distribution in a subset of the perturbative region,  TeV, compared to NLL predictions carried out with SCET [46,47] and corrected for missing neutralparticle energy and multi-parton interaction effects using PYTHIA8 [62] or HERWIG7 [79]. The shaded bands denote systematic uncertainty on the NLL calculations. The distributions are normalized such that the integral of the region defined by 0.7 < z r < z NP r is unity, where z NP r is denoted by the dashed vertical blue lines. The nonperturbative scale in Eq. 4 is taken to be Λ = 1 GeV/c. In determining the normalization, bins that overlap with the dashed blue line are considered to be in the non-perturbative (right) region.
In the inclusive subjet case, the measured z r distributions are generally in agreement with the SCET calculations within uncertainties in the intermediate region 0.1 z r 0.9. The calculations begin to diverge from the data at large z r in the non-perturbative regime, where the theoretical calculations are expected to break down. For r = 0.1, the calculations can, in fact, describe the data at large z r within the systematic uncertainties of the calculation, whereas for r = 0.2 this is no longer the case. At small z r , the calculations diverge from the data, with a large overestimate of the peak at z r < 0.02 relative to the magnitude of the theoretical uncertainties. The calculations do not include a resummation of large logarithms of small z r , which suggests that such a resummation is needed to describe the data. The measured distributions can serve to test future calculations that include small z r resummation, which is relevant for attempts to calculate hadron observables using perturbatively calculable jet functions in the r → 0 limit [84,85].
In the leading subjet case, we observe identical behavior to the inclusive subjet case in the region z r > 0.5, since the distributions coincide. For z r < 0.5, where the inclusive and leading distributions differ, the calculations involve nonlinear evolution of the jet fragmentation function. In this region, the NLL results are consistent with the data within the uncertainties of the calculation. Note that this comparison involves only a single bin due to the fact that the z r distributions fall steeply as z r decreases. More differential data would allow for stricter tests of the non-linear evolution of leading jet fragmentation functions. A test of related calculations of leading dijet energy spectra have recently been examined in e + e − collisions [8], finding good agreement. Finally, we note that the quantity z loss described in Eq. 3 is strongly affected by the non-perturbative region z r > z NP r , since a significant fraction of the integral is located in this interval. This presents a challenge to the prospects for theoretically calculating z loss , and calls for further theoretical studies. Figures 5 and 6 show the leading subjet z r distributions in pp and 0-10% central Pb-Pb collisions for r = 0.1 and r = 0.2, respectively, with their ratios displayed in the bottom panels. We report a restricted range in z r in Pb-Pb collisions due to the contamination from the underlying event at low-z r , as explained in Sec. 3, although we note that this excludes only a small portion of the leading subjet distribution. The reported distributions are accordingly normalized to the cross section of inclusive charged-particle jets conditioned with z r in the reported range. The relative uncertainties are assumed to be uncorrelated between pp and Pb-Pb collisions, and are added in quadrature in the ratio.

Subjet fragmentation in Pb-Pb collisions
For both r = 0.1 and r = 0.2, the distributions are consistent with no modification of the z r distribution in central Pb-Pb compared to pp collisions. However, the distributions are also consistent with a hardening effect in Pb-Pb compared to pp collisions that reverses as z r → 1. For r = 0.1, we observe consistency with stronger hardening effects than for r = 0.2. To understand the behavior of the z r distribution, we note that in vacuum there are significant differences in the parton-to-subjet fragmentation functions between quarks and gluons, with the fraction of quark-initiated jets increasing with z r [47]. If the QGP suppresses gluon jets more than quark jets, then a hardening of the z r distribution is expected -in line with previous measurements of hadron fragmentation [86]. On the other hand, medium-induced radiations will in general shift the distribution to smaller z r . The competition between these two distinct sources of jet substructure modification -quark vs. gluon suppression and medium-induced radiation -can result in non-trivial modification to the shape of the z r distribution in different intervals of z r . As z r → 1, the jet sample in vacuum becomes almost entirely dominated by quark jets -thereby rendering the quark vs. gluon fraction modification negligible. This presents an opportunity to expose a region of quark-initiated jets depleted by soft medium-induced emissions. Our measurements are qualitatively consistent with such a modification pattern: a hardening of the z r distribution due to the relative suppression of gluon vs. quark initiated jets, followed by a turnover of the distribution as z r → 1 due to medium-induced soft radiations.
We compare the ratio of the measurements in pp and Pb-Pb collisions with several theoretical models of    [46,50,80,81,[87][88][89]. jet quenching: -Medium jet functions [46,50] are a SCET-based calculation obtained by modifying the z r distributions from pp collisions according to the medium-modified parton-to-jet fragmentation functions extracted in Ref. [50]. The quark/gluon fractions in the extracted medium-modified jet function exhibit a relative suppression factor of approximately four between gluon jets and quark jets.
-JETSCAPE [87,90] consists of a medium-modified parton shower calculated with the MATTER model [89] controlling the high-virtuality phase and the Linear Boltzmann Transport (LBT) model describing the low-virtuality phase [88]. The version of JETSCAPE used for this calculation employs a jet transport coefficient,q, that includes dependence on parton virtuality, in addition to dependence on the local temperature and running of the parton-medium coupling. The calculation includes medium recoil particles, and a subtraction of the thermal component of the recoils is performed by summing the transverse momentum of the thermal particles within the jet (subjet) radius and subtracting this from the corresponding jet (subjet) transverse momentum.
-JEWEL [80,81] implements BDMPS-based medium-induced gluon radiation in a medium modeled with a Bjorken expansion. We use JEWEL 2.2.0 with an initial temperature T i = 590 MeV and initial quenching time τ i = 0.4, which provides an accurate description of a variety of jet quenching observables [82]. The impact of medium recoil is studied by displaying results both with and without recoils enabled. In the case with recoils included, the thermal component of the recoils is subtracted with the same method used in the JETSCAPE calculation (which is similar to the "4MomSub" method [82]) except randomly discarding 33% of the thermal particles (which JEWEL assigns to be neutral) in order to account for the fact that our measurement uses charged-particle jets.
For the JETSCAPE and JEWEL simulations, the width of the curves denotes statistical uncertainty. For the "Medium jet functions" calculation, systematic uncertainties are included but are smaller than the width of the plotted curve.
The comparison of our result to the "Medium jet functions" calculation provides a test of universality of jet fragmentation functions in the QGP, since the calculation uses a parton-to-jet function extracted from inclusive jet measurements, and employs it as the parton-to-subjet function in the z r calculation. We find a consistent description of the z r distribution, and therefore consistency with the universality of jet fragmentation in the QGP. While this does not exclude process-dependent effects or factorization breaking, it does place constraints on the magnitude of such effects, and establishes a new avenue to search for them. These measurements can be used to directly extract the parton-to-subjet function in future work and serve as input for global tests of factorization breaking in the QGP (see Ref. [50]).
The JETSCAPE model describes the data well within the precision of our measurement. The JEWEL model, on the other hand, describes the data well for r = 0.1 when recoils are included, but fails to describe the data for r = 0.2 or when recoils are not included. For both r = 0.1 and r = 0.2, there are large differences in the JEWEL predictions depending whether recoils are enabled, suggesting that this observable may be significantly impacted by medium response. In general, it is expected that medium response will soften the z r distribution since it tends to broaden reconstructed jets. We indeed observe this in the results of the JEWEL calculations, where the z r distribution in the largest z r bin is significantly suppressed when recoils are included compared to when recoils are disabled. However, it appears that for r = 0.2, this suppression is significantly stronger than the experimental data allows (noting that the large enhancement observed at smaller z r is necessitated by the suppression at large z r due to the self-normalization condition). This corroborates previous observations that the medium response implementation in JEWEL, which does not include rescattering of the medium response particles, overestimates the impact of medium recoil, but that the calculations with and without recoil generally bracket the experimental data (see e.g. Ref. [32,91]).  In order to test the r-dependence of the z r distribution with higher precision, in Fig. 7 we compare the z r distributions measured in 0-10% central Pb-Pb collisions with r = 0.1 and r = 0.2 for the p ch jet T interval between 100 and 150 GeV/c. The ratio of the two distributions is shown in the bottom panel of Fig. 7. Since the two measurements are correlated, the systematic uncertainty of the ratio partially cancels out. While separate values of σ z r >0.7 are used to normalize the r = 0.1 and r = 0.2 distributions, the normalization factors only differ by the integral of the z r < 0.7 tails and are therefore within a few percent. We refrain from constructing the corresponding pp ratio due to sizeable statistical uncertainties of the recorded pp data set.
We compare the r=0.1 r=0.2 ratio to the JETSCAPE and JEWEL models discussed above. We find that JEWEL fails to describe the ratio either with recoils on or recoils off, but that the two implementations bracket most of the data. However, the JEWEL calculations predict that at large z r , the r = 0.1 "core" of the jet contains more p T relative to r = 0.2 as compared to the experimental data -i.e. that the large-z r fragmentation is narrower in JEWEL than in experimental data.
The JETSCAPE model describes the ratio significantly better, however given the precision of the measurement, we find significant tension in the shape of the distribution with the JETSCAPE prediction. While this may not be immediately obvious by eye, we note that due to the self-normalization condition in the z r distributions in the top panel of Fig. 7, a shift in one point necessitates an opposite shift in the remaining points to ensure the distribution integrates to unity. For example, if the value of the ratio in the rightmost bin 0.98 < z r < 1 of the JETSCAPE calculation were to move down, the leftmost bins near 0.7 < z r < 0.8 would have to compensate by moving up, rendering the calculation unable to reproduce the experimental data. The distribution from JETSCAPE simulations for pp collisions, on the other hand, describes the r=0.1 r=0.2 ratio well; however, one must use caution in interpreting this, since in Sec. 5.1 we have discussed challenges in achieving an accurate description of the pp baseline at large z r . The JETSCAPE calculation exhibits a hint that the r = 0.1 "core" of the jet contains more p T relative to r = 0.2 as compared to the experimental data, similar to JEWEL although with less significant tension.

Conclusions
We have presented new measurements of subjet fragmentation with ALICE. In pp collisions, we find agreement of pQCD calculations with the data in the perturbative regime at intermediate z r and discrepancies at large z r which may imply that threshold resummation and hadronization play important roles as the distribution becomes increasingly non-perturbative. The PYTHIA8 event generator generally describes the data well, however some tension is observed at large z r , which is consistent with these findings given that threshold resummation is not directly included in PYTHIA8. In the inclusive subjet case, we find a disagreement of the pQCD calculations with the data at small z r , suggesting a need to include a small-z r resummation in order to describe the data. These measurements provide future opportunities to study threshold and small-z r resummations, and motivate new measurements extended to even smaller values of r, which are relevant for understanding parton-hadron duality and the interplay of jet observables and hadron observables as r → 0 [84,85].
In heavy-ion collisions, these measurements serve as a key ingredient to study the high-z region of jet quenching and test the universality of jet fragmentation in the QGP. By comparing the z r distributions for r = 0.1 and r = 0.2 to Monte Carlo jet quenching models, we find indications that quenched jets at large z r are narrower in JEWEL and JETSCAPE than in experimental data. By probing large z r , these measurements can isolate a region of quark-dominated jets with an inclusive jet sample alone, offering the potential to expose a sample of jets depleted by medium-induced soft radiation. Together, our measurements demonstrate that while the large-z r region is theoretically challenging to describe in pp collisions due to threshold resummation and hadronization effects, it is a particularly interesting region to study jet modification in heavy-ion collisions. This calls for theoretical investigation of the large-z r region in greater detail. Future measurements of z r in coincidence with other substructure observables such as the groomed jet radius [34] offer the potential to disentangle medium-induced soft radiation effects from differences in the suppression of gluon vs. quark jets. By comparing our measurements to perturbative calculations based on QCD factorization, we find consistency with the universality of jet fragmentation and no indication of factorization breaking in the QGP. These measurements can be used as input to extract the parton-to-subjet fragmentation function in future work and perform global tests of factorization breaking in the QGP.