Time reclustering for jet quenching studies

The physics program of ultra-relativistic heavy-ion collisions at the Large Hadron Collider (LHC) and Relativistic Heavy-Ion Collider (RHIC) has brought a unique insight into the hot and dense QCD matter created in such collisions, the Quark-Gluon Plasma (QGP). Jet quenching, a collection of medium-induced modifications of the jets' internal structure that occur through their development in dense QCD matter, has a unique potential to assess the time structure of the produced medium. In this work, we perform an exploratory study to identify jet reclustering tools that can potentiate future QGP tomographic measurements with jets at current energies. Our results show that by using the inverse of formation time to obtain the jet clustering history, one can identify more accurately the time structure of QCD emissions inside jets, even in the presence of jet quenching.


Introduction
Modifications of jets in heavy-ion collisions due to interactions with the dense medium created in such reactions lead both to a suppression of the jet cross-section and alterations of their inner structure, collectively denoted as jet quenching. Jets are complex objects that build up a characteristic structure through successive emissions while they propagate in and interact with the background medium. In the absence of such a medium, this scale evolution is calculable in perturbation theory. It is a e-mail: liliana@lip.pt b e-mail: andre.cordeiro@tecnico.ulisboa.pt c e-mail: korinna.zapp@thep.lu.se responsible for the fact that jets carry information from a broad range of scales (from the hard scale set by the hard scattering matrix elements down to the hadronic scale). As jets propagate through the evolving medium, they probe the medium also at different times and positions. The spatio-temporal structure of the background thus gets imprinted together with the scale information on the modified jets. Jets thus offer a unique chance of accessing information on the scale dependence as well as the spatial profile and time evolution of the background medium. Decoding this information from the observed jets is, however, a highly non-trivial task. The reasons for this are manifold and range from the lack of analytical control to the complications of large fluctuations in the radiation pattern. In this paper, we present a first step towards accessing time information by using a non-standard clustering algorithm and give a first example how this new technique can be used to gain information about the time evolution of the medium.
Jet algorithms are the attempt to identify the spray of final state hadrons that result from the fragmentation of a high energy parton (quark or gluon) (see [1] for a review). The sequential recombination jet algorithm family is massively used in both proton-proton (pp) and heavy-ion collisions for a multitude of studies. They start by identifying the pair of particles (i, j), that are closest in a distance measure given by: where with p T i(j) , y i(j) , φ i(j) the transverse momentum, rapidity and azimuthal angle, respectively, of the particle i(j). R is the jet radius that defines the maximum (y, φ) reach of the algorithm, and p a continuous parameter that specifies the jet clustering algorithm. We recover the Cambridge/Aachen (C/A) [2,3], anti-k T [4] and k T [5,6] jet algorithms by setting p = 0, −1 and 1 respectively. The C/A and anti-k T algorithms play a crucial role in jet studies. The latter, typically yielding circular jets and most sensitive to high transverse momentum particles, is usually employed in the identification of signal jets of interest. The former is a purely geometric jet algorithm, that starts clustering particles that are close in (y, φ), and progresses to larger distances. As QCD vacuum emission follows, at leading logarithmic accuracy, angular ordering, i.e., subsequent emissions are restrained from being emitted at larger angles than the previous one, the clustering sequence obtained with the C/A algorithm resembles the QCD radiation pattern.
Since the pioneering work of [7], which suggested to use the information about the internal structure of jets obtained from the jet clustering sequence to identify boosted Higgs bosons decaying hadronically, the study of jet substructure in a variety of contexts has become a very active topic. The general strategy is to identify jets using the anti-k T algorithm and, in a second step, recluster the particles forming jets with a different algorithm, typically C/A . This procedure is often combined with grooming techniques designed to reduce contamination from uncorrelated underlying event activity mostly manifest in soft large-angle structures inside jets. One such procedure is SoftDrop [8], which reclusters jets with the C/A algorithm and subsequently goes backwards through the clustering sequence discarding sub-jets that are soft and/or at large angle. In pp collisions the energy sharing in the first soft-drop unclustering step was shown to be a proxy to the QCD Altarelli-Parisi splitting functions [9].
These techniques, developed in pp collisions, were heavily imported to the heavy-ion field, opening a new window of exploration to deepen our understanding of jet quenching mechanisms (cf. e.g. [10, 11,12], and [13] for a recent review on jet substructure observables in heavy-ion collisions). Despite their usefulness, a common limitation of many methods developed for pp is that they are based on the assumption that QCD emissions are angular ordered, a property that holds in vacuum (see, e.g. [14]). In the presence of a medium, the phase space opens up to allow anti-angular ordered emissions [15,16]. It is thus natural to expect that reclustering algorithms different from C/A could be more suitable for jet quenching studies.
In this work, we explore a range of generalised k T jet algorithms, with a particular focus on the formation time algorithm, hereon denoted as τ algorithm, obtained by setting p = 0.5. In this case, Eq.(1) reduces approximately to where θ ≈ ∆R ij /R. The formation time of an emission, τ f orm , is the time it takes to behave as an independent source of radiation, and is given by [14]: .
The energy and virtuality of the incoming parent parton are denoted by E and Q 2 , the fraction of energy transported by the extra radiation by z, and θ 12 is the emission angle between the two outgoing partons. The rightmost expression of Eq.(4) is only valid in the high energy limit. If we focus further on the collinear and soft limits, we obtain the parametric inverse of Eq. (3): This manuscript is organised as follows: after general considerations on how jets are generated in Monte Carlo event generators (Section 2), we demonstrate that unclustering with p = 0.5 can be used to extract information about formation times from reconstructed jets in pp and PbPb collisions (Section 3). Section 4 shows an example of an application of the proposed τ algorithm in heavy-ion studies. The final conclusions appear in Section 5. Several details and further comparisons between the different reclustering algorithms are not included in the main manuscript, but can be consulted in the appendices (A to C).

