UvA-DARE (Digital Academic Repository) Renormalization group flows for track function moments

,


Introduction
The characterization of energy flow within jets, colloquially known as jet substructure, provides new ways to study QCD and search for potential new physics at the LHC [1,2].The remarkable advances in this area in the last decade have primarily focused on the calculation of infrared and collinear (IRC) safe observables that can be computed within perturbative QCD, up to power corrections.The famous theorems of Kinoshita, Lee and Nauenberg [3,4] state that this is only possible if one is completely inclusive over hadron species.As a consequence, such calculations can only be used to describe observables constructed from energy flow information, disregarding all the interesting information contained in other particle properties.Theoretically, these observables are therefore (combinations of) correlation functions of energy flow operators, ⟨E(⃗ n 1 )E(⃗ n 2 ) There is significant motivation to go beyond this energy flow paradigm, both for allowing more detailed tests of QCD, and for sharpening our tools in new physics searches.Such observables are inherently non-perturbative, as they require knowledge of the spectrum of hadrons in the theory.For example, at the LHC, many precision jet substructure measurements are made using tracks (charged particles), due to the improved angular resolution of the tracking system.This sensitivity to hadronization can of course also be viewed as a positive if the goal is to understand features of the hadronization process.For example, the study of energy flow on charged or strange particles provides insight into how these quantum numbers evolve in the confinement process.
The departure from IRC safety should not be done arbitrarily, and in particular, one should attempt to maintain the wealth of theoretical structures and advances of perturbative quantum field theory, but generalize them to a wider class of observables.In ref. [5], building on [6], it was shown that the natural way to extend the space of IRC safe observables to incorporate particle species information is to consider correlations of energy flow on subsets of particles.These are defined theoretically by considering an energy flow operator on a subset R of particles, E R (⃗ n 1 ), and enable a much more general class of correlations to be studied, where in general the subsets, R i , are distinct.As we will discuss, these observables exhibit a clean factorization into a non-perturbative component, and a perturbative component.The perturbative component shares many of the features of the standard energy correlators, and in particular can be computed at high perturbative orders using well-developed techniques from perturbative quantum field theory.
Although the correlators ⟩ cannot be directly computed in perturbation theory, they can be matched onto the standard energy flow correlators using nonperturbative track functions [7,8].These track functions were introduced to describe the fraction of energy deposited into charged hadrons from a perturbative quark or gluon, however, they can trivially be generalized to the study of any other quantum number.Unlike standard fragmentation functions, track functions incorporate correlations between particles, arising from the fact that quarks and gluons can fragment into an arbitrary number of charged hadrons.As such, their evolution with scale is substantially more complicated, since all the correlations mix under evolution.
In ref. [6], it was shown that, by restricting to correlation functions of energy flow measured on tracks, one is only sensitive to low moments of the track functions.These characterize the fluctuations in the hadronization process. 1 To describe N -th order fluctuations requires only a finite set of operators, which mix under renormalization.Furthermore, the full track function distributions seem well-described by a truncated Gaussian, whose form is fixed by the first two moments.In ref. [5] it was shown that energy conservation places severe constraints on the RG evolution of the fluctuations, fixing the evolution of the first three moments to be DGLAP, up to corrections proportional to powers of ∆ = T q (1) − T g (1).For track functions describing the production of electrically charged hadrons in QCD, ∆ ≪ 1, effectively suppressing these contributions by an order in the perturbative expansion.At the fourth moment and beyond the fluctuations in the hadronization process exhibit non-trivial RG flows describing the mixing between different cumulants, for example κ(4) and κ(2) 2 .
In this paper we discuss in detail the structure of the RG for the moments of the track functions.In dimensional regularization, the corrections for the track functions are scaleless thus linking the evolution (UV poles) and the IR poles needed for incorporating track functions in calculations.We derive general constraints on the structure of the evolution that hold to all orders in perturbation theory, and in generic theories.In QCD, we then analytically compute the first six moments at next-to-leading order (NLO), and study the structure of their RG flows, which exhibit interesting mixing.For the first three moments the mixing terms are all suppressed by powers of ∆ and smaller than the NNLO corrections, allowing us to extend our calculation to this order.We also argue, that due to the nonlinear nature of the track function evolution, it exhibits a UV fixed-point where the track functions become a delta function.Our explicit results enable the calculation of jet substructure observables sensitive to up to six point correlations in energy flow on tracks.
While the primary motivation for this work is practical, namely enabling higher point correlators to be precisely measured at the LHC, the study of track functions is also of more formal theoretical interest.Track functions, and related multi-hadron fragmentation, are intrinsically Lorentzian observables whose RG evolution goes beyond standard DGLAP evolution.Although there has been significant recent progress in understanding certain classes of Lorentzian operators using lightray operators [9], this has primarily been restricted to operators on the leading Regge trajectory (which includes DGLAP).Understanding how the more general class of track function observables fits into this picture is interesting, and could lead to a better understanding of the analytic structure of Lorentzian observables in conformal field theories (CFTs).While we will not address this issue directly in this paper, our perturbative calculations provide important theoretical data for future investigations.
The outline of this paper is as follows: We discuss the flow of energy on subsets of particles in sec.2, motivating the study of moments of track functions.In sec. 3 we review the fieldtheoretic definition of track functions, and derive all-orders constraints on the renormalization group evolution of their moments.We then restrict to NLO, and derive the specific constraints both for a pure gluon theory, as well as for QCD.In sec. 4 we present results for the first six moments of the track functions at NLO, and describe the techniques used in the calculation.More details of the calculation for Pure Yang-Mills are given in app.A, which include results up to ninth moment, and the time-like splitting functions entering our results are collected in app.B. In sec.5 we numerically study the structure of the RG flows.We first show that in QCD, ∆ ≪ 1, allowing us to extend our results for the evolution of the first three moments to NNLO.We then study the importance of non-linearities in the evolution of the fourth and fifth moments.We conclude in sec.6.

