  Jets and their substructure play a central role in many analyses at the Large
Hadron Collider (LHC). To improve the precision of measurements, as well as to
enable measurement of jet substructure at increasingly small angular scales,
tracking information is often used due to its superior angular resolution and
robustness to pile-up. Calculations of track-based observables involve
non-perturbative track functions, that absorb infrared divergences in
perturbative calculations and describe the transition to charged hadrons. The
infrared divergences are directly related to the renormalization group
evolution (RGE), and can be systematically computed in perturbation theory.
Unlike the standard DGLAP evolution, the RGE of the track functions is
non-linear, encoding correlations in the fragmentation process. We compute the
next-to-leading order (NLO) evolution of the track functions, which involves in
its kernel the full $1\rightarrow3$ splitting function. We discuss in detail
how how we implement the evolution equation numerically, and illustrate the
size of the NLO corrections. We also show that our equation can be viewed as a
master equation for collinear evolution at NLO, by illustrating that by
integrating out specific terms, one can derive the evolution for any $N$-hadron
fragmentation function. Our results provide a crucial ingredient for obtaining
track-based predictions for generic measurements at the LHC, and for improving
the description of the collinear dynamics of jets.