Jets in Monte Carlo event generators
In the absence of a background medium, the scale evolution of jets is calculable in perturbation theory and included in Monte Carlo event generators in the form of parton showers. These generate the additional radiation through an iterative procedure, where the emissions are strictly ordered in the evolution variable. Parton showers are typically implemented in the leading logarithmic approximation. To this accuracy, the evolution variable is not uniquely determined, and a whole class of variables is allowed. The most commonly used are transverse momentum (used for instance, by PYTHIA 8 [17] and SHERPA [18]), angle (HERWIG [19]) and virtuality (PYTHIA 6 [20] and JEWEL [21]). Inverse formation time also belongs to the class of possible evolution variables but is not widely used. A consequence of ordering in a variable other than (inverse) formation time is that the emissions are not necessarily ordered in time. In a virtuality ordered shower, like JEWEL, both energy and virtuality in Eq.(4) decrease in each emission step, but the ratio does not necessarily decrease as well. For typical kinematics, however, the occurrence of emissions un-ordered in time is rare. This justifies the common practice in jet quenching studies of using a transverse momentum or virtuality ordered parton shower with the assumption that the emissions are, in addition, ordered in time. This is the case of JEWEL, a medium modified parton shower that includes the concept of formation time to perform a veto on subsequent emissions.

Event samples
To understand how a change in the p parameter can be more sensitive to in-medium effects of subsequent QCD emissions, we will use two Monte Carlo event generators: PYTHIA 8 (v8.2.35, tune 4C [22]) to check how these jet algorithms correlate to the vacuum shower, and JEWEL, our reference for jet quenching effects. Anti-angular ordered emissions is a feature still missing in most of jet quenching Monte Carlo event generators. In JEWEL, for instance, emissions inside the medium can populate the full accessible phase space, i.e. they are neither angular nor anti-angular ordered. As a first case study to identify the correlation between the emission pattern and the reclustering history, we find JEWEL to be the best choice. We used an unofficial version of JEWEL , based on v2.2.0, that provides the output of the parton shower history. Both generators were set to produce dijet events at a centre-of-mass energy √ s N N = 5.02 TeV. For the medium-induced effects, we restricted ourselves to the medium toy model provided within this event generator (Bjorken model [23] for a boost-invariant longitudinal expansion of an ideal quark-gluon gas), with an initialisation time set to τ i = 0.4 fm/c and initial temperature T i = 0.44 GeV. These values are known to provide an R AA ≃ 0.4, in agreement with current experimental observations at the same centre-of-mass energy [24,25,26] (see also section 4). Two PbPb event samples are used: one with recoils from elastic scattering to include the effects of medium response, and one where these particles are discarded from the final event. When recoils are included, the subtraction of thermal momenta is performed using the new constituent subtraction algorithm [27], which is particularly suitable for jet substructure studies. Unless noted otherwise, the JEWEL results shown are without recoils, so as not to complicate the discussion.