Energy Flow on Tracks and Track Function Moments
To motivate the study of track function moments, we begin by reviewing the natural generalization of the study of correlations of energy flow, to the study of energy flow on subsets R of particles.Figure 1: (a) For a standard dijet event shape observable, which constrains the phase space of all emissions, a separate track function is needed for every emission, leading to a complicated structure of the hadronization process.(b) For energy correlators, matching can be performed at the level of the detectors, instead of for each parton.Since the number of detectors is fixed this leads to a much simpler description of the transition from quarks and gluons to hadrons.
Here we will see that the non-perturbative information required for this extension is precisely the moments of track functions, motivating our focus on these moments.
Energy flow in final states is characterized by the energy flow operator [9-16] The canonical observables of the theory are the k-point correlation functions These generalize the original two-point correlator introduced early on in the QCD literature [17].
Although these observables appear similar to more standard jet observables, which are typically called "jet shapes", they are in fact quite different.Jet shapes constrain radiation about some underlying hard process, can be thought of as infrared and collinear safe resolution variables for an S-matrix element of quarks and gluons.On the other hand, the correlation functions are statistical correlators defined as an ensemble average, and do not constrain the emitted radiation.While these correlators have been well studied in the formal CFT literature, that they can be useful phenomenologically to systematically probe the structure of QCD was emphasized in ref. [6].
The energy correlators are simpler perturbatively, which has enabled a number of remarkable calculations in both QCD [19,20] and N = 4 SYM [18,21].However, for phenomenological applications to QCD, it is perhaps their non-perturbative simplicity that is even more important, due to the poor current understanding of the hadronization process in QCD.Standard jet or event shape observables are sensitive to the complete structure of emissions.This makes their extension to charged particles (or other subsets R of particles) extremely complicated, since it requires a description of the hadronization process for every single perturbative particle.This is illustrated in fig.1a.Furthermore, in addition to having additional track functions at each perturbative order, the observable also depends on the complete functional form, T a (x), of these non-perturbative functions.On the other hand, for correlation functions of energy flow operators, the fragmentation process should be thought of as a matching between detector operators in the perturbative and non-perturbative theory.Since the number of detectors is fixed (and in practical applications only low numbers of detectors are considered), this leads to a simple theoretical description of the fragmentation process, that is unchanged order by order in perturbation theory, see fig.1b.It is this simple property of the energy correlators that allows them to be naturally extended to a description of energy flow on subsets of particles.
We now formalize this in a factorization theorem involving moments of track functions.This will motivate the study of the renormalization group structure of these moments, which will be the focus of the remainder of this paper.To understand the energy correlators on tracks, we begin by introducing an energy flow operator that only measures energy flow on a restricted set of states, E R .This is a fundamentally non-perturbative object, which does not admit a perturbative expansion about free asymptotic quark and gluon states.This restricted energy flow operator admits an OPE onto partonic energy flow operators, E R (⃗ n 1 ) = T q(1)E q(⃗ n 1 ) + T q (1)E q (⃗ n 1 ) + T g (1)E g (⃗ n 1 ) . (2. 2) The matching coefficients are given by first moment of the track function T a (1), describing the average momentum fraction of the subset R, whose formal definition and RG structure will be given in the next section.(Note that track functions can differ between quark flavors, which we ignore here for notational simplicity.)To study multi-point energy correlators on tracks, one will therefore need to perform the perturbative calculations of the matrix elements These are more general than what has been studied in the literature, but the same calculational techniques can be used, as will be discussed in sec.4.
We are now able to present the general form of the factorization formula for a k-point correlator in terms of these partonic correlators and moments of track functions (2.4) The contact terms arise when any two detectors are in the same direction, introducing dependence on higher moments of the track functions.We will now explicitly show the structure of the contact terms for the two-and three-point correlator.For the two-point correlator, we have while for the three-point correlator, we have The extension to higher point correlators should be clear.These contact terms introduce dependence on higher track function moments T a (n).The precise operator definition of the corresponding lightray operators, , will not be important here, but in perturbation theory these simply weight the state by E n , where n is the number of 1 in the exponent.The precise notation is chosen due to their relation to multi-hadron fragmentation functions.
One appealing aspect of this factorization formula is that for an N -point correlator, it contains a finite sum over the different track function structures.This structure is fixed by the properties of the detectors, and independent of the order in perturbation theory, as visualized in fig.1b.This follows the general philosophy arising from CFTs, namely that one should study the space of detectors rather than the states, which leads to significant simplifications here.

Track Function Moments and their Renormalization Group Evolution
Having shown how moments of track functions naturally appear in the study of energy flow, in this section we study in detail their renormalization group structure.

Definition and Sum Rules
The track function describes the momentum fraction x of an initial parton i that is converted to a subset R of the final-state hadrons specified in terms of some particular quantum number, e.g.charge, strangeness, etc.Its definition in terms of a matrix element in quantum field theory is in light-cone gauge given by [7,8] In general covariant gauges, Wilson lines are required to maintain gauge invariance, as is standard for fragmentation functions.The Fourier transform of y + fixes the large light-cone momentum of the initiating field to be k − , and the y ⊥ -integral sets its transverse momentum to zero.The delta function encodes the measurement of the momentum fraction x of the subset R of the final-state X.Finally, the matrix elements encode the probability of a quark or gluon to produce a finalstate X, averaged over its color and spin (with d the number of space-time dimensions, used as regularization).
We will often work in terms of the moments of the track functions, defined as Note that this differs by one unit from the standard convention, which is why the evolution of T (n, µ) will involve the DGLAP anomalous dimensions γ(n + 1) in the standard convention.The zeroth moment satisfies the sum rule implying that the track function is normalized.

Comparison to Fragmentation Functions
The difference between the definition of the track function in eq.(3.1) and the fragmentation function so instead the momentum fraction x of a hadron h (e.g.h = π + ) is measured.Because a single parton can produce multiple hadrons, the fragmentation function is not normalized, in contrast to eq. (3.3).Instead, it satisfies the momentum sum rule where the sum on h is over all hadron species.Note that this is consistent with eqs.(3.3) and (3.4) because In grouping h and X ′ together in X, the factor P − h /k − is necessary to get the correct symmetry factor, because X ′ may also contain another hadron h.This is discussed in sec.2.5 of ref. [37].
The first moment of the track function and fragmentation function are related However, for the second moment where D a→h 1 h 2 (1, 1, µ) is a moment of the dihadron fragmentation function.This arises because x = i x i where x i is the momentum fractions of the i-th hadron in R, and x 2 = i x 2 i + i̸ =j x i x j .(For the corresponding discussion in the context of jet charge, see ref. [38].)This can be extended to the n-th moment of the track function, which involves n-hadron fragmentation functions, clearly demonstrating that the track function is sensitive to correlations between finalstate hadrons.

Renormalization Group Evolution and Shift Symmetries
The track function evolution has the following general form where we suppressed the argument µ for brevity.There is a sum over all possible splittings of a parton a into partons a f with momentum fractions z f , and for each of these parton there is a track function T a i .The total momentum fraction x is obtained by summing over the x i of these partons, which are rescaled because these fractions are with respect to the parton a i who carry a momentum fraction z i of the initial parton a.The sum on N goes up to the order α N −1 s that one is working to in perturbation theory.E.g. at order α 2 s we need at most N = 3, corresponding to 1 → 3 collinear splittings.The explicit expression for P is only known at order α s , for which N = 2.
We note that this evolution equation is invariant when the arguments of all track functions are shifted T a (x) → T a (x + b) and T a i (x i ) → T a i (x i + b).This follows because x − i z i x i = (x + b) − i z i (x i + b) due to momentum conservation i z i = 1.Track functions must satisfy 0 ≤ x, x i ≤ 1, and thus for a generic track function this shift cannot physically be performed.However, the evolution equation is independent of the functional form of the track function, so that one can choose to consider a compactly supported track function on which the shift does make physical sense.This allows shifts to be used to constrain the form of the evolution.
Converting eq.(3.9) to moment space for integer n, we can use the multinomial expansion to obtain (3.10) The sum of the moments of the track functions on the right-hand side must equal n, i.e. i m i = n.
The aforementioned shift symmetry of the evolution is particularly convenient for moments: Explicitly, for the first few moments, (3.12) In the next subsections we will work out the consequences of this, starting with the case of a pure Yang-Mills theory that allows us to ignore flavors.The evolution of the fragmentation function can be derived from the same P in eq.(3.9) In moment space this becomes Here we have used the standard conventions for the timelike twist-two spin-n, anomalous dimensions, γ(n).A comparison of eqs.(3.10) and (3.14) reveals that the coefficient of the anomalous dimension of T a (n) involving T b (n) is the same as that entering in the evolution of the moment D a→h (n) of the fragmentation function,

Constraints from Shift Symmetry: Pure Yang-Mills theory
We will now demonstrate how the shift-symmetry determines the structure of the evolution equation for a pure Yang-Mills theory. 2 From the form of eq.(3.10) we know that etc. Since we have only a gluon, we suppress flavor labels.The notation γ 1 , γ 2 , γ 11 , . . .for the anomalous dimensions is only used in the pure gluon case described here and in app. A. From the perspective of the shift symmetry alone, these anomalous dimensions are arbitrary.We will later relate them to the timelike twist-2 spin-n anomalous dimensions γ(n) (note the differing notation).
Applying the shift to these equations, we obtain which leads to and thus γ 1 = 0 in this case (this is not true when there are other parton species), as well as implying A more economical approach to deriving these equations is to directly use shift-invariant central moments where the average ⟨x⟩ is simply the first moment T (1, µ).Note that this can simply be thought of as a change of basis.Now we immediately have since no other terms can appear on the right-hand side.Inserting σ(2) = T (2) − T (1) 2 , we then again obtain γ 11 = −γ 2 .As we will see, in the case of multiple flavors one can form shift invariant first moments, T i (1) − T j (1).Extending this to higher moments, we obtain the general structure of the renormalization group evolution of the central moments of the track functions etc.Because the evolution of T (n) can involve at most 3 track functions at order α 2 s , the form of these equations are further restricted at this order.Thus, up to order α 2 s , This structure for the evolution is fixed entirely by shift symmetry alone.However, this does not fix the values of the anomalous dimensions.To further fix the anomalous dimensions, we note that from their definition, the diagonal anomalous dimensions γ n are related to the timelike twist-2 anomalous dimensions (moments of the gluon fragmentation function), γ gg (n), by These anomalous dimensions are known to NNLO [39][40][41][42].Therefore up to σ 5 all anomalous dimensions are constrained in terms of the DGLAP splitting functions, for σ 6 only one new anomalous dimension needs to be calculated and no new one is needed for σ 7 .Beyond σ 7 , one (or more) new anomalous dimensions need to be calculated for every moment.
An alternate approach is to exploit the symmetry of the matrix elements.This is in practice equivalent to the shift symmetry, though restricted to a specific order in perturbation theory.For example, at order α 2 s for which N = 3, we can express the γ in the equations above to that in eq.(3.10), using momentum conservation z 1 + z 2 + z 3 = 1.Similarly, using the symmetry under permutations of z 1 , z 2 , z 3 .In the final steps we used that under the integral the following identities hold Clearly the use of shift-symmetric central moments is much simpler.