Vacuum parton shower
From hadronic PYTHIA 8 dijet events we reconstruct anti-k T jets with R = 0.5. After identifying the leading jet with a minimum transverse momentum p T,min = 300 GeV that falls within |η jet | < 1.0, we recluster the jet particles with the generalised k T algorithm, for several p values (setting R = 1.0 for the reclustering). These steps are performed within FastJet v3.3.0 [28].
We then iteratively uncluster the reclustered jet, following the leading branch. At each step, we evaluate the formation time as provided by the rightmost expression in Eq.(4), using the information from the two obtained subjets 1 . This procedure yields the unclustering history. The obtained sequence is not necessarily ordered in formation time, when a distance measure different from (inverse) formation time is used for the reclustering. In practice, however, the first unclustering step will usually be the one with the shortest formation time. The probability that this is not the case was found to be a few percent in this study (as expected, it is larger for C/A than for the τ algorithm 2 . Therefore, we always assume that the first unclustering step has the shortest formation time (instead of walking through the whole sequence to find the step with the shortest formation time).
To correlate the unclustering sequence with the parton shower, we identify the initial parton (coming directly from the matrix element) that generated this jet 3 and, following the leading branch, calculate in each step the splitting kinematics as the rightmost expression in Eq. (4). We also studied the impact of applying the high-energy limit, for which we refer the reader to Appendix C.
The correlation of the obtained τ f orm between the first unclustering step and the first parton shower emission in PYTHIA 8 is shown in Fig. 1. The top (Fig. 1a) and bottom (Fig. 1b) panels are for jets reclustered with the C/A (p = 0) and τ algorithm (p = 0.5), respectively. Two features are present: (i) the diagonal 1 Note that the angle θ 12 is evaluated as the angle between the three-dimensional momentum vectors of the two subjets. 2 It is not zero even with the τ algorithm because its distance measure equals the inverse of Eq.(4) only in the soft and collinear approximation. 3 To overcome the difficulties in following the parton shower history when hadronisation effects are taken into account, we identify all the final state particles produced by the two ancestors. We reclustered them with anti-k T jet algorithm and R = 1. The ancestor that generated the jet is considered to be the one whose jet axis is within ∆R < 1.0. For more information, we refer the reader to Appendix A. elements that represent the events in which the correlation is positive (true) (ii) the vertical elements that correspond to large angle emissions that are not captured by the final jet, which leads to a mismatch between the first emission of the parton shower and the first emission captured by the reconstructed jet. Naturally, this vertical component is reduced if we increase the jet radius.  Without any jet grooming technique, the correlation between the unclustering and the parton shower is more pronounced for the τ algorithm. However, this is not an entirely fair comparison, as the C/A algorithm is particularly sensitive to uncorrelated soft large-angle emissions, but is able to provide a proxy for the Altarelli-Parisi splitting function if a SoftDrop procedure is used. Following [8], we apply such a condition to select only emissions (during the parton shower) and/or unclustering steps (during the jet unclustering process) with:  Compared with the results without SoftDrop (Fig. 1), there is a significant decrease of the dispersion in both cases. The vertical events aligned at τ P artonShower f orm ≃ 0 are also drastically reduced. These correspond to relatively soft emissions during the parton shower that would fall outside of the jet cone. With the introduction of a z cut , these emissions are discarded and the match between parton shower and unclustering history improves. Overall, the correlation factor increased from 0.26 (C/A) and 0.38 (τ ) to 0.65 (C/A) and 0.66 (τ ). As such, the use of a SoftDrop procedure is recommended in order to increase the correlation between the jet clustering history and the parton shower.
The effect of SoftDrop grooming is also reflected in the overall τ f orm distribution. In Fig. 3, we show log(τ f orm ) of the first emission as obtained from the PYTHIA 8 parton shower. The blue line represents the distribution when grooming is introduced, and the green when z cut = 0.0. Without grooming, the τ f orm distribution is dominated by (early) large angle and soft emission. With z cut = 0.1, these emissions are discarded and, consequently, the selected first emission has a significantly larger τ f orm . In the following, we will keep If the correlation is successful, the resulting distribution should be narrow and peaked at zero. As an illustration, we show, in Fig. 4, the distribution for the difference defined by Eq. (7) obtained for the first emission/unclustering step with the C/A algorithm. From these distributions, we calculate the median Q 2 , and the Q 1 and Q 3 quartiles as a measure of the width of distribution. These results are collected in Fig symmetric and narrow distribution. As we proceed further through the primary parton shower branch/unclustering history, we start to see deviations in the median value. In particular, the second emission only has a peaked distribution around 0 for p ≃ 0.5. It is also substantially broader, even with SoftDrop grooming. These effects are also visible when considering all emissions along the primary branch. However, we do not see an increasing effect with respect to the second emission, as this distribution is dominated by the first emission (many jets consist of only two subjets after grooming and thus have only one unclustering step). In addition, the average formation time also increases for the second splitting, naturally leading to a larger interquartile range 4 (IQR = Q 3 − Q 1 ). When comparing the relative differences instead of absolute differences the resolution for the second splitting is identical to the first splitting. Nonetheless, in this manuscript, we will restrict ourselves to illustrate an application of using the first unclustering step as a tool for jet quenching studies.
To understand if the apparent success of the correlation of the τ algorithm with the parton shower history is not dominated by details of the Monte Carlo event generator, we repeated the exercise in JEWEL without medium effects (see Fig. 6). The magnitude of the All Emissions 1st Emission 2nd Emission  7)) obtained for different p parameters defined in the generalised k T jet algorithms for the first (orange), second (green) or all emissions along the primary branch (purple) in JEWEL+PYTHIA (pp) events. The asymmetric error bars correspond to +Q 3 and −Q 1 quartiles. The right panel shows a zoom of Q 2 alone. observed shifts and dispersions is similar to the ones obtained in PYTHIA 8 . The main differences are observed for large values of p, where the shift in the distribution of the second emission is more compatible with zero. While the results differ in details, the conclusions are still qualitatively the same: intermediate values of p ∈ [0, 1] yield a better correlation between jet algorithm and parton shower history, with the τ algorithm (p = 0.5) showing the best performance. In particular, it has the most symmetric distribution, which is important for obtaining an unbiased estimate of the formation time. Such an outcome was expected as the τ algorithm uses (inverse) formation as distance measure in the clustering. The fact that the findings are so similar for two rather different parton showers (in particular with different ordering variables) gives us confidence in the robustness of the procedure and why thus turn to study its performance in central PbPb collisions.
The IQR extracted from the inclusive ∆τ f orm distribution does therefore not represent an average resolution.