Constraints from Shift Symmetry: Multi-Flavor
Having described in detail how shift symmetry constrains the form of the evolution in the case of a pure gluon theory, we here extend the discussion to the case of multiple parton species, which is needed for QCD.We will consider the case of one quark species and assume that the track functions for quarks and anti-quarks are the same, to keep the discussion simple and highlight the new features.The extension to multiple quarks is straightforward, and our final results do not use this assumption.
The simplifying feature of the pure gluon evolution is that the mean, T (1), is not shift invariant, and therefore cannot appear in the evolution equations.Shift symmetry, combined with the uniqueness of the shift invariant second and third moments, then fixes to all orders in perturbation theory the evolution equations for the second and third moments When moving to multiple flavors there are two new features that appear.The first is a trivial extension, namely that we must extend the evolution equations to be matrix equations in flavor space, as is familiar from DGLAP.Focusing for simplicity on the case of one quark and one gluon, we define as well as the standard matrix of anomalous dimensions The second extension that appears in the case of multiple flavors is a more non-trivial modification, namely the appearance of a new shift invariant quantity, constructed from the difference of first moments.This object can appear in the evolution equations, leading to additional complexity.Focusing on the first five moments, which makes the general structure clear, shift invariance then implies that to all orders in perturbation theory, The presence of ∆ significantly complicates the form of the evolution compared with the pure gluon case, and in particular, the first three moments are no longer uniquely fixed by the shift symmetry.Note that the anomalous dimensions γσ 2 σ 2 , γσ 3 σ 2 and γσ 2  2 ∆ are rank 3 tensors, taking a matrix as input and returning a vector.
The additional complexity arising from the presence of quarks can be thought of in the two different ways discussed in sec.3.4: From the shift-symmetry perspective, the complexity arises purely from the presence of the new invariant ∆.From the perspective of the calculation from matrix elements (discussed briefly at the end of sec.3.4 and made more concrete in sec.4.1.2),the presence of quarks implies that one can no longer symmetrize over the final state particles when using momentum conservation arguments to reduce integrals.The differences that arise from this lack of ability to symmetrize are then captured by powers of ∆.The integrals for these residual ∆-dependent pieces turn out to be simpler to compute.
Despite the fact that the terms proportional to ∆ are not fixed in terms of the DGLAP kernels, we will see that this organization still proves extremely useful, particularly for the case of track functions describing the momentum fraction of charged particles in QCD.In the high energy limit, where the energy cost to produce pions is negligible, one expects that the average properties of the track functions are fixed by isospin, namely T g (1) ≃ T q (1) ≃ 2/3, and ∆ ≃ 0. This intuition is born out by the evolution equation for ∆ in eq.(3.32),where the positivity of γ qq (2) + γ gg (2) drives ∆ → 0 at asymptotic energies.This behavior is already well born out at moderate energies, where one finds the approximate numerical relation ∆ 2 /σ 2 ∼ a 3/2 s , showing that its contribution to the evolution of the second moment is suppressed in the perturbative expansion of the evolution.We will show in sec.5.1, the NLO terms proportional to ∆ in the evolution of the second moment are irrelevant even compared to the NNLO DGLAP corrections.For the third moment, the corrections in ∆ are effectively suppressed by one order in the perturbative expansion.This allows us to extend our results for the first three moments to NNLO, which is the most important practical application of the shift symmetry.
The shift symmetry also forces the evolution of the first moments to be proportional to ∆, namely This result also follows from energy conservation in the one point function ⟨E(⃗ n 1 )⟩, further emphasizing the connection between the shift symmetry and energy conservation.This result shows that the evolution of the first moments of the track functions is numerically suppressed by a factor of ∆/T (1), as compared to the naive expectation.The inclusion of tracks in factorization formulas for energy correlators will therefore have an extremely minor effect, explaining the observation of [6].Finally, one appealing feature of the structure of the equations in eq.(3.32) is that it is known that the eigenvalues of the γ(N ) are positive.This allows us to immediately see that the cumulants (or central moments) of the track functions decay to zero.In the high energy limit, they converge to a delta function with ∆ = 0, which is the unique attractive fixed point of the evolution.The limiting value of T q (1) = T g (1), corresponding to the position of the delta function, is the only nonperturbative parameter that remains.