Medium-modified parton shower
We now include jet quenching effects, and we turn our attention to the results provided by JEWEL with the medium model settings as described in section 2.2. To illustrate the change in the ∆τ distribution, we show, in Fig. 7, the results for JEWEL+PYTHIA (pp) (in blue) and JEWEL+PYTHIA (PbPb) (in green), for the first emission 5 /unclustering step obtained with the τ algorithm. An immediate consequence of including medium effects is the overall smearing of the ∆τ distribution. This is due to further interactions of the partons coming from the first splitting in the medium. They lead to a deterioration of the correlation between the partons' and the subjets' momenta, and limit the feasibility of reconstructing the kinematics at the splitting vertex from the final state particle distribution. Energy loss (energy transported to large distances) and transverse momentum broadening are thus a limiting factor in the endeavour to estimate the formation time of the first splitting in the presence of a background medium (this comes on top of effects like vacuum-like emissions and contamination from initial state radiation governing the resolution in pp). In addition, the average formation time increases in PbPb compared to pp (see discussion in section 4), and the correlation between reconstructed and true formation time deteriorates for larger τ f orm (cf. fig. 2).
The characteristics of the ∆τ distributions obtained from the unclustering history and the parton emission history in JEWEL (PbPb) for different exponents of the generalized k T algorithm are shown in Fig. 8. As expected, the smearing observed in Fig. 7 is also present All Emissions 1st Emission 2nd Emission  7)) obtained for different p parameters defined in the generalised k T jet algorithms for the first (orange), second (green) or all emissions along the primary branch (purple) for JEWEL (PbPb) events. The asymmetric error bars correspond to +Q 3 and −Q 1 quartiles. On the right, a zoom of the Q 2 alone is shown.
for subsequent emissions and across all p exponents. Nonetheless, the same qualitative features appear as in Fig. 5: p ≃ 0.5 provides the best correlation (less biased and narrower ∆τ ) between parton shower and jet unclustering history.
When recoils are considered, part of the energy lost by the jet is recovered (in the form of additional soft activity inside the jet cone). Including medium response thus mitigates the mismatch between subjet kinematics and kinematics at the splitting vertex. At the same time, it also introduces as additional smearing: as the additional activity from medium response has a broad distribution, subjets will, in general, also include some level of activity from other sources than its parent parton. Overall, there remains a net improvement of the correlation reducing the width of the distribution, independently of the p exponent in the generalized k T algorithms (see Fig. 9). This effect is more noticeable in the second emission: being softer, the effect of the medium recoils is, proportionally, larger than for the first unclustering step. Given that the nature of medium recoils is still under debate, we will focus on JEWEL without considering this medium recoil component, thus simplifying following discussions.
We also investigated the role of the jet radius and Soft-Drop z cut parameter, with focus on C/A and τ clustering algorithms (see Appendix B and Appendix C for the corresponding results without the high-energy approximation). As expected, an increase of the jet radius increases the correlation and reduces the smearing of the ∆τ distribution, for both first and second emission. Large radius jets, however, pose a challenge in a heavy-ion environment due to the high event activ-  Fig. 9 Median value of the ∆τ distribution (Eq. (7)) obtained for different p parameters defined in the generalised k T jet algorithms for the first (orange), second (green) or all emissions along the primary branch (purple) for JEWEL (PbPb) events with recoils. The asymmetric error bars correspond to +Q 3 and −Q 1 quartiles. On the right, a zoom of the Q 2 alone is shown.
ity [26,29]. As for the effect of varying z cut , there is a trade-off between reducing contamination from the uncorrelated soft activity and discarding correlated structures. Overall, we found that z cut ≃ 0.1 seems to be a good compromise.