Track Function Moments at NLO
Having discussed the general structure of the RG evolution of track function moments in sec.3, we now move on to their calculation in QCD.We describe our calculational technique in sec.4.1, and present the full results for the first six moments in sec.4.2.For simplicity, throughout this section we use the language of track functions for charged particles, as opposed to a generic subset of particles.However, our calculations are completely generic, and can be applied to any general subset, R, of hadrons.

Calculational Technique
To verify the universality of the renormalization of the moments of the track functions, we compute it in two different ways: First we use an IRC safe observable that is directly sensitive to the track function moments, namely the EEC and projected EECs.When computed on tracks, this observable is no longer IRC safe, and the infrared poles directly determine the RG evolution of the track function moments.Second, we compute the moments of the track function by computing a jet function on tracks.This approach is computationally much simpler since it only requires the integration of splitting functions instead of complete matrix elements, but it assumes collinear factorization, and hence the universality of the track functions.The agreement between these two approaches provides a strong check both on our calculations and on the universality of the track functions.The universality of the first three moments of the track functions was tested at NLO in this same manner in [5].Here we extend this to the sixth moment.In the following two subsections we detail these two approaches.

Using Projected Energy Correlators
We begin by computing the RG for the track functions from the structure of infrared poles in energy-energy correlators, which was briefly described in [5] for the case of the two-point correlator.Here we describe it in some detail, as well as its extension to projected energy correlators, which is necessary to extract the RG of higher moments of the track functions.
The standard two-point energy correlator [17,43,44] is defined as This can be extended to a projected N -point energy correlator [6], which is sensitive to higher point correlations, but is only differential in the longest side, z L .It is defined as where X m denotes a m-particle final state and The projected correlators are IRC safe observables.However, when computed on tracks, they have collinear divergences.These collinear divergences must be absorbed by the track functions.Therefore by computing these collinear divergences, we can obtain the RG of the track functions.To simplify the notation, we combine all the products of track functions of a fixed total weight n (see (3.10)) into a vector ⃗ T n (e.g. for n = 2, ⃗ T 2 = {T g (2), T q (2), T q (1)T q (1), T g (1)T q (1), T g (1)T g (1)}).For notational simplicity, throughout this section we consider the case of a single flavor of quarks, and make the assumption T q = T q.However, we have performed the complete calculation without this assumption.Writing the renormalization group evolution of ⃗ T n as then n R (1) where a s = α s (µ)/(4π).
In terms of the tree-level track functions T (0) , we can write the two-point track EEC as dΣ dz tr = a,b∈{q j ,q j ,g} dΣ ab dz + c∈{q j ,q j ,g} The perturbatively calculable components entering this formula are given by Here f i 1 , f i 2 , f i denote the flavors of the final-state partons with the four-momenta p µ i 1 , p µ i 2 , p µ i , δ a,i 1 , δ b,i 2 and δ c,i are Kronecker deltas in flavor space, dΦ m denotes m-body phase space and M m is the corresponding matrix element.
Using that in dimensional regularization the loop corrections to the track function are scaleless, T (0) = T bare , we can employ (4.4) to rewrite (4.9) in terms of the renormalized track functions, (1) The UV poles of the track function renormalization must cancel against the IR poles in ⃗ Σ to yield a finite result, allows us to extract the RG evolution of the first and second moments of the track function.
To have access to the higher moments of the track functions, we must consider the higher point projected correlators.These proceed in a similar manner.Focusing on the three-point projected correlators, we have dΣ dz L tr = a,b,c∈{q j ,q j ,g} + a,b∈{q j ,q j ,g} + c∈{q j ,q j ,g} The perturbatively calculable components entering this formula are These have the same structure as for the two-point correlator, with the only difference being the higher energy weights.They can therefore be computed using the same techniques.The integrals Σ abc are more complicated, but fortunately the shift symmetry can be used to reconstruct the full answer from just Σ ab 2 and Σ c 3 (at least to the order at which we are currently working).More generally, for the evolution of the higher moments of the track functions, we consider the integrals and then use the shift symmetry to reconstruct the full result.These integrals can be computed using the same approach as was used to compute the standard energy correlator in ref. [19], and subsequently in refs.[20,45].This approach is an extension of the reverse unitarity method [46], which expresses delta functions from phase space constraints in terms of propagators allowing more standard loop integration techniques to be used.Using the Cutkosky rules [46,47], we express the on-shell delta functions as the cut propagators and the measurement function as , where we set the center-of-mass energy Q = (1, 0, 0, 0) for simplicity (the dependence on Q can be restored by dimensional analysis).The phase-space integrals can then be reduced to master integrals (MIs) using techniques from the study of multi-loop integrals.In particular, integration by parts and Lorentz invariance identities were generated with LiteRed [48,49] and the reduction to master integrals was performed using FIRE6 [50].The MIs are the same as that for the standard EEC and can be evaluated by the method of differential equations (DEs).The canonical forms of the DE systems are obtained by CANONICA [51].The solutions of the DEs are written in terms of harmonic polylogarithms, which can then be simplified to classical polylogarithms using HPL [52].
The calculation of Σ c p is equivalent to the calculation of cut bubble integrals, and the master integrals can be found in refs.[53,54].