Summary
To present a summary of the main properties of the ∆τ distribution, we focus on the first and second emissions of the leading anti-k T jet with R = 0.5 and p T > 300 GeV, reclustered with the generalized k T algorithm and Soft-Drop procedure (z cut = 0.1, β = 0).
The median of the distribution for the first emission is compatible with 0 and so we only focus on the interquartile range (IQR = Q 3 − Q 1 ). To provide a better visualisation given the different magnitudes between pp and PbPb, we identify the minimum IQR value obtained from the different p exponents used in the generalized k T algorithm and subtract this amount, such that the minimum is zero for each of different samples. These are shown in Fig. 10 for PYTHIA 8 (purple/circle markers), JEWEL (pp) (orange/square markers) and JEWEL (PbPb) (green/triangle markers). All models consistently indicate that p = 0.5 yields the smallest dispersion. When C/A is used instead, the width of the distribution increases by an additional 1 fm/c (pp) up to 2.2 fm/c (PbPb).
For the distribution of the second emission we focus on the median, as there are significant differences between the clustering algorithms. Again, we take the difference with respect to the minimum value obtained for the different p exponents. The results for the different Monte Carlo parton showers are shown in Fig. 11. In this case, p = 0.5 seems to be preferred by JEWEL (PbPb), but vacuum parton showers do not agree on the exponent that maximises the correlation: with PYTHIA 8 intermediate values yield a better correlation, while JEWEL (pp), which is based on PYTHIA 6, shows a less biased distribution when larger values of p are used instead. These differences might be related to the details of the implementation of the parton shower, including small violations to the time ordering of the emissions. Nonetheless, the differences between the two vacuum parton showers are up to 0.5 fm/c for the k T algorithm. As remarked earlier, we here focus on the first splitting. These results suggest that the use of unclustering tools others then standard C/A for jet quenching studies might potentiate an accurate evaluation of the QGP properties. In particular, the use of the τ algorithm for reclustering allows to directly extract an estimate of the formation times in a more meaningful way. The use of formation time from splitting processes is already under exploration by the community [30,31]. It will, for sure, bring a novel approach to probe the QGP evolution, following recent efforts [32,33]. By using the τ algorithm for reclustering, we can increase the precision of such studies and access more differential QGP timescales. While such a feasibility study will be left for a future communication, we will proceed with an example to compare the use of τ and C/A for jet quenching studies.

Comparison of τ and C/A reclustering algorithms for jet quenching studies
In Fig. 12, we compare the formation time of the first groomed emission obtained in JEWEL pp (blue) and PbPb (green) from the parton shower. There is a clear increase in the average formation time in PbPb with respect to pp. Several effects contribute to a change of the resulting τ f orm distribution when medium-induced effects are considered. The most immediate consequence -energy loss (both elastic and inelastic) -leads to a decrease of τ f orm . In the case of elastic energy loss this occurs through the degradation of the energy in eq. (4) (the virtuality does not change), while in the case of inelastic energy loss it happens via the occurrence of an extra splitting (with formation time shorter than the vacuum-like splitting). However, both effects are small for the first emission. Elastic energy loss is not important since the parton is very energetic, and inelastic scattering is extremely rare due to the very short formation time of the first vacuum-like splitting. On top of energy loss, QGP interactions also induce jet collimation, an effect that has already been observed in other observables [34,35,13]. The surviving jets in PbPb collisions are biased towards having a hard (vacuum-like) fragmentation pattern, which is correlated to larger formation times. This effect is dominant for the first emission and contributes to the overall shift of the formation time distribution observed in Fig. 12.
We now proceed to evaluate the nuclear modification factor of leading jets provided by JEWEL, but extending this study to lower transverse momentum (p T > 100 GeV). Using C/A and τ to recluster the leading jet, we identify the first unclustering step and create two populations: early jets, whose first unclustering step has τ f orm < 1 fm/c, and late jets, with a first unclustering step with τ f orm > 3 fm/c. We make this selection in JEWEL (PbPb) and JEWEL (pp) events to obtain the leading jet transverse momentum spectrum in both cases. The corresponding medium-over-vacuum ratio (nuclear modification factor, R AA , for leading jets) is shown in Fig. 13. For reference, we also include the inclusive leading jet ratio, in solid back. The purple lines refer to reclustering with τ algorithm for late (solid line) and early (dashed line) jets, while the orange refers to the C/A algorithm. For reference, we also include the results directly read from the parton shower, in green. There is a clear difference in the leading jet suppression when, instead of using the full sample, we select jets whose fragmentation starts shortly after its production.
These jets are, as expected, strongly suppressed, and both C/A , and τ algorithms provide similar results. Taking the results from section 3, we do not expect to see much deviations between the two. However, as we move towards late times, the two algorithms show some differences. In particular, if we use τ to recluster the jet particles, the obtained R AA is compatible with 1. As discussed earlier, these jets have a hard fragmentation pattern and are therefore not so susceptible to modifications due to medium interactions as those with a soft fragmentation, and thus early first splitting. In particular, the late jets consist of only one effective colour charge with high momentum (in this case ∼ 300 GeV) for the first 3 fm of the evolution. This object loses little energy through elastic scattering, and, when it finally splits, the medium density is already diluted (ǫ < 5 GeV/fm 3 for the medium settings used here and the simple medium model). At relatively low p T , both algorithms yield a similar difference with respect to the Monte Carlo truth (one suppressed, the other enhanced), but at high p T , the results using the τ algorithm approaches the Monte Carlo. This is in line with the observations of the previous section.
When medium recoils are considered, we see the same behaviour, see Fig. 14 (same colors and line settings as in Fig. 13). The early (and inclusive) leading jet R AA are now slightly larger, as part of the energy is recovered by the presence of recoils. Both reclustering algorithms continue to yield the same results. However, for late jets, there are sizable differences between the two reclustering algorithms. The difference between C/A and τ is related to the fact that the τ f orm values extracted with C/A are systematically too small, particularly for large τ f orm (cf. figs 2, 6 and 7). The result obtained with the τ algorithm comes closer to the one extracted directly from the parton shower. This highlights the importance of an unbiased estimate of the true formation time and is a strong argument for using the τ algorithm to obtain the formation times.