Using Splitting Functions
While the calculation of the track function RG from the energy correlators provides a robust check on the universality of the track functions, it becomes computationally expensive at higher moments.Indeed, the main advantage of that approach, is that one also gets the full EEC distribution on tracks, which is itself a physically interesting observable.However, if one just wants the renormalization of the track functions, which is purely collinear in nature, it is easier to directly take advantage of collinear factorization, and obtain the RG from the splitting functions.
Here we give a general description of this approach, with more details for the case of pure Yang-Mills given in app A. Although we focus in this paper on deriving moments, this approach has the added advantage that it can be generalized to allow a derivation of the full RG of the track functions in x-space.
To obtain a non-scaleless integral in the collinear limit, one must consider the measurement of some additional observable.We consider the measurement of the jet mass of all particles and the energy fraction on charged particles, encoded in the jet function J a (s, x).The measurement of the jet mass renders the integrals non-scaleless, but importantly, the renormalization of J a (s, x) is identical to the standard J a (s) (see e.g.[55]).After performing this renormalization, as well as the standard renormalization of the strong coupling constant, the remaining poles determine the renormalization of the track functions.Unlike the pure gluon case considered in app.A, where all terms in the NLO evolution can be related to those involving three track functions, in the multi-flavor case, one must also consider terms involving two track functions.Therefore one must properly incorporate both the 1 → 3 triple collinear splitting functions [56,57], as well as the NLO corrections to the 1 → 2 splitting functions [58][59][60].
We will now provide a bit more detail for each of these steps, starting with the calculation of the jet function J a (s, x): Here Φ c N is the N -particle collinear phase space with total invariant mass s ′ and σ c a→{a f } is the squared collinear matrix element for a → a 1 a 2 • • • a N .At LO, J (0) bare,f (s, x) = δ(s)T (0) f (x).The NLO calculation of the jet function gives rise to the LO RG evolution of the track functions.To derive the NLO RG for the track functions, we must consider the NNLO calculation of J bare (s, x).
At NNLO, we have both the NLO corrections to the two-particle final state (real-virtual corrections) and the three-particle final state (real-real corrections).Explicitly, J a,bare (s, x) where σ c a→bc and σ c a→bcd are the NLO 1 → 2 splitting and LO 1 → 3 splitting functions respectively.
Taking moments of this equation and using the sum rule for the track functions, one finds that J a (s, n) is expressed in terms of integrals of the 1 → 2 and 1 → 3 splitting functions weighted by a polynomial of weight n, as is done explicitly in app.A for the pure gluon case.These integrals can be performed explicitly using the approach of [61] (many integrals relevant for the quark case can be found in [55]).
For each value of n, the renormalization of J a (s, n) in the variable s is the same as for J a (s).Renormalizing the coupling using and expanding the bare jet function and the renormalization factor in terms of the renormalized coupling, Ja , the two loop renormalization for the jet function is then a,bare + Z (0) a,bare . (4.17) The explicit form of the renormalization factors can be found in [55] (for a = q) and [62] (for a = g) up to order a 2 s .After performing this renormalization in s, the RG for the track functions can be directly read off, as for the EEC based calculation in sec.4.1.1.Explicitly, rewriting the treelevel track functions in (4.14) in terms of the renormalized track functions, using ⃗ T (0) n = ⃗ T n,bare and (4.4), the UV poles from the renormalization in (4.4) should cancel against the IR poles from the direct integration in (4.14).This should be compared with the approach in app.A, which starts from the matching of the jet function onto renormalized track functions, where the matching coefficient is finite and the IR poles are contained in the track functions.Here, instead by expressing T (0) in terms of renormalized track functions, one automatically gets something of the form of a matching relation and the resulting coefficient is therefore the finite matching coefficient.Compared to the full EEC calculation, the integrals over the splitting functions are much easier (and mostly known).However, the fact that identical results are obtained from both approaches provides a strong check on our results.

Results
In this section we present results for the first six moments of the track functions.The results for the first three moments were presented in [5] and those for the fourth through sixth moments are new.These results are provided in electronic format accompanying this paper.We write the evolution equations for the central moments, whose definition can be found in (3.20), in terms of a perturbative expansion At a given order in perturbation theory there are constraints to which combinations of track functions can appear in the evolution equations.These constraints arise from the fact that in the evolution equation of T a , a term involving the combination T b T c originates from a a → bcX splitting contribution.The constraints from the possible splittings at a given order in perturbation theory results in linear dependencies between different terms in the evolution of central moments.

Numerical Studies of Track Function Evolution
In this section we numerically study the structure of the evolution equations for the track function moments.The goal of this section is two-fold.First, we show that ∆ is sufficiently small in QCD, that corrections to DGLAP for the first three moments are effectively suppressed by (at least) an order in the perturbative expansion, allowing us to extend their RG evolution to NNLO.Second, we show that for the fourth moment and beyond, non-linearities in the evolution give rise to genuinly new behaviour beyond DGLAP.

The Size of ∆ in QCD and Extension to NNLO
We begin by studying the numerical impact of ∆ for the first three-moments.The evolution of the first three central moments is constrained by shift symmetry to be of the form where the evolution of ∆ is fixed by DGLAP to all orders.For the second and third moment the evolution can be split into two parts: a linear term fixed by DGLAP and corrections proportional to powers of ∆.Recall that ∆ = T q (1)−T g (1), or more generally in the multi-flavor case is give by differences between the first moments of the track functions of different flavors.Since QCD final states at high energies are dominated by large numbers of nearly massless pions, the average values of the track functions are largely fixed by isospin, and hence satisfy T g (1) ≃ T q (1) ≃ 2/3, and ∆ ≃ 0. Small corrections to this pictures give rise to ∆ ≪ 1 in real world QCD.This suppression of ∆, combined with the shift symmetry is particularly convenient, since it effectively suppresses the corrections to DGLAP by (at least) an order in the perturbative expansion.Indeed, we will see that this allows us to include the NNLO corrections to the DGLAP evolution while keeping the terms involving ∆ at NLO.In our numerical studies we use the following initial conditions [8], T g (1) = 0.624 , T g (2) = 0.417 , T g (3) = 0.293 , T q (1) = 0.611 , T q (2) = 0.425 , T q (3) = 0.319 , ( at µ = 10 GeV, and α s (M Z ) = 0.116 with n f = 5.
To demonstrate that the effect of ∆ on the evolution is much smaller than that of DGLAP, we study the following ratio In this ratio we compare the effect of including ∆ with the effect of including the NNLO corrections to the DGLAP evolution.The notation σ i | (N)NLO,∆=0 means setting the ∆ terms in the (N)NLO evolution to zero, but not in the lower order terms of the evolution.We note that this ratio is scale dependent, and furthermore depends strongly on the value of ∆.Since this ratio is meant to illustrate the approximate size, we have for simplicity kept the initial conditions the same for all scenarios, using the values in eq.(5.2). Figure 2 shows this ratio for a range of values of µ, which is much smaller than 1 for the second moment, as it only involves ∆ 2 terms.For the third moment, which involves terms linear in ∆, the ratio is of order 1, indicating that the ∆ terms at NLO are of the same size as the NNLO correction to the DGLAP evolution.The (unknown) ∆ terms at NNLO are of course much smaller.We further investigate the various contributions to the third moment in figure 3.Here we show the size of the NLO evolution, the ∆ term in the LO and NLO evolution, and the NNLO evolution (without ∆ term) by taking The ratio defined in eq. ( 5.3) for the quark (darker) and gluon (lighter) second (blue dashed) and third (orange dotted) central moments as a function of the renormalization scale µ.Note that the ratio for the second moment has been amplified by a factor 100 such that it is visible in this plot.The effect of ∆ on the evolution of the second central moment is much smaller than for the higher moments because ∆ appears only squared in the evolution for σ(2), while for the other moments terms linear in ∆ are also allowed.
(a) (b) Figure 3: The difference in renormalization group evolution for (a) σ g (3) and for (b) σ q (3) for the initial conditions in (5.2).Shown are the effect of the ∆ terms at LO (blue dotted), NLO (orange solid), the effect of the NLO evolution (red dot-dashed) and the NNLO evolution (green dashed).Note that two curves are multiplied by 10 for better visibility.
appropriate differences, demonstrating that the ∆ terms are effectively suppressed by one order in the perturbative expansion.The ∆ terms at NNLO can therefore safely be neglected.This allows us to immediately extend the evolution of the first three central moments of the track function to NNLO using known results for the timelike spin-n anomalous dimensions [39].This simplification is quite convenient, as it allows us to immediately consider NNLO evolution for up to three-point correlators.For convenience, we provide the DGLAP anomalous dimensions for the first three moments up to NNLO in Appendix B.

Non-Linearities in the Fourth and Fifth Moments
Although the evolution of the first three moments are DGLAP up to correction in ∆, this is not the case for higher moments.This is because the evolution of higher moments can contain non-linear terms that are not proportional to ∆ and are therefore not suppressed, even in a pure gluon theory.For example, the evolution of the fourth and fifth central moment is constrained by shift symmetry to be of the form While the terms involving ∆ are suppressed, the terms involving products of σ(2) and σ(3) are not.These non-linear terms are not constrained by DGLAP and require additional calculational techniques.Therefore extending the evolution of higher track function moments to NNLO is beyond the scope of this paper.
Let us continue to study the effects of the non-linear terms in the evolution equations.For simplicity we consider the evolution of the fourth and fifth cumulant in pure Yang-Mills theory, where the evolution of these moments simplifies to (5.5) These simplified expressions allow us to study the non-linearity of the evolution by means of a two-dimensional RG flow plot, shown in figure 4.This figure shows the RG flow for the fourth and fifth cumulant in the κ(4) − κ 2 (2) and κ(5) − κ(3)κ(2) planes respectively.From these plots it is clear that there is a single fixed-point in the evolution at the origin, corresponding to the trivial fixed point where all cumulants vanish.In addition to this fixed point, the flow lines are attracted to a common valley before flowing to the fixed point.Note that the range of the axes on these plots are somewhat arbitrary, as the figure is invariant under a simultaneous rescaling of both axes.While it is clear that the trivial fixed point is an attractive fixed point, these plots give interesting insight into the behavior of the track function.For example, we can consider a Gaussian track function for which all higher cumulants vanish.In this case, the track function  will first generate a non-zero value of κ(4) through the non-linear mixing, after which the DGLAP anomalous dimensions drive it back to zero.In this case, which is a good approximation to real world QCD, the mixing anomalous dimensions dominate the behavior of the track function evolution.Since physically the distribution must eventually collapse to a delta function under RG evolution, this suggests that there should be a positivity bound on γ κ 2 κ 2 .This provides further evidence that it may have a direct interpretation as an anomalous dimension of some generalized lightray operator, and it would be interesting to understand this better.
The RG flow of the fifth cumulant, κ(5), is interesting in that it illustrates the structure of odd moments.The RG of the track functions preserves symmetry/anti-symmetry properties under RG flow.This is manifest in the κ(5) → −κ(5), κ(3) → −κ(3) symmetry of the RG flow in the figure.For higher moments, additional non-linear terms in the evolution appear and a visualization of the RG flow can only be realized in higher-dimensional RG flow plots.
Due to the dominance of mixing terms beyond the third moment, we are not immediately able to extend our calculation to NNLO.While the complete calculation of the NNLO evolution of higher moments is beyond the scope of this paper, we briefly comment on what would be required to do so.The constraints from shift symmetry hold to all orders in perturbation theory.Focusing on pure Yang-Mills theory for simplicity, one can show that to all orders in perturbation theory the fourth moment takes the form Here we see that only one anomalous dimension, γ 1→4 , beyond the standard DGLAP anomalous dimension, appears.Interestingly, this particular contribution does not involve any soft singularities, since it has one energy weighting on each parton.Its calculation is therefore much simpler than calculations of the NNLO DGLAP kernels.It could be computed, for example, using the known 1 → 4 splitting functions [63,64].

Conclusions
Track functions characterize the fluctuations in the fragmentation process of quarks and gluons into charged hadrons (or some other subset of hadrons), and its moments are essential for the description of track-based measurements of higher-point correlation functions in jet substructure.Although they are fundamentally non-perturbative objects, the track function evolution is perturbative and exhibits interesting renormalization group structure involving mixings between different moments.
In this paper we have derived the all-orders structure of the RG for the moments of track functions, using the action of energy conservation as a shift symmetry.This highlights the remarkably constrained structure of the evolution, implying that the RG can be expressed in terms of cumulants (or equivalently, central moments), and differences of first moments.
We performed an explicit calculation of the first six moments of the quark and gluon track functions in QCD.At the fourth moment and beyond one finds interesting RG flows describing the mixing with products of cumulants, for example between κ(4) and (κ(2)) 2 .We studied the structure of these RG flows, finding that these mixing terms dominate the evolution.These higher cumulants of the track functions therefore probe evolution in the fragmentation process that goes beyond the standard DGLAP evolution, and it would be interesting to better understand the structure of these mixing terms in terms of anomalous dimensions of the underlying field theory, and study them experimentally.
Finally, we showed that for the first three moments cumulants of the track function, shift symmetry constrains any evolution beyond DGLAP to be proportional to ∆.For track-based measurements in QCD, ∆ ≪ 1, making the corrections proportional to ∆ suppressed by an effective order in the perturbative expansion.This allows us to extend the evolution to NNLO, enabling up to three-point correlators to be studied on tracks at this order.We also outlined the missing ingredients for a similar extension to NNLO beyond the third moment, where genuinely new ingredients are required.
Although we have primarily focused in this paper on the experimental utility of track functions, we believe that better understanding the evolution of the moments of the track functions could be of more formal theoretical interest.The DGLAP anomalous dimensions have a deep connection to the twist-2 operators of the theory, which has recently recieved renewed attention in the study of lightray operators in CFTs [9].Track functions are another class of intrinsically Lorentzian observables, that probe features of the theory beyond the leading twist trajectory.It would be interesting if they could be put on a similarly firm theoretical footing, and if one could more precisely understand what features of the theory they are probing, and how they are related to its operator content.
Our results allow the calculation of up-to six point energy correlators on tracks, which have recently been investigated with CMS open data [65][66][67] providing a view on the hadronization transition, non-Gaussianities and quantum scaling dimensions.The three-point energy correlator has also been proposed as a new way to extract the top quark mass [68], with the potential to reduce the theoretical uncertainty, particularly from nonperturbative effects.The angular resolution offered by tracks is essential to carry out these measurements.This is also the case for the azimuthal decorrelation in vector-boson plus jet production [69], which however requires knowledge of (the evolution of) the full track function.In conclusion, we believe that our work will be of significant interest for precision studies at the LHC, and we look forward to their application in phenomenology in the near future.
where a s = α s /(4π).The Mellin moments of timelike splitting functions are Note that this is shifted by one from the definition of the moments of the track function in eq.(3.2).All the results of P ij (z) up to order-a 3 s are e.g.listed in the ancillary file, "PT.txt", of [39], and P ij (z) corresponds to PT["ij"] in that file.At LO, P (0) qq and P (0) Qq vanish while the non-vanishing moments up to the 7th moment are given by For Q ̸ = q we have γ Qq = γ Qq and up to the 7th moment we have γ For the EEC evolution to NNLL, we need the N = 3 moment LO, NLO and NNLO, which can be obtained from refs.[40][41][42]71].(Note that we include the pure singlet term in the qq element.)At NNLO, we have  Up to three-loop order in the MS scheme, the coefficients of the β function are [72,73]

Figure 2 :
Figure 2:The ratio defined in eq.(5.3) for the quark (darker) and gluon (lighter) second (blue dashed) and third (orange dotted) central moments as a function of the renormalization scale µ.Note that the ratio for the second moment has been amplified by a factor 100 such that it is visible in this plot.The effect of ∆ on the evolution of the second central moment is much smaller than for the higher moments because ∆ appears only squared in the evolution for σ(2), while for the other moments terms linear in ∆ are also allowed.

Figure 4 :
Figure 4: Renormalization group flow in pure Yang-Mills theory at fixed µ for (a) the fourth cumulant and (b) for the fifth cumulant.The arrows denote the direction of the derivatives with respect to ln µ and their color reflects their strength.The black line indicates the eigenvector of the evolution equation.