Conclusions
In this work, we explore the use of generalised k T algorithms to correlate the parton shower history to the unclustering sequence for the primary branch. We focus on the formation time, τ f orm , to have a space-time picture of the parton shower emissions. This study is based on Monte Carlo event generators with different ordering variables (PYTHIA 8 and JEWEL) for subsequent emissions. The aim was to check the resilience of our results. As exponents for the generalized k T algorithm we focus on p = 0 (C/A) and p = 0.5 (τ ). We found that the formation times extracted with the τ algorithm generally correlate better with the parton shower values, even in the absence of SoftDrop grooming. When a z cut is introduced, both algorithms yield relatively small differences in τ f orm , but C/A shows consistently a broader and more asymmetric distribution. Naturally, the fraction of unclustering sequences that are not ordered in formation time is larger for C/A . This could, to some extent, be overcome by going through the entire sequence to find the unclustering step with the shortest formation time (instead of always taking the first one). While in pp collisions the angular ordered sequence obtained with the C/A algorithm is often preferred due to its similarity with the QCD radiation pattern, this is not necessarily the case in heavy ion collision. As an example, this study shows that to estimate formation times it is advantageous to recluster the jets with the τ algorithm, that interprets the fragment distribution in terms of inverse formation time. As an example, we found that the τ f orm obtained with C/A yields a leading jet nuclear modification factor of around 0.8 for jets whose fragmentation starts after τ f orm = 3 fm/c, while the τ algorithm does not yield any suppression. This observation is consistent with the event generator results for jets with p T > 300 GeV.
Time differential measurements of the QGP are increasingly needed. Jet reclustering tools that allow to increase the precision with which such short timescales can be extracted will undeniably be most helpful in unlocking the use of jets as precise tools for QGP tomographic measurements. In the future, we plan to use these tools to perform a feasibility study in a heavy-ion environment.
Acknowledgments The authors would like to thank to R. Conceição, L. Cunqueiro, G. Milhano

Appendix A: Identification of the jet initiating parton
The identification of the jet initiating particle can, in PYTHIA 8, be done easily by looking at the ancestor list of the jet constituents (denoted as v1 in Fig. 15a  and 15b). The initiating parton is taken to be the outgoing parton of the matrix element that contributes most of the jet constituents. In this way, at parton level the initiating parton is uniquely defined. At hadron level, however, this is not always the case, as hadrons can't always be associated with only one of the initiating partons. It is thus expected that in some events it will not be possible to make a correspondence between the outgoing parton of the matrix element and the final reconstructed jet.
To overcome such difficulty, we adopted the following procedure (denoted as v2 in Fig. 15a and 15b): -Identify the outgoing particles of the matrix element -Tag all final state particles produced by such initiating partons 6 -Use them as input to reconstruct anti-k T jets with R = 1. -Compare the ∆R distance between the leading jet and the two jets obtained in this way. -Identify the initiating parton as the one whose final jet is the closest one in ∆R (with up to ∆R < 1).
In PYTHIA 8 at hadron level this procedure yields a loss of 9% of the selected events (with at least one jet with p T > 300 GeV and η < 1), compared to the 19% when adopting the approach of going through the ancestor list.  We compared the obtained ∆τ distributions, at parton and hadron level, following these two procedures, and the results are statistically consistent (see Fig. 15). At parton level, the number of lost events is about 1% for both v1 and v2 2% at hadron level. Both tagging procedures (v1 and v2) fail to identify 1% of the events at parton level and 2% of the events at hadron level.
In JEWEL , however, such procedures are not possible, as the hadronization vertex links all particles of the hard scattering event. Therefore, we instead reconstruct R = 1 jets from all final state partons (instead of hadrons) originated from the matrix element. This reduces the effects of the hadronization and increase the success rate of identifying the jet initiating parton. In addition, for JEWEL (PbPb), we also follow the outgoing medium particles that participate in the elastic scattering process. Overall, these modifications intro-duce a loss of 2% of the selected events in JEWEL (pp) and JEWEL (PbPb). As such, we expect the obtained final distributions of ∆τ to be unbiased by the tracing procedure.
Appendix B: Correlation as a function of the jet clustering radius and z cut The equivalent of Fig. 8, but fixing p = 0 (C/A) or p = 0.5 (τ ) while varying the jet radius, R, and the SoftDrop parameter z cut , are shown in Fig. 16 and Fig. 17.

All Emissions
1st Emission 2nd Emission   7)) for different jet radii for the first (orange), second (green) or all emissions along the primary branch (purple) in JEWEL (PbPb). The asymmetric error bars correspond to the +Q 3 and −Q 1 quartiles. On the right, a zoom of Q 2 alone is shown.
As expected, an increase of the jet radius allows to capture more particles coming from the parton shower, leading to a more complete reconstruction of the jet fragmentation pattern. The matching between parton shower and reclustering algorithms is thus improved, up to R = 1.0. From there, the results seem to stabilize, and there is no significant gain from increasing the jet radius further. We would like to note that increasing the jet radius in PbPb events is experimentally challenging given the large multiplicity [26,29]. Moreover, one would also expect an increasing contamination from initial state radiation (ISR). In our case, we do not expect a large effect from ISR given the central rapidity range (|η| < 1), but this effect should be taken also into account if a larger η range is desired.

All Emissions
1st Emission 2nd Emission    7)) obtained with different values of the SoftDrop parameter z cut for the first (orange), second (green) or all emissions along the primary branch (purple) in JEWEL (PbPb). The asymmetric error bars correspond to the +Q 3 and −Q 1 quartiles. On the right, a zoom of Q 2 alone is shown.
The dependence on z cut is the less obvious. By increasing this parameter, the resulting ∆τ distribution for the first emission is shifted towards values near zero, but the overall dispersion is increased. Two competing effects are at play. Introducing a z cut allows to select only hard QCD emissions, that are more easily matched between clustering and parton shower history (the correlation between the two procedures improves, as observed in section 3). On the other hand, even though the transverse momentum selection is made on the ungroomed jet p T , grooming also introduces a selection bias favouring jets with a harder fragmentation, and thus larger formation times. When there is a mismatch between the estimated τ f orm obtained from reclustering and parton shower history, the magnitude of the possible ∆τ is thus also increased. Overall, we find that a z cut broadens the distribution, but also shifts it to ∆τ ∼ 0.
Appendix C: Effect of the τ f orm definition for jet quenching studies When calculating τ f orm from Eq. (4) (rightmost expression) an approximation is used. We found that if instead we use τ f orm ≈ p T,1 + p T,2 p T,1 p T,2 ∆R 2 (C.1) where p T,i refers to the transverse momentum of subjet 1 and 2 and ∆R 2 to the (y, φ) distance between the two, the results are very similar to the ones obtained from Eq. (4). Both expressions assume high-energy limit, and neglect the virtuality (mass) of the parton (jet). In the parton shower, the assumption of massless particles might be unrealistic, in particular, if we focus on the first emission. Therefore, we have studied the impact on our distributions if, instead, we use where m 1 and m 2 now refers to the virtuality (mass) of the 2 outgoing partons (subjets), and p i to their respective 4-momenta. Focusing only on the first emission, we show in Fig. 18 the comparison between the ∆τ distribution when using τ f orm as defined by equation (4) (named τ 1 , in purple) and when we use the definition in (C.2) (named τ 2 , in orange). For both the τ algorithm is used in the reclustering procedure. By comparison, τ 2 (eq. (C.2)) produces a more asymmetric distribution, typically yielding a smaller τ with respect to the parton shower. This is not entirely unexpected, since τ 2 involves the sub-jet masses, a quantity that is known to receive large nonperturbative corrections. The sub-jet invariant mass obtained at hadron level is therefore not a reliable estimate of the parton level quantity, i.e. the parton's virtuality. Since the virtualities are typically small compared to the parton energy, it is safer to assume the sub-jets to be massless. Under this assumption equation (C.2) reduces to equation (4). The same effect is visible when using the C/A algorithm (Fig. 19), but the shift seems to be less pronounced. Based on these results, we con-