Generalized Fragmentation Functions for Fractal Jet Observables

TL;DR

This paper introduces a new class of fractal jet observables based on generalized fragmentation functions, which are capable of capturing complex, correlated properties of hadrons in jets, with applications in quark-gluon discrimination.

Contribution

It develops a formalism for nonlinear evolution of generalized fragmentation functions for fractal jet observables, including a numerical code and phenomenological applications.

Findings

Nonlinear evolution equations for generalized fragmentation functions are derived.

A numerical code for evolution of fractal jet observables is developed.

Fractal jet observables can effectively discriminate quark from gluon jets.

Abstract

We introduce a broad class of fractal jet observables that recursively probe the collective properties of hadrons produced in jet fragmentation. To describe these collinear-unsafe observables, we generalize the formalism of fragmentation functions, which are important objects in QCD for calculating cross sections involving identified final-state hadrons. Fragmentation functions are fundamentally nonperturbative, but have a calculable renormalization group evolution. Unlike ordinary fragmentation functions, generalized fragmentation functions exhibit nonlinear evolution, since fractal observables involve correlated subsets of hadrons within a jet. Some special cases of generalized fragmentation functions are reviewed, including jet charge and track functions. We then consider fractal jet observables that are based on hierarchical clustering trees, where the nonlinear evolution equations also exhibit tree-like structure at leading order. We develop a numeric code for performing this evolution and study its phenomenological implications. As an application, we present examples of fractal jet observables that are useful in discriminating quark jets from gluon jets.

Equations

$$x=\hat{x}(z,x_1,x_2),$$

$$z\equiv \frac{E_1}{E_1+E_2}$$

$$\mu \frac{\text{d}}{\text{d}\mu }{\mathcal{F}}_i(x,\mu )=\frac{1}{2}\sum _{j,k}\int \text{d}z\text{d}{x}_1\text{d}{x}_2\frac{{\alpha }_s(\mu )}{\pi }{P}_{i\to jk}(z){\mathcal{F}}_j(x_1,\mu ){\mathcal{F}}_k(x_2,\mu )\delta [x-\hat{x}(z,x_1,x_2)],$$

$$x=\sum _{a\in \mathrm{jet}}{w}_a{z}_a^{\kappa },z_a\equiv \frac{E_a}{E_{\mathrm{jet}}},$$

$$x=\sum _{\mathrm{nodes}}{\left(4z_L{z}_R\right)}^{\kappa /2},$$

$$\sum _h{\int }_0^1\text{d}xx{D}_i^h(x,\mu )=1.$$

$$\displaystyle D_{i}^{h}(x,\mu)=\frac{1}{x}\int\text{d}^{2}p_{h}^{\perp}\!\int\frac{\text{d}y^{+}\text{d}^{2}y_{\perp}}{2(2\pi)^{3}}e^{ip^{-}y^{+}}\sum_{X}\frac{1}{2N_{C}}\mbox{Tr}\bigg{[}\frac{\gamma^{-}}{2}\langle 0|\psi_{i}(y^{+},0,y_{\perp})|hX\rangle\langle hX|\overline{\psi}_{i}(0)|0\rangle\bigg{]},$$

$$w^{\mu }=w^{-}\frac{n^{\mu }}{2}+w^{+}\frac{{\bar{n}}^{\mu }}{2}+w_{\perp }^{\mu },w^{-}=\bar{n}\cdot w,w^{+}=n\cdot w,$$

$$\frac{1}{{\sigma }^{(0)}}\frac{\text{d}\sigma }{\text{d}x}(e^{+}{e}^{-}\to hX)=\sum _i{\int }_x^1\frac{\text{d}z}{z}{C}_i(z,s,\mu )D_i^h(x/z,\mu ),$$

$$\mu \frac{\text{d}}{\text{d}\mu }{D}_i^h(x,\mu )=\sum _j{\int }_x^1\frac{\text{d}z}{z}\frac{{\alpha }_s(\mu )}{\pi }{P}_{ji}(z)D_j^h(x/z,\mu ).$$

$$P_{ji}(z)=P_{ji}^{(0)}(z)+\frac{{\alpha }_s}{2\pi }{P}_{ji}^{(1)}(z)+\dots ,$$

$$P_{ji}^{(0)}(z)\equiv P_{i\to jk}(z),$$

$$\mu \frac{\text{d}}{\text{d}\mu }{D}_i^h(x,\mu )$$

$$\displaystyle\quad\times\Bigl{(}D_{j}^{h}\left(x_{1},\mu\right)\,\delta[x-zx_{1}]+D_{k}^{h}\left(x_{2},\mu\right)\,\delta[x-(1-z)x_{2}]\Bigr{)}.$$

$$\int \text{d}x{\mathcal{F}}_i(x,\mu )=1.$$

$${\mathcal{F}}_i(x,\mu )$$

$$\displaystyle\quad\times\mbox{Tr}\bigg{[}\frac{\gamma^{-}}{2}\langle 0|\psi_{i}(y^{+},0,y_{\perp})|SX\rangle\langle SX|\overline{\psi}_{i}(0)|0\rangle\bigg{]},$$

$${\mathcal{F}}_g(x,\mu )$$

$$\times ⟨0|G_{\lambda }^{-,a}(y^{+},0,y_{\perp })|SX⟩⟨SX|G^{-,a,\lambda }(0)|0⟩,$$

$$\frac{1}{{\sigma }^{(0)}}\frac{\text{d}\sigma }{\text{d}{z}_J\text{d}x}(e^{+}{e}^{-}\to \mathrm{jet}+X)=\sum _i\int \frac{\text{d}{y}^{\prime }}{y^{\prime }}{C}_i(z_J/y^{\prime },E_{\mathrm{cm}},\mu )$$

$$\displaystyle\quad\times\Big{\{}\delta(1-y^{\prime})\,\mathcal{F}_{i}(x,\mu)+\sum_{j}\mathcal{J}^{(1)}_{i\to j}(y^{\prime},E_{\rm jet}R,\mu)\,\mathcal{F}_{j}(x,\mu)$$

$$+\delta (1-y^{\prime })\frac{1}{2}\sum _{j,k}\int \text{d}z\text{d}{x}_1\text{d}{x}_2{\mathcal{J}}_{i\to jk}^{(1)}(z,E_{\mathrm{jet}}R,\mu ){\mathcal{F}}_j(x_1,\mu ){\mathcal{F}}_k(x_2,\mu )\delta [x-\hat{x}(z,x_1,x_2)]$$

$$\displaystyle\quad+\frac{1}{2}\sum_{j,k}\int\text{d}z\,\text{d}x_{1}\,\text{d}x_{2}\,\mathcal{J}^{(2)}_{i\rightarrow jk}(y^{\prime},z,E_{\rm jet}R,\mu)\,\mathcal{F}_{j}(x_{1},\mu)\,\mathcal{F}_{k}(x_{2},\mu)\,\delta[x-\hat{x}(z,x_{1},x_{2})]+\dots\Big{\}},$$

$$\mu \frac{\text{d}}{\text{d}\mu }{\mathcal{F}}_i(x,\mu )=\frac{1}{2}\sum _{j,k}\int \text{d}z\text{d}{x}_1\text{d}{x}_2\frac{{\alpha }_s(\mu )}{\pi }{P}_{i\to jk}(z){\mathcal{F}}_j(x_1,\mu ){\mathcal{F}}_k(x_2,\mu )\delta [x-\hat{x}(z,x_1,x_2)],$$

$$\mu \frac{\text{d}}{\text{d}\mu }\int \text{d}x{\mathcal{F}}_i(x,\mu )=\frac{1}{2}\sum _{j,k}\int \text{d}z\frac{{\alpha }_s(\mu )}{\pi }{P}_{i\to jk}(z)\int \text{d}{x}_1{\mathcal{F}}_j(x_1,\mu )\int \text{d}{x}_2{\mathcal{F}}_k(x_2,\mu )=0,$$

$${\mathcal{F}}_i^{(1)}(x)$$

$$\times \int \text{d}{x}_1\text{d}{x}_2{\mathcal{F}}_j^{(0)}(x_1,\mu ){\mathcal{F}}_k^{(0)}(x_2,\mu )\delta [x-\hat{x}(z,x_1,x_2)].$$

$$\hat{x}(z,x_1,x_2)=z{x}_1+(1-z)x_2,$$

$$d_{ij}=\min [E_i^{2p},E_j^{2p}]{\Omega }_{ij}^2,$$

$$\underset{z\to 1}{\lim }\hat{x}(z,x_1,x_2)=x_1,\underset{z\to 0}{\lim }\hat{x}(z,x_1,x_2)=x_2,$$

$$\mu \frac{\text{d}}{\text{d}\mu }{\mathcal{F}}_i=\frac{{\alpha }_s}{\pi }{P}_{i\to jk}\otimes {\mathcal{F}}_j\otimes {\mathcal{F}}_k+{\left(\frac{{\alpha }_s}{\pi }\right)}^2{P}_{i\to jk\ell }\otimes {\mathcal{F}}_j\otimes {\mathcal{F}}_k\otimes {\mathcal{F}}_{\ell }+\dots ,$$

$$\hat{x}(z,x_1,x_2)=x_1{z}^{\kappa }+x_2(1-z{)}^{\kappa },$$

$$x=\sum _{a\in \mathrm{jet}}{w}_a{z}_a^{\kappa },z_a\equiv \frac{E_a}{E_{\mathrm{jet}}},$$

$$x_{\mathrm{A}}$$

Full text

††preprint:

MIT-CTP 4887

NIKHEF 2017-002

CERN-TH-2017-064

1 Introduction

Fragmentation functions (FFs) have a long history in QCD for calculating cross sections for collinear-unsafe observables. Ordinary FFs are process-independent nonperturbative objects that describe the flow of momentum from a fragmenting quark or gluon into an identified final-state hadron Georgi:1977mg ; Mueller:1978xu ; Ellis:1978ty ; Curci:1980uw ; Altarelli:1981ax ; Collins:1981uk ; Collins:1981uw . Since the momentum of a single hadron is not collinear safe, cross sections for single-hadron observables have singularities beginning at $\mathcal{O}(\alpha_{s})$. These collinear singularities are absorbed by the FFs order by order in $\alpha_{s}$. From this singularity structure, one can derive the renormalization group (RG) evolution for FFs, leading to the well-known DGLAP equations Lipatov:1974qm ; Gribov:1972ri ; Altarelli:1977zs ; Dokshitzer:1977sg . This evolution is linear, since FFs depend only on the momentum of a single hadron in the final state.

In this paper, we present a formalism for generalized fragmentation functions (GFFs), which describe the flow of momentum from a fragmenting quark or gluon into subsets of final-state hadrons. Because GFFs depend on correlations between final-state hadrons, their evolution equations are nonlinear and therefore more complicated than in the ordinary FF case. Motivated by the structure of the DGLAP equations, we define fractal jet observables where the evolution, albeit nonlinear, takes a special recursive form that is well-suited to numerical evaluation.111This should not be confused with “extended fractal observables” recently introduced in LABEL:Davighi:2017hok, which are based on determining the fractal dimension of a jet.

Specifically, we focus on observables defined using hierarchical binary clustering trees that mimic the leading-order tree-like structure of the evolution equations. A fractal jet observable $x$ can then be defined recursively according to fig. 1 as

[TABLE]

where $x_{1}$ and $x_{2}$ are the values of the observable on the branches of a $1\to 2$ clustering tree, and $z$ is the momentum sharing between branches, defined by

[TABLE]

with $E_{i}$ the energy of branch $i$.222While it would be more accurate to call eq. (2) the “energy fraction”, we use momentum fraction since that is more common in the fragmentation function literature. With these definitions, the leading-order evolution equation of the corresponding GFF takes the simplified form

[TABLE]

where $\mathcal{F}_{i}(x,\mu)$ is the GFF for parton $i=\{u,\bar{u},d,\dots,g\}$, $P_{i\rightarrow jk}(z)$ is the $1\to 2$ QCD splitting function, and $\mu$ is the $\overline{\rm MS}$ renormalization scale. This evolution equation has the same structure as a $1\to 2$ parton shower, which is sufficiently straightforward to implement numerically. Although we mostly restrict ourselves to lowest order in perturbation theory, our framework allows for the systematic inclusion of higher-order corrections, in contrast to the semi-classical parton shower approach.

The class of fractal jet observables described by eq. (1) is surprisingly rich, allowing for many collinear-unsafe observables to be calculated with the help of GFFs. For example, eq. (3) describes the evolution of weighted energy fractions,

[TABLE]

where $w_{a}$ is a weight factor that depends on non-kinematic quantum numbers such as charge or flavor, $\kappa>0$ is an energy weighting exponent, and the sum extends over all jet constituents. These observables are defined by associative recursion relations, such that their value is independent of the choice of clustering tree. Examples of weighted energy fractions include weighted jet charge FIELD19781 , whose nonlinear evolution was first studied in LABEL:Waalewijn:2012sv; track functions which characterize the fraction of a jet’s momentum carried by charged particles Chang:2013rca ; Chang:2013iba ; and the observable $p_{T}^{D}$ used by the CMS experiment for quark/gluon discrimination Pandolfi:1480598 ; Chatrchyan:2012sn , whose nonlinear evolution was first studied in LABEL:Larkoski:2014pca. While we focus on the case of $e^{+}e^{-}$ collisions with jets of energy $E_{\rm jet}$, our formalism easily adapts to hadronic collisions with jets of transverse momentum $p_{T}^{\rm jet}$.

In addition to performing a more general analysis of weighted energy fractions, we also present examples of fractal observables with non-associative recursion relations. These quantities depend on the details of the clustering tree used to implement eq. (1), providing a complementary probe of jet fragmentation. In particular, while eq. (1) does not involve any explicit angular separation scales, the clustering tree does introduce an implicit angular dependence. Remarkably, the details of the clustering do not affect the leading-order RG evolution in eq. (3) considered in this paper, beyond the requirement that particles are appropriately clustered in the collinear limit. An example of a non-associative fractal observable is given by node-based energy products,

[TABLE]

where the observable depends on the momentum fractions carried by the left and right branches at each node in the clustering tree. We also study observables defined entirely in terms of eq. (1), with no obvious simplification. This sensitivity to the tree structure allows non-associative observables to probe parton fragmentation from a different perspective than previously-studied jet observables. As one application, we consider the discrimination between quark- and gluon-initiated jets (see e.g. Gallicchio:2011xq ; Gallicchio:2012ez ; Larkoski:2013eya ; Larkoski:2014pca ; Bhattacherjee:2015psa ; Badger:2016bpw ; FerreiradeLima:2016gcz ; Komiske:2016rsd ; Gras:2017jty for recent studies). We find that fractal observables are effective for this purpose, in some cases yielding improved quark/gluon separation power compared to weighted energy fractions.

For clustering trees obtained from the Cambridge/Aachen (C/A) algorithm Dokshitzer:1997in ; Wobisch:1998wt , the depth in the tree is directly related to the angular separation scale between subjets. This opens up the possibility of modifying the recursion relation $\hat{x}$ in eq. (3) to be a function of angular scale. For example, starting from a jet of radius $R$, one can introduce a subjet radius parameter $R_{\rm sub}\ll R$ such that evolution equation takes a different form below and above $R_{\rm sub}$. A particularly simple case is if the weighted energy fraction with $\kappa=1$ is measured on the branches below $R_{\rm sub}$, since this effectively amounts to defining fractal observables in terms of subjets of radius $R_{\rm sub}$. In this case, the initial conditions for the GFF leading-order evolution is simply given by $\mathcal{F}_{i}(x,\mu_{\rm sub})=\delta(1-x)$ at the initial scale $\mu_{\rm sub}=E_{\rm jet}R_{\rm sub}\gg\Lambda_{\rm QCD}$, such that no nonperturbative input is needed. By evolving the GFFs to $\mu=E_{\rm jet}R$, we achieve the resummation of leading logarithms of $R_{\rm sub}/R$. Related evolution techniques have been used to resum logarithms of the jet radius $R$ in inclusive jet cross sections Dasgupta:2014yra ; Kang:2016mcy ; Dai:2016hzf .

The formalism of GFFs is reminiscent of other multi-hadron FFs in the literature. This includes dihadron fragmentation functions which describe the momentum fraction carried by pairs of final-state hadrons Sukhatme:1980vs ; Vendramin:1981te , and fracture functions which correlate the properties of one initial-state and one final-state hadron Trentadue:1993ka ; Graudenz:1994dq . In all of these cases, the RG evolution equations are nonlinear. The key difference here is that fractal jet observables are not based on a fixed number of hadrons, but rather allow for arbitrary hadron multiplicities. Depending on the observable, this may require that all hadrons can be consistently labeled by non-kinematic quantum numbers (e.g. charge). As discussed in LABEL:Waalewijn:2012sv for the case of weighted jet charge, the $n$-th moment of GFFs can sometimes be related to moments of $n$-hadron FFs. At the level of the full distribution, though, GFFs are distinct from multi-hadron FFs, and thereby probe complementary aspects of jet fragmentation.

The rest of this paper is organized as follows. In sec. 2, we review the theoretical underpinnings of ordinary parton fragmentation and explain how to extend the formalism to generalized fragmentation and fractal observables. We then construct generic fractal jet observables using clustering trees in sec. 3. In sec. 4, we treat the case of weighted energy fractions, exploring their RG evolution for a range of parameters. We introduce two new sets of non-associative fractal observables in sec. 5—node products and full-tree observables—and motivate their application in quark/gluon discrimination in sec. 6. We briefly explain how our formalism also applies to fractal observables based on subjets rather than hadrons in sec. 7. We conclude in sec. 8, leaving calculational details and a description of the numerical RG implementation to the appendices.

2 Formalism

To motivate the definition of fractal jet observables, it is instructive to first review the formalism of standard fragmentation and then generalize it to arbitrary collinear-unsafe observables. We give a general definition of fractal jet observables at the end of this section, which serves as a preamble to the explicit constructions in sec. 3.

2.1 Review of Standard Fragmentation

Ordinary FFs, denoted by $D_{i}^{h}(x,\mu)$, are nonperturbative objects that describe the number density of hadrons of type $h$ carrying momentum fraction $x$ among the particles resulting from the fragmentation of a parton of type $i$. They are the final-state counterpart to parton distribution functions (PDFs). For any parton flavor $i$, they satisfy the momentum conservation sum rule

[TABLE]

At leading order, the FFs are independent of the factorization scheme (see e.g. Furmanski:1981cw ).

The field-theoretic definition of the bare unpolarized quark FF is given by Collins:1981uk ; Collins:1981uw

[TABLE]

where we are working in a frame with quark transverse momentum $\vec{p}_{\perp}=0$ and using the gauge choice $A^{-}=0$. The jet-like state $|hX\rangle$ contains an identified hadron $h$ of momentum $p_{h}$ with $p_{h}^{-}\equiv xp^{-}$, and $X$ refers to all other hadrons in that state. The factor $1/(2N_{C})$, where $N_{C}=3$ is the number of colors, accounts for averaging over the color and spin of the quark field $\psi$ of flavor $i$. Here and in the rest of the paper, we adopt the following convention for decomposing a four-vector $w^{\mu}$ in light-cone coordinates:

[TABLE]

where $n^{\mu}$ is a light-like vector along the direction of the energetic parton, and $\bar{n}$ is defined such that $n^{2}=\bar{n}^{2}=0$ and $n\cdot\bar{n}=2$. Thus at leading order, $p^{-}=2E_{\rm jet}$. Gauge invariance requires adding eikonal Wilson lines in eq. (7) (see e.g. Collins:1989gx ), which we suppress here for notational convenience. An analogous definition applies for the gluon FF.

In the context of $e^{+}e^{-}$ annihilation, FFs are crucial ingredients in the factorization formula for the semi-inclusive cross section at leading power in $\Lambda_{\rm QCD}/\sqrt{s}$,

[TABLE]

where $x=2E_{h}/\sqrt{s}\leq 1$ is the hadron energy fraction, $\sigma^{(0)}$ is the tree-level cross section and $X$ represents all other final state particles in the process.333In the literature (see e.g. Olive:2016xmw ), the cross section $1/\sigma^{(0)}\,\text{d}\sigma/\text{d}x(e^{+}e^{-}\rightarrow hX)=F^{h}(x,\mu)$ is sometimes referred to as the total FF, in which case $D_{i}^{h}(x,\mu)$ is called the parton FF. The coefficients $C_{i}(z,s,\mu)$ are process-dependent perturbative functions that encode the physics of the hard subprocess. The FFs $D_{i}^{h}(x,\mu)$ are universal, process-independent functions, which appear (with appropriate PDF convolutions) in related channels such as $ep\rightarrow hX$ or $pp\rightarrow hX$. Since the coefficients $C_{i}$ contain logarithms of $s/\mu^{2}$, in order to avoid terms that could spoil perturbative convergence in eq. (9), the renormalization scale $\mu$ should be chosen close to $\sqrt{s}$.

While computing the FFs themselves requires nonperturbative information about the hadronic matrix elements in eq. (7), their scale dependence is perturbatively calculable. This allows us to, for example, take FFs extracted from fits to experimental data at one scale and evolve them to another perturbative scale. The RG evolution of FFs is described by the DGLAP equations Lipatov:1974qm ; Gribov:1972ri ; Altarelli:1977zs ; Dokshitzer:1977sg ,

[TABLE]

Here, the splitting kernels $P_{ji}(z)$ can be calculated in perturbation theory,

[TABLE]

and are at lowest order the same as the splitting kernels for PDF evolution. The next-order splitting function $P^{(1)}_{ji}$ arises from $1\to 3$ splittings as well as loop corrections to $1\to 2$ splittings.

In order to motivate the transition to generalized fragmentation, it is convenient to rewrite the lowest-order splitting function explicitly as a $1\to 2$ process:

[TABLE]

where the parton $j$ carries momentum fraction $z$, e.g. $P_{g\rightarrow gg}(z)$ or $P_{q\rightarrow qg}(z)=P_{q\rightarrow gq}(1-z)$. With this notation, we can rewrite the leading-order DGLAP equation in a suggestive form444Because the splitting functions are divergent as $z\rightarrow 1$ and as $z\rightarrow 0$, plus-function regulators are required at both endpoints when integrating over the entire range $0\leq z\leq 1$.

[TABLE]

Though we have written eq. (2.1) as an integral over both $x_{1}$ and $x_{2}$, corresponding to the two final state branches from the $i\to jk$ splitting, the FFs only require information about one single final-state hadron in each term, so the evolution simplifies to the linear form in eq. (10). This will no longer be the case with generalized fragmentation, which depends on correlations between the final-state hadrons.

2.2 Introducing Generalized Fragmentation

We now extend the FF formalism to handle the distribution of quantities $x$ carried by a subset $S$ of collinear particles, where $x$ can be more general than the simple momentum fraction and $S$ is defined by non-kinematic quantum numbers. For example, we will consider observables defined on all particles within a jet, but also on charged particles only. For a given observable $x$, there is a GFF for each parton species $i$, which we denote by $\mathcal{F}_{i}(x,\mu)$. At lowest order in $\alpha_{s}$, the GFF is the probability density for the particles in $S$ to yield a value of the observable $x$ from jets initiated by a parton of type $i$. The GFF automatically includes information about hadronization fluctuations. Being a probability density, the GFFs are normalized to unity for each parton type,

[TABLE]

For any collinear-unsafe (but soft-safe) observable $x$, we can give an operator definition for GFFs analogous to that for fragmentation functions. A (bare) quark GFF for the gauge choice $A^{-}=0$ is defined as

[TABLE]

to be compared with eq. (7). Here, $|SX\rangle$ is the asymptotic final state divided into the measured subset $S$ and unmeasured subset $X$, and $\tilde{x}(p^{-},S)$ is the functional form of the quantity being observed, which can depend on the overall jet momentum and any information from $S$. We stress that, in contrast to the standard FFs, a GFF involves a sum over polarizations and a phase-space integration over all detected particles in $S$; if the measured set $S$ consists of a single hadron, then eq. (2.2) reduces to eq. (7) for a quark FF. The definition for gluon-initiated jets is

[TABLE]

where $G^{-,a}_{\lambda}=\overline{n}^{\mu}G_{\mu\lambda}^{a}$ is the gluon field strength tensor for generator $T^{a}$, the factor of $1/(d-2)$ comes from averaging over the gluon polarizations in $d$ space-time dimensions, and the factor of $1/(N_{C}^{2}-1)$ comes from averaging over the color of the gluon.

The definitions in eqs. (2.2) and (2.2) extend the ones introduced in LABEL:Chang:2013rca for track functions. In the track function case, $x$ is the momentum fraction carried by the charged particles in the final states, irrespective of their individual properties or multiplicities. As mentioned in the introduction, GFFs are reminiscent of multi-hadron FFs Sukhatme:1980vs ; Vendramin:1981te , with the key difference that multi-hadron FFs describe a fixed number of identified final-state hadrons (i.e. two in the case of dihadron FFs), whereas GFFs allow for a variable number of final-state hadrons in the subset $S$.

With these GFFs in hand, we can calculate the cross section differential in the fractal observable $x$ for an inclusive jet sample with radius parameter $R\ll 1$. Letting $z_{J}$ be the fraction of the center-of-mass energy carried by the measured jet ($z_{J}\equiv 2E_{\rm jet}/E_{\rm cm}$), we have

[TABLE]

where the ellipsis includes further terms at next-to-next-to leading order and $\sigma^{(0)}$ denotes the tree-level cross section. There is a similar version of eq. (17) for $pp$ and $ep$ collisions with the inclusion of PDFs, where the jet rapidity would appear in the $C_{i}$ coefficients. As in eq. (9), the effects of the hard interaction producing a parton $i$ are encoded in the coefficients $C_{i}$, which can be expanded perturbatively and depend on $z_{J}$ and $E_{\rm cm}$. At leading order, the jet only consists of parton $i$, thus $C_{i}^{(0)}(z_{J})=\delta(1-z_{J})$ and the dependence on the fractal observable $x$ arising from parton production and hadronization is described simply by $\mathcal{F}_{i}$. For most of the paper, we restrict ourselves to leading order, though we stress that eq. (17) provides the tools to interface our GFF formalism with fixed-order calculations and to extract GFFs beyond leading order.

At next-to-leading order in eq. (17), the parton $i$ can undergo a perturbative splitting into partons $j$ and $k$. If only $j$ is inside the jet then $z_{J}<1$, as described by the perturbative coefficient $\mathcal{J}_{i\to j}^{(1)}$ that can be derived from ref. Kang:2016mcy , and the $x$-dependence is described by $\mathcal{F}_{j}$. If both partons belong to the jet then again $z_{J}=1$, but the observable $x$ now follows from combining the values $x_{1}$ and $x_{2}$ of the GFFs for partons $j$ and $k$ with the momentum fraction $z$ of the perturbative splitting described by the $\mathcal{J}_{i\to jk}^{(1)}$ from ref. Waalewijn:2012sv . At next-to-next-to-leading order, there are even more contributions, including one with three partons in the jet involving $\mathcal{J}_{i\to jk\ell}^{(2)}$. In eq. (17), we displayed only the term with two partons belonging to the jet, since it is the first term that directly correlates $z_{J}$ and $z$. The natural scale of the coefficients $\mathcal{J}_{i\to j},\mathcal{J}_{i\to jk},\dots,$ is the typical jet invariant mass $E_{\rm jet}R$, so we conclude that the GFFs should be evaluated at $\mu\simeq E_{\rm jet}R$ to minimize the effect of higher-order corrections. If $R\gtrsim 1$, then $C_{i}$ and $\mathcal{J}$ can be combined, and the natural scale to evaluate the GFF would be $\mu\simeq E_{\rm jet}$.

It is important to note that eq. (17) really combines two different formalisms. The first is the formalism for GFFs discussed initially in refs. Waalewijn:2012sv ; Chang:2013rca for track-based observables and further developed here. The second is the formalism for fragmentation in inclusive jet production of refs. Dai:2016hzf ; Kang:2016ehg , which builds upon work on fragmentation in exclusive jet samples refs. Procura:2009vm ; Jain:2011xz ; Liu:2010ng ; Jain:2011iu . Both of these formalisms are needed to perform higher-order jet calculations, though at leading order, the GFF formalism alone suffices. For the interested reader, we provide all details of the matching for $e^{+}e^{-}\to{\rm jet}+X$ at next-to-leading order in app. A. As in refs. Waalewijn:2012sv ; Chang:2013rca , we expect that the absorption of collinear divergences by GFFs can be carried out order-by-order in $\alpha_{s}$ due to the universality of the collinear limits in QCD.

2.3 Introducing Fractal Observables

The above generalized fragmentation formalism works for any collinear-unsafe (but soft-safe) observable. The RG evolution for a generic $\mathcal{F}_{i}(x,\mu)$, however, can be very complicated. In order to deal with numerically tractable evolution equations, we focus on observables whose RG evolution simplifies to a nonlinear version of eq. (2.1). Specifically, we want to find the most general form of the function $\tilde{x}(p^{-},S)$ in eqs. (2.2) and (2.2) such that the RG evolution of $\mathcal{F}_{i}(x,\mu)$ depends only on itself and other GFFs for the same observable, and does not mix with other functions. An example of an observable that involves GFF mixing is given in app. B, where the evolution equation is considerably more complicated than considered below.

We define fractal observables as those whose GFFs obey the (leading-order) RG equation in eq. (3), repeated here for convenience:

[TABLE]

where $\hat{x}(z,x_{1},x_{2})$ is a function related to $\tilde{x}(p^{-},S)$, which now depends on the momentum $p$ only through the momentum sharing $z$. As advertised, the evolution of $\mathcal{F}_{i}(x,\mu)$ depends only on GFFs for the same observable $x$, and no other nonperturbative functions. We leave a detailed discussion of higher-order evolution to future work, and focus primarily on the leading-order evolution here. As a consistency check, the $\delta$ function in eq. (18) ensures that the RG evolution automatically preserves the GFF normalization,

[TABLE]

where we used the fact that $\sum_{j,k}\int\text{d}z\,P_{i\rightarrow jk}(z)=0$.

As a simple example of a fractal observable, consider the momentum fraction $x$ carried by a subset $S$ of hadrons of a common type. This case has already been studied in the context of track functions Chang:2013rca ; Chang:2013iba , where $S$ corresponded to charged particles. Treating the states $|SX\rangle$ in eqs. (2.2) and (2.2) partonically, the next-to-leading-order bare GFF in dimensional regularization with $d=4-2\epsilon$ satisfies

[TABLE]

Here, the function $\hat{x}(z,x_{1},x_{2})$ is the form of $\tilde{x}(p^{-},S)$ written in terms of two subjets,

[TABLE]

where $x_{1}$ and $x_{2}$ are the momentum fractions carried by particles belonging to subjets 1 and 2 within $S$, and $z$ is the momentum fraction carried by subjet 1, as defined in eq. (2). Renormalizing the UV divergences in eq. (2.3) in the $\overline{\rm MS}$ scheme leads directly to the RG equation in eq. (18). Thus, the momentum fraction $x$ carried by the final-state subset $S$ is indeed a fractal observable.

In the above analysis, we implicitly assumed massless partons, since otherwise the parton mass $m$ would regulate the $1/\epsilon_{\rm IR}$ divergence. As long as $m\ll E_{\rm jet}R$, it is consistent to take the $m\to 0$ limit, which resums the large logarithms of $E_{\rm jet}R/m$ in the cross section for the fractal observable. At the scale $\mu=m$, one has to match the GFF evolution onto the appropriate heavy-quark description.

3 Fractal Observables via Clustering Trees

We now present a straightforward way to build a broad class of fractal observables that have the desired RG evolution in eq. (18). The idea is to use recursive clustering trees that mimic the structure of the leading-order RG evolution equations. Our construction is based on the following three ingredients, as shown in fig. 2:

Weights $w_{a}$ for each final-state hadron; 2. 2.

An IRC-safe binary clustering tree; 3. 3.

The recursion relation $\hat{x}(z,x_{1},x_{2})$.

By implementing the function $\hat{x}$ directly on recursive clustering trees, the resulting observable is guaranteed to have fractal structure.

3.1 Construction

For this discussion, we start with a collection of hadrons from an identified jet, found using a suitable jet algorithm, e.g. anti-$k_{t}$ Cacciari:2008gp in the studies below. As the initial boundary condition for the observable, each final-state hadron within the jet is assigned a weight $w_{a}$ (possibly zero) based on some non-kinematic quantum number associated with that hadron. This weight controls how much each type of hadron contributes to the value of the jet observable. For example, to construct an observable that only depends on the charged particles in the jet, all charged particles would be given weight 1 and all neutral particles weight 0. It is crucial that $w_{a}$ is independent of the energy and direction of the hadron, otherwise the NLO GFF would not take the form in eq. (2.3).

These final-state hadrons are then used as inputs to an IRC-safe binary clustering tree, which is in general different from any clustering algorithm used to determine the identified jet. For our studies, we use the generalized-$k_{t}$ family of jet clustering algorithms Cacciari:2008gp , which are designed to follow the leading-order structure of the parton shower. In the context of $e^{+}e^{-}$ collisions, these algorithms have the pairwise clustering metric

[TABLE]

where the exponent $p$ parametrizes the tree-dependence of the observable, with $p=\{-1,0,1\}$ corresponding to the \{\text{anti-k_{t}~{}\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{Cacciari:2008gp}{\@@citephrase{(}}{\@@citephrase{)}}}},\text{C/A~{}\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{Dokshitzer:1997in,Wobisch:1998wt}{\@@citephrase{(}}{\@@citephrase{)}}}},\text{k_{t}~{}\cite[cite]{\@@bibref{Authors Phrase1YearPhrase2}{Catani:1993hr,Ellis:1993tq}{\@@citephrase{(}}{\@@citephrase{)}}}}\} clustering algorithms, and $\Omega_{ij}^{2}$ is a measure of the angular separation between two constituent’s momenta scaled by the jet radius parameter $R$.555Since we start with the constituents of an identified jet, all of the particles are (re)clustered into a single tree. For this reason, the single-particle distance measure and the jet radius parameter $R$ in the (re)clustering algorithm are irrelevant. For any value of $p$, generalized-$k_{t}$ provides a pairwise clustering structure that directly mimics eq. (18). For $pp$ collisions, one insteads use a form of eq. (22) based on transverse momenta $p_{T}$ and distance $\Delta R_{ij}$ in azimuthal angle and rapidity.

From this clustering tree, one can determine the observable $x$ by applying the recursion relation $\hat{x}(z,x_{1},x_{2})$ at each stage of the clustering. Specifically, the value of $x$ at each node depends on the momentum fraction $z$ given by the $2\to 1$ merging kinematics as well as on the $x_{1}$ and $x_{2}$ values determined from the corresponding daughter nodes (which might be the initial weights $w_{a}$). When all nodes are contained in a single connected tree, the root node represents the entire jet, and the root value of $x$ determines the final observable.

Even though the clustering tree is IRC safe, the resulting fractal observable $x$ is generally collinear unsafe. These collinear divergences are absorbed into the GFFs, and are in fact responsible for the evolution in eq. (3).

3.2 Requirements

There are a few fundamental limitations on the choice of $\hat{x}(z,x_{1},x_{2})$ dictated by the fact that this same function will appear in eq. (18). First, the recursion relation must be symmetric under the exchange $z\leftrightarrow 1-z$, $x_{1}\leftrightarrow x_{2}$, since the assignment of these labels is unphysical.666In the case of jets with heavy flavor, one could use heavy-flavor tags to define asymmetric recursion relations (see e.g. Ilten:2017rbd ). We do not give a separate treatment of heavy-flavor GFFs in this work, and instead assume to always work in the $m_{b,c}\ll E_{\rm jet}R$ limit. Second, the recursion relation has to be IR safe, since the GFF formalism only regulates collinear (and not soft) divergences. In order that an emission with $z\to 0$ does not change the observable, IR safety translates into the conditions

[TABLE]

such that an arbitrarily soft branch in the clustering tree has no impact on the values of $x$. Third, the recursion relation has to have unambiguous limits. As a counterexample, $\hat{x}(z,x_{1},x_{2})=x_{1}^{z}x_{2}^{1-z}$ satisfies eq. (23) when $x_{1}$ and $x_{2}$ are non-zero, but not when they vanish. Apart from these limitations, any choice of $\hat{x}(z,x_{1},x_{2})$ (along with starting weights and a clustering tree) defines a fractal observable.

The tree traversal prescription, along with the requirement in eq. (23), helps ensure IR safety to all $\alpha_{s}$ orders. As a counterexample, consider the sum over all tree nodes of some function $f(z)$ which vanishes as $z\rightarrow 0$ or $z\rightarrow 1$. In that case, the resulting observable would receive no contribution from a single infinitely soft splitting, but subsequent finite $z$ splittings that followed the soft one would not be suppressed, violating IR safety. By contrast, eq. (23) requires the contribution from an entire soft branch to be suppressed, as desired.

In this paper, we mainly focus on recursion relations that do not depend explicitly on the opening angle $\theta$ between branches in the clustering tree. In sec. 7, we do discuss how the recursion relation gets modified if a threshold value for $\theta$ is introduced (i.e. $\theta_{\rm thr}=R_{\rm sub}\ll R$). Of course, fractal observables depend indirectly on angular information through the structure of the clustering tree, but as discussed below, the leading-order evolution equations do not depend on the clustering algorithm. When explicit $\theta$-dependence is included in the $\hat{x}$ function, this sometimes results in a fully IRC-safe observable, requiring a different type of evolution equation that is beyond the scope of the present work (see e.g. HarvardInProgress ).

3.3 Evolution Equations

The generalized-$k_{t}$ clustering tree has an obvious mapping to a parton branching tree, such that at order $\alpha_{s}$, the RG evolution is given precisely by eq. (18), with the flavor of the GFF matching the flavor of the jet’s initiating parton. More formally, as discussed in sec. 2.3, the NLO calculation of the bare GFF shows that the same recursion relation $\hat{x}(z,x_{1},x_{2})$ appears in eq. (2.3), as desired.

In fact, to order $\alpha_{s}$, the evolution in eq. (18) is insensitive to the clustering tree, as long as it is IRC safe, even if the fractal observable itself depends on the clustering order. We explicitly test this surprising feature in sec. 5. Note that if the clustering tree is not collinear safe, in the sense that particles with collinear momenta are not clustered with each other first, then the collinear divergences in the GFF will not cancel against the collinear divergences in the hard matching coefficients of eq. (17). If the clustering tree is not IR safe, then the observable $x$ is not IR safe, and the GFF formalism does not apply.

We stress that the evolution in eq. (18) is only valid to lowest order in $\alpha_{s}$. At higher orders in $\alpha_{s}$, the evolution of fractal observables is more complicated, but, as discussed more in the paragraph below, still satisfies the property that the evolution of $\mathcal{F}_{i}(x,\mu)$ depends only on GFFs of the same observable. Schematically, this can be written as

[TABLE]

where $\otimes$ represents a convolution. This equation includes $1\rightarrow n$ splittings at order $\alpha_{s}^{n-1}$. There is no longer a one-to-one correspondence between pairwise clustering trees and GFF evolution trees, and one has to explicitly carry out the calculation in eq. (2.3) to higher orders to determine the evolution. In particular, there will be different clusterings of the $1\to n$ splitting into a binary tree when integrating over phase space, which depend on the choice of clustering algorithm. Because our specific realization of fractal observables in this section is based on recursive clustering trees, this guarantees that eq. (24) depends only on GFFs of the same type as $\mathcal{F}_{i}$ at all perturbative orders.

To justify the structure of eq. (24) in a bit more detail, it is instructive to take a closer look at the $1/\epsilon_{\rm UV}$ poles of $\mathcal{F}_{i}$. As usual, the anomalous dimension of the GFFs is determined by the single $1/\epsilon_{\rm UV}$ poles. At order $\alpha_{s}$, we get $(1/\epsilon_{\rm UV})P_{i\to jk}$, as shown in eq. (2.3). At order $\alpha_{s}^{2}$, the $1\to 3$ splitting factorizes into a sequence of two $1\to 2$ splittings when the angles of the splittings are strongly ordered. This leads to a term like $(1/\epsilon_{\rm UV}^{2})P_{i\to jk}\otimes P_{j\to\ell m}$ which does not contribute to the GFF’s anomalous dimension. However, it does justify attaching $\mathcal{F}_{j}$ and $\mathcal{F}_{k}$ to the external splittings in eq. (18), as it corresponds to the cross term between a one-loop renormalization factor and one-loop $\mathcal{F}_{j}$ (and tree-level $\mathcal{F}_{k}$). Away from the strongly-ordered limit, the $1\to 3$ splitting does have a genuine $1/\epsilon_{\rm UV}$ divergence, contributing to the second term in eq. (24). The precise structure of this term depends on how the clustering algorithm maps the three partons to a binary tree. The justification for attaching GFFs to each of the three external partons follows again by considering higher-order corrections with some strong ordering. For example, consider a $1\to 5$ splitting that is strongly ordered such that it factorizes in a $1\to 3$ splitting, in which two partons undergo $1\to 2$ splittings. Such a term would have a $1/\epsilon_{\rm UV}^{3}$ divergence, corresponding to the cross term of the renormalization factor for the $1\to 3$ splitting term at order $\alpha_{s}^{2}$ with two one-loop $\mathcal{F}$’s and one tree-level $\mathcal{F}$. Finally, the $1/\epsilon_{\rm UV}$ from the one-loop virtual contribution to the $1\to 2$ splitting gives a higher-order correction to the first term in eq. (18). For the remainder of this paper, we focus on the leading-order evolution, leaving an analysis at higher orders to future work.

4 Weighted Energy Fractions

The procedure outlined in sec. 3 is very general, but for special choices of $\hat{x}(z,x_{1},x_{2})$, the definition of a fractal observable can simplify greatly. In this section, we consider the recursion relation

[TABLE]

where $\kappa>0$ is an energy exponent. As we will see, for any choice of pairwise clustering tree, the resulting observable simplifies to a sum over the hadrons in a jet,

[TABLE]

where $\kappa$ is the same as in eq. (25), and $w_{a}$ is the hadron weight factor. We call these observables weighted energy fractions.

Several examples of weighted energy fractions have already been studied in the literature. The weighted jet charge is defined for any $\kappa>0$ and weights given by the electric charges of final-state hadrons FIELD19781 ; Krohn:2012fg ; Waalewijn:2012sv . This quantity has, for example, been used in forward-backward asymmetry measurements at $e^{+}e^{-}$ experiments Stuart:1989db ; Decamp:1991se , as well as to infer the charge of quarks Berge:1979qg ; Buskulic:1992sq ; Abazov:2006vd . Recently, the scale dependence of the average jet charge was observed in $pp\to$ dijets Aad:2015cua . Track fractions correspond to the case of $\kappa=1$, where charged particles are given weight 1 and neutral particles given weight 0 Chang:2013rca ; Chang:2013iba . Jet $p_{T}^{D}$ is a weighted energy fraction with $\kappa=2$ and all particles given weight 1 Chatrchyan:2012sn ; Pandolfi:1480598 . Weighted energy fractions with arbitrary $\kappa>0$ and $w_{a}=1$ for all particles were studied in LABEL:Larkoski:2014pca for applications to quark/gluon discrimination.

4.1 Associativity

Weighted energy fractions have an associative recursion relation, meaning that the order of the clustering tree does not affect the final observable. To see this, consider the case of just three particles with weights $\{w_{1},w_{2},w_{3}\}$ and respective momentum fractions $\{z_{1},z_{2},z_{3}\}$. As shown in fig. 3, there are three clustering trees that can be built using only $1\rightarrow 2$ splittings, labeled as A, B, and C.777Of course, for a specific choice of kinematics, not all of these trees will be possible from generalized-$k_{t}$ clustering, particularly in the collinear limit. The corresponding observables are

[TABLE]

Using eq. (25) and the fact that $z_{1}+z_{2}+z_{3}=1$, it is straightforward to prove that

[TABLE]

owing to the fact that the recursion relation has homogenous scaling with $z$. This argument generalizes to an arbitrary numbers of particles, so the weighted energy fractions are indeed independent of the clustering tree.888Remember that this tree is one obtained from reclustering the particles in the jet. The value of a jet observable of course depends on the choice of initial jet algorithm, which may itself be a clustering algorithm.

Of course, there are other observables that have non-associative recursion relations, where the observable does not simplify to a sum over final-state hadrons and the full tree traversal is necessary. We explore some non-associative observables in sec. 5.

4.2 Extraction of GFFs

In general, to extract GFFs, one has to numerically match the cross section in eq. (17) using perturbatively calculated values for the coefficients $C_{i}$, $\mathcal{J}_{i\to j}$, $\mathcal{J}_{i\to jk}$, …. For the parton shower studies in this paper, we limit ourselves to leading order where $C_{i}^{(0)}(z_{J})=\delta(1-z_{J})$, and we use parton-shower truth information to assign the parton label $i$. To generate pure samples of quark- and gluon-initiated jets, we use the $e^{+}e^{-}\to\gamma/Z^{*}\to q\bar{q}$ and $e^{+}e^{-}\to H^{*}\to gg$ processes in Pythia 8.215 Sjostrand:2014zea , switching off initial-state radiation. We find jets using FastJet 3.2.0 Cacciari:2011ma , with the $ee$-generalized $k_{t}$ algorithm with $p=-1$ (i.e. the $e^{+}e^{-}$ version of anti-$k_{t}$ Cacciari:2008gp ) and then determine the various weighted energy fractions on the hardest jet in the event. At leading order, the normalized probability distributions for the weighted energy fractions directly give the corresponding GFF $\mathcal{F}_{i}(x,\mu)$.

As discussed in sec. 2.2, for jets of a given energy $E_{\rm jet}$ and radius $R$, the characteristic scale for GFFs is expected to be

[TABLE]

which is roughly the scale of the hardest possible splitting in the jet. By varying $E_{\rm jet}$ and $R$ but keeping $\mu$ fixed, we can estimate part of the uncertainty in the extraction of the GFFs. In addition, we assess the uncertainty from using different parton shower models. Here, since our primary interest is in the perturbative uncertainty in different shower evolution equations, we test the native Pythia parton shower along with the Vincia 2.0.01 Giele:2007di and Dire 0.900 Hoche:2015sya parton shower plugins. A further source of uncertainty would be given by the hadronization model, which enters the boundary conditions used for GFF evolution. This is not included in our present study, since we decided to interface all of the showers above with the Lund string model. In the context of an experimental analysis, one would also have statistical and systematic uncertainties from the extraction of GFFs from data.

For each observable $x$, there are 11 GFFs, corresponding to 5 quark flavors $\{u,d,s,c,b\}$, 5 anti-quark flavors, and the gluon. To avoid a proliferation of curves, it is convenient to define singlet (denoted by $\langle\rm Quark\rangle$ in the figures below) and non-singlet combinations for the quark GFFs, respectively,

[TABLE]

For the observables we study, the anti-quark GFFs are either identical to the quark GFFs or simply involve the replacement $x\to-x$, due to charge conjugation symmetry. We start by showing numerical results for the gluon GFF and the quark-singlet combination, postponing a discussion of the non-singlet case to sec. 4.5.

In fig. 4, we show the extracted gluon and quark-singlet GFFs at $\mu=E_{\rm jet}R=100~{}$ for the weighted energy fractions with $w_{a}=1$, comparing $\kappa=0.5$ and $\kappa=2$. Since gluon jets have roughly a factor of $C_{A}/C_{F}$ larger hadron multiplicity than quark jets, the mean of the gluon GFF is roughly a factor of $(C_{A}/C_{F})^{1-\kappa}$ higher than the mean of the quark-singlet GFF. In the left column, we show the impact of changing the jet radius $R=\{0.3,0.6,0.9\}$, leaving $\mu$ fixed. The envelope from changing $R$ is very small, indicating that $\mu=E_{\rm jet}R$ is an appropriate definition for the RG scale. In the right column, we show the impact of switching between the Pythia, Vincia, and Dire parton shower models. The envelope is larger, but still reasonably narrow, giving us confidence in the extraction of the GFFs, at least as far as changing the perturbative shower model is concerned. Though not shown here, we checked that the GFFs for the $\kappa\to 1$ and $\kappa\to\infty$ limits behave sensibly as well (see sec. 4.4 below).

4.3 Evolution of GFFs

We now use these extracted GFFs as boundary conditions for the RG evolution in eq. (18). In app. C, we describe in detail the numeric implementation of the evolution. Formally, the evolution equations work equally well running up or down in $\mu$, but in practice downward evolution is numerically unstable, as discussed further in app. D. As a proof of principle for our RG evolution code, we show upward evolution from $\mu=100~{}$ to $\mu=4~{}$, comparing our RG evolution in eq. (18) to that obtained from parton showers.

In figs. 5 and 6, we present the evolution results for gluon and quark-singlet GFFs respectively, for the weighted energy fractions with $\kappa=\{0.5,1.0,2.0\}$. We test three different choices for the particle weights: $w_{a}=1$ for all particles, $w_{a}=1$ ($w_{a}=0$) for charged (neutral) particles, and $w_{a}=Q_{a}$ with $Q_{a}$ being the particle’s electric charge. The initial conditions extracted from the parton showers at $\mu=100$ GeV are the same as those shown in fig. 4, with the same color scheme of red for gluon GFFs and blue for quark-singlet GFFs. As described in sec. 4.2, the uncertainty bands are given by the envelope of values obtained both from varying the jet radius/energy (keeping $\mu$ fixed) and from using different parton showers. The evolved distributions to $\mu=4$ TeV are shown in orange for the gluon GFFs and light blue for the quark-singlet GFFs, where the uncertainty bands show the spread in final values due to the spread in initial conditions.

For comparison, we show in dashed lines the GFFs extracted at $\mu=4$ TeV, averaged over the three parton showers and three $R$ values.999The uncertainties from varying the jet radius/energy and changing parton showers at $\mu=4$ TeV are similar to the ones shown at $\mu=100$ GeV. Overall, our numerical GFF evolution agrees well with parton shower evolution, with both methods giving the same shift in the peak locations. As previously seen in LABEL:Waalewijn:2012sv, the two evolution methods agree best for $\kappa\geq 1$, with larger differences seen in the widths of the distributions when $\kappa<1$. This is likely because $\kappa<1$ is more sensitive to collinear fragmentation, with larger expected corrections from higher-order perturbative effects. Note the expected $\delta$-function when $\kappa=1$ and $w_{a}=1$ for all particles, since the sum of the energy fractions for all particles in the jet equals 1. The $\kappa\to 1$ limit of weighted energy fractions is discussed in sec. 4.4 below.

4.4 Limits

There are a few interesting limits of the weighted energy fractions. For the case of $\kappa=0$, the energy fractions $z_{a}$ drop out, so $x$ simply counts the hadrons in the final state, weighted by $w_{a}$. Although hadron multiplicity is IR unsafe, it is possible to calculate the evolution of the average hadron multiplicity using fragmentation functions, see e.g. refs. Malaza:1985jd ; Lupia:1997in ; Bolzoni:2012ii . This case requires special care, however, because of the soft gluon singularity of the splitting functions. IR-safe variants of multiplicity that have only collinear singularities are explored in a forthcoming paper HarvardInProgress .

For the case of $\kappa=1$ with all hadrons assigned weight 1, the weighted energy fraction simply becomes $x=\sum_{a}z_{a}=1$. Still, we can expand around the $\kappa\to 1$ limit to find a non-trivial observable Larkoski:2014pca . Consider the modified weighted energy fraction and its limit,

[TABLE]

In the limiting case, the recursion relation becomes

[TABLE]

with initial hadron weights of $w_{a}=0$ (due to the $-1$ in eq. (31)). This is easy to verify by testing the three clustering trees in fig. 3.101010Amusingly, the recursion relation in eq. (32) is associative for any choices of initial hadron weights, leading to the fractal observable $x=\sum_{a\in\text{jet}}z_{a}(w_{a}+\ln z_{a})$.

The behavior of the evolved GFFs in the $\kappa\rightarrow 1$ limit offers a non-trivial cross check of our evolution code. Away from the limiting value, the RG evolution can be implemented using the recursion relation in eq. (25). At the limiting value, we have to use a different RG evolution based on the recursion relation in eq. (32). The smooth convergence of the evolved distributions as $\kappa\rightarrow 1$ is illustrated in fig. 7a, showing the modified weighted energy fraction from eq. (31). The solid curves show the extraction of the corresponding GFFs at $\mu=100$ GeV with $\kappa=0.99$ and $\kappa=1.01$, which correctly bracket the $\kappa\to 1$ limit.111111In practice, we first extract the $\kappa=0.99$ and $\kappa=1.01$ distributions for the unmodified weighted energy fraction, and then do a simple change of variables to match the definition in eq. (31). The dashed curves show the evolution to $\mu=4$ TeV, where there is again a smooth approach to $\kappa\to 1$.

In the limit that $\kappa\rightarrow\infty$, the most energetic hadron in the jet dominates the sum in eq. (26). We can then take the $\kappa$-th root of the weighted energy fraction to have a smooth $\kappa\to\infty$ limit:

[TABLE]

where the maximum is only taken over particles with non-zero weights. The corresponding recursion relation is

[TABLE]

with modified initial hadron weights of $\tilde{w}_{a}=|\text{sign}(w_{a})|=\{0,1\}$. For these modified weights, it is easy to verify that eq. (34) gives an associative recursion relation using fig. 3.121212It is also possible to keep track of the signs of hadron weights by using an alternative recursion relation $\hat{x}(z,x_{1},x_{2})=\text{signed-max}(zx_{1},(1-z)x_{2})$, where the signed-max function takes the term (positive or negative) with the largest absolute value. Again, this only yields an associative recursion if the hadron weights are a constant multiple of $\{-1,0,1\}$.

In fig. 7b, we show the approach to $\kappa\rightarrow\infty$ for the gluon GFFs, considering the case of all particles with equal weight $w_{a}=1$. Here, the finite-$\kappa$ evolution equations use the recursion relation in eq. (25) while the $\kappa\to\infty$ limit uses eq. (34), and we plot the $\kappa$-th root of the weighted energy fractions as given in eq. (33). Both the extracted distributions at $\mu=100$ GeV and the evolved distributions to $\mu=4$ TeV show a smooth transition from $\kappa=4$ to $\kappa=6$ to the final $\kappa\to\infty$ limit. This is again a non-trivial cross check of our evolution code.

4.5 Moment Space Analysis

To gain further insight into the evolution of the GFFs, it is instructive to examine the evolution equations for the first two moments, which are related to averages and widths of the distribution for the fractal observable. In general, the moments of a GFF are defined as

[TABLE]

with $N\geq 0$. For the specific case of the weighted energy fractions, it is convenient to introduce a transformed version of the splitting functions

[TABLE]

Integrating eq. (18) against $x^{N}$, the moment space evolution equation for a weighted energy fraction is

[TABLE]

where it is crucial that $N$ is an integer. A derivation of this expression is given in app. E.

These evolution equations are more compact in the color singlet/non-singlet basis introduced in eq. (4.2). For the quark-non-singlet pieces, the evolution of the first moment (i.e. the mean) is given by

[TABLE]

Since $\overline{P}_{q\rightarrow qg}(\kappa)<0$ for all positive $\kappa$, eq. (38) implies that the averages of the different (anti-)quark GFFs functions converge to a common value as $\mu$ evolves upward. This behavior is expected, since QCD branchings only depend on the parton’s color charge, so the low-scale differences between the (anti-)quark flavors, due to e.g. electric charge, get washed out at high scales.

The quark-singlet combination mixes with the gluon GFF. For the first moment this is given by

[TABLE]

As shown in fig. 8, the matrix in eq. (39) always has one negative eigenvalue for all $\kappa$, which implies that the first moment of the quark-singlet GFF tries to track the first moment of the gluon GFF. For example, in the case of $\kappa=1$, the combination $2C_{F}\,\overline{\mathcal{S}}(1,\mu)-n_{f}T_{F}\,\overline{\mathcal{F}}_{g}(1,\mu)$ asymptotes to zero at high $\mu$. The second eigenvalue has different signs depending on the value of $\kappa$. For $\kappa<1$, it is positive, so the first moments of both the quark-singlet and gluon GFF increase with $\mu$. For $\kappa>1$, the second eigenvalue is negative, so the first moments decrease with $\mu$. For the special case $\kappa=1$, the second eigenvalue is zero, and the corresponding eigenvector $\overline{\mathcal{S}}(1,\mu)+\overline{\mathcal{F}}_{g}(1,\mu)$ stays constant with $\mu$. These broad features agree with the behaviors already seen in figs. 5 and 6.

Turning to the second moments, the non-singlet evolution is

[TABLE]

Since the splitting function in the first term is negative for all values of $\kappa$, this term pushes the second moment of the non-singlet GFFs towards zero as well. Note, however, that the splitting function in the second term has the opposite sign. For the weighted energy fractions with $\kappa>1$, which have $\overline{\mathcal{F}}_{g}(1,\mu)\rightarrow 0$ as $\mu\rightarrow\infty$, this second term is not important, so the different quark GFFs asymptote to the same second moment. For the weighted energy fractions with $\kappa\leq 1$, however, this is not the case. As shown below in fig. 9d for $\kappa=0.5$, the growth of $\overline{\mathcal{F}}_{g}(1,\mu)$ outpaces the decrease in $\overline{\mathcal{N}}_{ij}(1,\mu)$ from the first term, which leads to differences in the widths (but not the means) of the different quark GFFs.

Assuming the asymptotic behavior $\overline{\mathcal{N}}_{ij}(1,\mu)\to 0$ for simplicity, the evolution of the second moments of the quark-singlet and gluon GFF can be written as

[TABLE]

where the assumption allows us to write the nonlinear term as a function of $\overline{S}(1,\mu)$ instead of individual (anti-)quark contributions.131313For the weighted energy fractions in this study, this approximation is very accurate, giving corrections at the per-mille level. Due to this nonlinear behavior, we now resort to a numerical analysis.141414Alternatively, if one assumes that $\overline{\mathcal{F}}_{g}(1,\mu)$ and $\overline{S}(1,\mu)$ have also reached their asymptotic behavior, the equation becomes linear again. This approximation turns out to be even more accurate than the assumption $\overline{\mathcal{N}}_{ij}(1,\mu)\rightarrow 0$, though for our numerical studies, we make no simplifications.

In fig. 9, we show an example of the RG evolution of the first and second moments of the gluon GFFs, quark-singlet GFFs, and $u$-$d$ quark-non-singlet GFFs. Here, we consider weighted energy fractions where charged particles have weight 1 and neutral particles have weight zero, comparing $\kappa=0.5$, $1$, and $2$. The evolution starts from GFFs extracted at $\mu=100$ GeV, as described in sec. 4.2. The GFF moments are then evolved up to $\mu=10^{7}$ GeV using the equations above.151515We checked that this agrees with first evolving the full binned distributions and then calculating the first and second moments. To connect with the plots in figs. 5 and 6, we also indicate the first (second) moments extracted from the parton shower average at $\mu=4$ TeV with dots (diamonds).

As expected, the first moments evolve in the direction predicted by the eigenvalues in fig. 8, with the $\kappa<1$ first moment moving to larger values as $\mu$ increases, and the $\kappa>1$ first moment moving to smaller values. For the boundary case of $\kappa=1$, the first moment of the gluon and quark singlet GFFs move toward each other, leaving their sum fixed. The second moments roughly evolve in the same direction as first moments, though with different rates. The exception is the $\kappa=1$ second moment, where both the gluon and quark singlet values decrease (very slowly), as seen already in figs. 5e and 6e. The first moment of the non-singlet GFFs approaches zero, as indicated by $\overline{P}_{q\rightarrow qg}(\kappa)<0$. The second moments behave as discussed above, decreasing for $\kappa=1$ and $\kappa=2$, and increasing for $\kappa=0.5$ since $\overline{\mathcal{F}}_{g}(1,\mu)$ grows very large.

We could continue our analysis to third and higher moments, which is a standard way to efficiently solve the DGLAP equations. An interesting difference with the evolution of the ordinary FFs is that we only get the simple expression in eq. (37) for integer moments. In addition, the simple form of eq. (37) does not hold for general fractal observables with more complicated recursion relations. For these reasons, we only show the evolution of the first two moments here. Brief moment-space analyses for the non-associative observables in sec. 5 are given in app. E.

5 Tree-Dependent Observables

We now study fractal jet observables that do depend on the choice of clustering tree. These are also called non-associative observables, since $x_{\rm A}\not=x_{\rm B}\not=x_{\rm C}$ in the notation of eq. (27). We start in sec. 5.1 with node-product observables, where the recursion relation simplifies to a sum over internal nodes of the tree. We then turn to a more general family of non-associative observables in sec. 5.2.

5.1 Node Products

Node-product observables are based on the recursion relation

[TABLE]

Note that the last term in eq. (51) is independent of $x_{1}$ and $x_{2}$, and the factor of 4 is added for convenience, to normalize the contribution of a balanced splitting with $z=1/2$ to be 1. It is straightforward to check that this recursion relation is not associative for generic values of $\kappa$, by considering the three-particle trees in fig. 3. For the special case of $\kappa=2$, the recursion relation is associative, yielding an observable closely related to $p_{T}^{D}$ (i.e. the weighted energy fraction with $\kappa=2$),

[TABLE]

For generic values of $\kappa$, this recursion relation simplifies to a sum over the leaves and nodes in the binary tree,

[TABLE]

where $z_{L,R}$ are the momentum fractions carried by the two branches at this node, relative to the whole jet (i.e. $z_{L}+z_{R}\not=1$).161616If $z_{L}$ and $z_{R}$ had been relative to the node instead, this observable would not be IR safe, as the contribution from an arbitrary soft gluon that subsequently splits collinearly would not be suppressed; see eq. (23). To see how this simplification arises, note that the $(4z(1-z))^{\kappa/2}$ term in eq. (51) adds the product of branch energy fractions to the observable; the $x_{1}\,z^{\kappa}$ and $x_{2}\,(1-z)^{\kappa}$ terms then rescale the energy product to the whole jet momentum. In this way, node products have intermediate complexity between the weighted energy fractions (with no tree dependence) and more general observables (where the full tree recursion is required).

For simplicity, we focus on the case with starting weights of $w_{a}=0$, such that the node-product observable only depends on non-leaf nodes, as advertised in eq. (5). In fig. 10, we show the distributions for the gluon GFFs for the node products extracted from Vincia at a jet scale of $\mu=100$ GeV. Here, we take $\kappa=\{1,2,4\}$, testing three different values of the generalized-$k_{t}$ clustering exponent $p=\{-1,0,1\}$. The tree dependence of this observable for $\kappa=1$ and $\kappa=4$ is evident. This is particularly true for $\kappa=4$, where the spikes near $x=1.1$ (and $x=0.8$) come from balanced splittings that are more prevalent in $k_{t}$ trees than C/A or anti-$k_{t}$ trees. For $\kappa=2$, the node-product observable is independent of $p$, since it is identical to the associative observable $2(1-p_{T}^{D})$, as shown in eq. (52).

Observables measured on anti-$k_{t}$ clustering trees tend to be qualitatively distinct from observables measured on $p\geq 0$ trees. This is expected, because C/A and $k_{t}$ trees are constructed according to angular and $k_{t}$ ordering, respectively, so these observables more directly mirror the singularity structure of QCD and the expected dynamics of the parton shower. By contrast, anti-$k_{t}$ trees have a hybrid ordering where angles tend to go from small to large, but energies tend to go from large to small. Indeed, this reversal in the energy ordering is reflected in fig. 10, where the product $z_{L}z_{R}$ tends to be smaller for anti-$k_{t}$ trees, leading to larger (smaller) values of node-product observable for $\kappa=1$ ($\kappa=4$). Because of this hybrid anti-$k_{t}$ ordering, one might expect higher-order perturbative corrections to be more important for $p<0$ when evolving the GFFs, but this can only be confirmed by doing an explicit calculation, which is beyond the scope of the present work.

Despite the fact that different values of $p$ lead to different observables, the leading-order evolution equations are independent of $p$. To check whether this is a sensible feature, we evolve the gluon GFFs in fig. 11 for node products with $\kappa=\{1,4\}$ and $p=\{-1,0,1\}$. The uncertainty bands in fig. 11 are obtained from the variation of jet radius $R=\{0.3,0.6,0.9\}$ and parton shower PS $=\{\textsc{Vincia},\textsc{Pythia},\textsc{Dire}\}$, as described in sec. 4.2. If the evolution from 100 GeV to 4 TeV would perfectly agree with the extraction at 4 TeV, this would confirm that the evolution is independent of $p$ and all $p$ dependence resides in the initial conditions. Although the agreement is not perfect, the amount of agreement between the evolution from 100 GeV to 4 TeV and the extraction at 4 TeV seems to be fairly independent of $p$, suggesting that this is a reasonable first approximation. Given the interesting features in the node-product observables as a function of scale, this motivates both higher-order calculations of their RG evolution, as well as measurements in data.

For completeness, we show the evolution of the first and second moments for the node-product observables in app. E.2.

5.2 Full-Tree Observables

As our final example of a fractal observable, we present a recursion relation that depends on the full structure of the clustering tree,

[TABLE]

This recursion relation satisfies the requirements in eq. (23), making this observable IR (but not collinear) safe. Eq. (54) defines a family of fractal observables which depend on the initial particle weights $w_{a}$, the generalized-$k_{t}$ clustering exponent $p$, and the parameters $\kappa$ and $\xi$. We know of no alternative way to calculate this observable apart from performing the full leaf-to-root recursive traversal of the clustering tree. Of course, for the special value of $\xi=0$, these observables become weighted energy fractions.

The tree dependence of this observable is illustrated in fig. 12 for $\kappa=2$ and $\xi=\{-2,0,2\}$, where charged particles are given weight 1 and neutral particles weights 0. For nonzero $\xi$, we see that the GFFs depend on the choice of $p$, with rather different behaviors for anti-$k_{t}$ compared to $k_{t}$ and C/A. The (associative) observables plotted in fig. 12b are equivalent to the weighted energy fraction with the same weights and $\kappa=2$, shown on this plot for comparison. Corresponding results for the evolution of the gluon GFFs are shown in fig. 13. In this case, it is much clearer that the amount of agreement between the evolution from 100 GeV to 4 TeV and the extraction at 4 TeV is independent of $p$. Thus, the fact that the leading-order RG evolution is independent of $p$ seems reasonable, even though the GFFs themselves are tree dependent. This is highlighted by fig. 13d, where the double hump structure at 100 GeV is smoothed out both by the RG evolution equations and the parton shower.

Again for completeness, we discuss the evolution of the first two GFF moments for these full-tree observables in app. E.3.

6 Application in Quark/Gluon Discrimination

Robust and efficient discrimination between quark- and gluon-initiated jets is a key goal of the jet substructure community Abdesselam:2010pt ; Altheimer:2012mn ; Altheimer:2013yza ; Adams:2015hiv , with applications both in searches for physics beyond the SM and precision tests of QCD (see further discussions in Gallicchio:2011xq ; Gallicchio:2012ez ; Larkoski:2013eya ; Larkoski:2014pca ; Bhattacherjee:2015psa ; Badger:2016bpw ; FerreiradeLima:2016gcz ; Komiske:2016rsd ; Davighi:2017hok ; Gras:2017jty ). Weighted energy fractions are already used for quark/gluon discrimination, specifically the $p_{T}^{D}$ observable Pandolfi:1480598 ; Chatrchyan:2012sn used by CMS in its quark-gluon likelihood analysis CMS:2013kfa . Here, we explore the potential discrimination power of non-associative fractal jet observables, corresponding to non-associative variants of $p_{T}^{D}$. An alternative application of the GFF formalism to quark/gluon discrimination will be presented in LABEL:HarvardInProgress.

It is not immediately obvious that non-associativity should be a valuable feature to help distinguish quark- from gluon-initiated jets. Compared to $p_{T}^{D}$, non-associative observables are of course sensitive to the angular structure of the jet through the clustering tree. Then again, discriminants like the (generalized) angularities Berger:2003iw ; Almeida:2008yp ; Ellis:2010rwa ; Larkoski:2014pca and energy correlation functions Larkoski:2013eya also encode angular information about particles in the jet, either their angular distance to the jet axis or their pairwise angular distance to each other. As we will see, there are non-associative observables that do exhibit better performance than $p_{T}^{D}$, at least in the context of a parton shower study, but we do not (yet) understand the origin of that improvement from first principles.

Here, our primary interest in non-associative observables is for testing the evolution of quark/gluon discrimination power as a function of RG scale $\mu$. As recently studied in refs. Badger:2016bpw ; Gras:2017jty , different parton showers exhibit different quark/gluon discrimination trends as a function of jet energy. Therefore, the study of fractal jet observables might help identify which higher-order effects in the parton shower are most important for correctly modeling the radiation patterns of quarks and gluons.

As an initial investigation into non-associative fractal observables for quark/gluon discrimination, we consider some examples of the node-product and full-tree observables from sec. 5. In fig. 14, we show two good quark/gluon discriminants, comparing the gluon GFF distribution to the quark-singlet GFF distribution. We also show the down-quark and bottom-quark GFFs as a cross check. An example of a node-product observable from eq. (53) is shown in fig. 14a, where we take $\kappa=1$ and $w_{a}=0$ on a C/A tree. An example of a full-tree observable from eq. (54) is shown in fig. 14b, where we take $\kappa=2$ and $\xi=4$ on a C/A tree with all particles given weight 1. There are noticeable differences between the gluon and quark-singlet GFFs which can be exploited for the purposes of discrimination. Among the observables we tested, these two performed among the best, outperforming, for example, variants using only charged particles.

To evaluate the potential quark/gluon discrimination power more quantitatively, we show ROC (receiver operating characteristic) curves showing the efficiency of identifying quark jets against the mistag rate for gluon jets. These plots are obtained from Vincia, comparing the discrimination performance at $\mu=100$ GeV to $\mu=4$ . In fig. 15, we show variants of the node-product observables defined on C/A trees for $\kappa=\{1,2,4\}$, recalling that $\kappa=2$ is the same as $2(1-p_{T}^{D})$. The node product with $\kappa=1$ exhibits much better discrimination power than $\kappa=2$, especially at $\mu=4$ . The discrimination power does continue increasing (slowly) with lower $\kappa$, but approaching the $\kappa\rightarrow 0$ limit, the observable becomes IR unsafe and the GFF formalism no longer applies.

We can check whether this jet-energy dependence is reasonable using the RG evolution equations, as shown in fig. 16. For $\kappa=1$, the discrimination power does indeed increase with increasing $\mu$, but not as much as predicted by the parton showers. This could have already been anticipated from the results in fig. 11b, where the RG-evolved gluon GFF does not shift as dramatically as predicted in the parton showers. This could either be a sign that the parton showers are too aggressive in their evolution, or that higher-order terms in the evolution equation are important for getting the proper shape of the $\kappa=1$ distribution. For $\kappa=2$, the evolution of the ROC curves according to eq. (18) does match the evolution in the parton shower, but this evolution is very slight, less than the spread in the ROC curves at either scale from varying $R$ and the parton shower. For $\kappa=4$, the discrimination power is poor at all scales, but the evolution matches well between eq. (18) and the parton showers.

We next turn to the full-tree observables in fig. 17, using a C/A tree with $\kappa=2$ on all particles. We compare $\xi=\{0,2,4,6\}$, where $\xi=0$ is identical to $p_{T}^{D}$. The $\xi=4$ observable yields comparable performance to $p_{T}^{D}$ at $\mu=100$ GeV, but performs somewhat better than $p_{T}^{D}$ at $\mu=4$ TeV. Note that the quark/gluon discrimination power is not monotonic as a function of $\xi$. We can again check whether this evolution is reasonable using the RG equations, as shown in fig. 18. For all three $\xi$ values, the evolution of the ROC curves in eq. (18) matches the parton shower, but the evolution is extremely slow.

As emphasized in LABEL:Larkoski:2014pca, predicting the quark/gluon discrimination power from first principles is a much more challenging task than predicting the distributions themselves. Because the ROC curve shapes depend sensitively on the overlap between the quark and gluon distributions, small changes in the distribution shapes can lead to large changes in the predicted discrimination power. This is especially evident in fig. 18, where the uncertainties in the ROC curves at the same scale are generally larger than the evolution between scales. This highlights the importance of precision calculations for correctly predicting quark/gluon discrimination behavior.

7 Fractal Observables from Subjets

As our final investigation into the structure of fractal jet observables, we now consider the possibility that the recursion relation in eq. (1) is modified to depend on the angular scale of the clustering. For simplicity, we only consider observables defined on angular-ordered C/A clustering trees, since in that case the depth in the C/A tree is directly associated with an angular scale $\theta$. This opens up the possibility to define a modified recursion relation with $\theta$ dependence, for example,

[TABLE]

As shown in fig. 19, the nodes as defined by $\hat{x}_{1}$ become the starting weights for the subsequent nodes defined by $\hat{x}_{2}$.

It is straightforward to implement the leading-logarithmic resummation of an observable defined by eq. (57). Starting from a low-energy boundary condition, this involves an initial evolution to the scale

[TABLE]

using eq. (18) with the recursion relation $\hat{x}_{1}$, followed by an evolution to $\mu=E_{\rm jet}R$ using $\hat{x}_{2}$ instead. The discontinuity in anomalous dimensions of the evolution equations across the threshold $\mu_{\rm sub}$ will be compensated by a fixed-order correction at that scale, but this only enters at next-to-leading-logarithmic order.

One interesting case is when the observable defined at small angular scales $\theta<R_{\rm sub}$ is the weighted energy fraction of all particles with $\kappa=1$. This observable is simply 1 for each of the branches, so the GFFs at the scale $\mu_{\rm sub}$ are

[TABLE]

which are then the input for the fractal observable $\hat{x}_{2}$ for $\theta>R_{\rm sub}$. This effectively removes the sensitivity to nonperturbative physics, allowing us to calculate fractal observables analytically, as long as the scale $\mu_{\rm sub}$ is perturbative. An example of this kind of observable is shown in fig. 20, where the observable is clustered using the recursion relation eq. (25) with $\kappa=1$ for angles $\theta<R_{\rm sub}$ and $\kappa=2$ for $\theta>R_{\rm sub}$. The spike at $x=1$ persists in the numerical evolution, even with very fine bins and a large amount of computing time.171717The generating functional approach (see e.g. ref. Ellis:1991qj ) provides an alternative implementation of the evolution in eq. (3) that can be used to resum (sub)jet radius logarithms Dasgupta:2014yra . This approach may be more amenable to an initial condition with a delta function. This feature is not seen in the Vincia evolution, which at every stage in the parton shower uses a scale closer to $\mu\simeq z\,E_{\rm jet}\,\theta$, where $z$ and $\theta$ are the momentum fraction and opening angle of the splitting. Compared to our choice of $\mu=E_{\rm jet}R$ for the shower as a whole, we would expect the Vincia scale, which corresponds to a larger coupling, to accelerate the depletion of the $\delta$ function in the evolution. It will be interesting to see if this behavior persists with higher-order evolution equations.

An alternative way of viewing the above prescription is that we can build fractal jet observables not just out of hadrons but also out of subjets of radius $R_{\rm sub}$, thus enlarging the range of applicability of the GFF framework. By taking $R_{\rm sub}$ not too small, the observable becomes perturbative. On the other hand, we still want $R_{\rm sub}\ll R$, such that the leading logarithms of $R/R_{\rm sub}$ dominate the observable and eq. (18) gives a reliable description of its behavior.

8 Conclusions

To date, the bulk of analytic jet physics studies are based on either single-hadron fragmentation functions or IRC-safe jet shapes. In this paper, we emphasized the intermediate possibility of IR-safe but collinear-unsafe jet observables defined on a subset of hadrons. We started by introducing the framework of Generalized Fragmentations Functions (GFFs), which are applicable to general collinear-unsafe jet observables. The GFFs are universal functions that absorb collinear singularities order by order in $\alpha_{s}$, which not only restores calculational control, but also implies that the GFFs evolve under a nonlinear version of the DGLAP equations. We then discussed fractal jet observables, defined recursively on an IRC-safe clustering tree with certain initial hadron weights, which satisfy a self-similar RG evolution at leading order given by eq. (18). The higher order evolution is no longer universal, but still self-similar, and has the schematic form in eq. (24).

The simplest fractal jet observables are those with associative recursion relations, whose value does not depend on the choice of clustering tree. This is indeed the case for the weighted energy fractions, studied in sec. 4, which include several observables already in use at colliders, including $p_{T}^{D}$, weighted jet charge, and track fractions. More exotic fractal jet observables depend on the clustering sequence, including the node-product and full-tree observables studied in sec. 5. Remarkably, the structure of the RG evolution for these observables is independent of the clustering tree at leading order.

As one potential application of fractal observables, we studied whether non-associative observables could be useful for quark/gluon discrimination. Indeed, we found examples in sec. 6 which do perform better than the weighted energy fraction $p_{T}^{D}$ currently used by CMS. Though the GFF formalism does not allow us to predict the absolute discrimination power of collinear-unsafe observables, it does allow us to predict the RG evolution of the discrimination power, a feature that will be further exploited in LABEL:HarvardInProgress. To gain more perturbative control, one can work with fractal observables defined on subjets (instead of hadrons), as briefly discussed in sec. 7.

Looking to the future, the next step for fractal jet observables is pushing beyond the leading-order evolution equations. This will require computing the bare GFFs to higher orders in $\alpha_{s}$, as well as extracting GFFs using the matching scheme sketched in eq. (17), and presented in detail at next-to-leading order for $e^{+}e^{-}$ collisions in app. A. More ambitiously, one would like to study correlations between two or more fractal jet observables, which would require multivariate GFFs. Such correlations are known to be important for improved quark/gluon discrimination Gallicchio:2012ez ; Larkoski:2014pca ; Komiske:2016rsd , though even for IRC-safe jet shapes, there are relatively few multivariate studies Larkoski:2013paa ; Larkoski:2014tva ; Procura:2014cba . Together with the work in this paper, higher-order and correlation studies would facilitate a deeper understanding of jet fragmentation, with important consequences for analyses at the LHC and future collider experiments.

Acknowledgements.

We thank Christopher Frye and Andrew Larkoski for helpful comments on quark/gluon discrimination. The work of B.E., J.T., and K.Z. is supported by the DOE under grant contract numbers DE-SC-00012567 and DE-SC-00015476. M.P. is supported by a Marie Curie Intra-European Fellowship of the European Community’s 7th Framework Programme under contract number PIEF-GA-2013-622527. W.W. is supported by the European Research Council under grant ERC-STG-2015-677323, and the D-ITP consortium, a program of the Netherlands Organization for Scientific Research (NWO) that is funded by the Dutch Ministry of Education, Culture and Science (OCW).

Appendix A Generalized Fragmentation in Inclusive Jet Production

In this appendix, we explicitly verify eq. (17) at $\mathcal{O}(\alpha_{s})$. We first calculate the left-hand side of this equation for the measurement of the fractal variable $x$ together with the fraction of the center-of-mass energy carried by the jet, $z_{J}\equiv 2E_{\rm jet}/E_{\rm cm}$. Assuming that $R$ is not so large that all final-state partons get clustered into one jet, we get

[TABLE]

Here, ${i,j}=1,2,3$ and $y_{i}$ is the parton momentum fraction normalized such that $y_{1}+y_{2}+y_{3}=2$. In the following calculations, we identify parton 1 with the quark, 2 with the antiquark, and 3 with the gluon. The angle $\phi_{ij}$ between partons $i$ and $j$ is given by

[TABLE]

and $k$ denotes the parton different from $i$ and $j$. Although the angle $\phi_{ij}$ becomes ambiguous when $y_{i}$ or $y_{j}$ is zero, IR safety ensures that the measurement is not. The term in eq. (A) with $\phi_{ij}<R$ describes the situation where partons $i$ and $j$ are clustered in a jet but parton $k$ is in a separate jet. The final term, where all $\phi_{ij}>R$, corresponds to the situation where all partons are in separate jets. Each of the three partons has a leading-order GFF attached to it. The squared matrix element that enters in eq. (A) is given up to $\mathcal{O}(\alpha_{s})$ by

[TABLE]

where

[TABLE]

Let us now focus on the right-hand side of eq. (17). In our case, the coefficients $C_{i}$ are the standard ones for inclusive fragmentation in $e^{+}e^{-}$ collisions Curci:1980uw ; Floratos:1981hs ; Altarelli:1979kv since the only kinematic variable appearing on the left-hand side of eq. (A) is the jet energy fraction $z_{J}$:

[TABLE]

The coefficients $\mathcal{J}^{(1)}_{q\to qg}$ and $\mathcal{J}^{(1)}_{q\to gq}$ for an $e^{+}e^{-}$ $k_{T}$-like jet algorithm were calculated using the $\overline{\rm MS}$ scheme in ref. Waalewijn:2012sv ,

[TABLE]

while $\mathcal{J}^{(1)}_{q\to q}$ and $\mathcal{J}^{(1)}_{q\to g}$ are given by the finite terms of eq. (2.34) and eq. (2.35) in ref. Kang:2016mcy

[TABLE]

where

[TABLE]

The coefficients for anti-quarks are identical. Note that the relation between $\mathcal{J}^{(1)}_{q\to q}$ and $\mathcal{J}^{(1)}_{q\to g}$ is not simply $z\leftrightarrow 1-z$, because the jet energy $E_{\rm jet}$ rather than the energy of the initiating parton is held fixed. Since $\mathcal{J}^{(1)}_{q\to q}$ and $\mathcal{J}^{(1)}_{q\to qg}$ describe the same splitting in complementary regions of phase space (in-jet versus out-of-jet), their sum vanishes in dimensional regularization,

[TABLE]

The final ingredient we need is the renormalized one-loop expression for the GFF (see eq. (2.3)),

[TABLE]

Let us first verify the cancellation of IR divergences between left- and right-hand sides in eq. (17). On the latter, these solely come from $C_{q}^{(0)}(z_{J},E_{\rm cm},\mu)[\mathcal{F}_{q}^{(1)}(x,\mu)+\mathcal{F}_{\bar{q}}^{(1)}(x,\mu)]$. On the left-hand side, we find

[TABLE]

which demonstrate the cancellation of the IR divergences. Note that the term on the first line of eq. (A) proportional to $\mathcal{F}_{i}^{(0)}$ does not contribute here because it is $y_{2}$-independent and

[TABLE]

To verify that also the finite terms match in eq. (17), we expand the angular constraint in the small $R$ limit as

[TABLE]

which implies $y_{k}\approx 1$ and $y_{j}\approx 1-y_{i}$. We first consider the $\theta(R-\phi_{13})$ term in eq. (A), which gives

[TABLE]

As the integral over $y_{2}$ yields a $\ln(1-y_{1})$, the resulting $\ln(1-y_{1})/(1-y)_{+}$ is not properly regularized, leaving the coefficient of $\delta(1-z)$ undetermined. As we will see, however, this ambiguity cancels exactly against the one arising from $\mathcal{J}_{q\to q}^{(1)}$, due to eq. (68). The $\theta(R-\phi_{23})$ term gives the corresponding contribution with quark and anti-quark interchanged, whereas the $\theta(R-\phi_{12})$ term is $\mathcal{O}(R^{2})$ suppressed due to the $e^{+}e^{-}\to q\bar{q}g$ squared matrix element.

For the last contribution in eq. (A), we rewrite

[TABLE]

where the first term in the sum corresponds to the calculation of the matching coefficients for inclusive fragmentation, thus yielding the $C_{i}(z_{J},E_{\rm cm},\mu)\mathcal{F}_{i}(x,\mu)$ contribution on the right-hand side of eq. (17). For the remaining terms, we can follow the same strategy as in eq. (A). For example, the $-\theta(R-\phi_{13})$ term gives

[TABLE]

The similarity with the calculation in eq. (A) and the relationship between $\mathcal{J}^{(1)}_{q\to q}$, $\mathcal{J}^{(1)}_{q\to g}$ in eq. (A) together with eq. (68) make this straightforward to verify. The $(\dots)_{13}$ term cancels in the sum with eq. (A). The $-\theta(R-\phi_{23})$ term corresponds to the term with quark and anti-quark interchanged and the $-\theta(R-\phi_{12})$ contribution is again suppressed by $\mathcal{O}(R^{2})$. This completes the check of eq. (17) at $\mathcal{O}(\alpha_{s})$.

Appendix B A Non-Fractal Example: Sums of Weighted Energy Fractions

While eq. (1) is rather general, there are of course many collinear-unsafe observables that are not fractal jet observables. In this appendix, we give an explicit example of an observable that does not satisfy the requirements in sec. 2.3.

Consider two weighted energy fractions

[TABLE]

for particle weights $w_{i}$ and $v_{i}$, and energy exponents $\kappa$ and $\lambda$. Individually, $x$ and $y$ are described by the evolution equation in eq. (18). On the other hand, their sum

[TABLE]

is not a fractal jet observable, though it still can be described by a GFF.

To see this, consider the GFF for $t$, $\mathcal{F}_{i}(t)$, which can be written in terms of a joint GFF for $x$ and $y$ as

[TABLE]

The evolution equation for the joint GFF follows from the analysis in eq. (2.3), leading to

[TABLE]

Plugging eq. (79) into eq. (78), we can insert a factor of

[TABLE]

to perform the integrals over $y_{1}$ and $y_{2}$. The resulting equation is

[TABLE]

As written, this is a valid GFF evolution equation, but the GFF for $t$ explicitly involves the joint GFF for $x$ and $y$, so we do not get an evolution equation of the form of eq. (18).

If and only if $\kappa=\lambda$, can we cancel the $x_{1}$ and $x_{2}$ terms inside of the $\delta$ function in eq. (81). In that case, we can rewrite the joint probabilities as probability densities for the sums $t_{1}=x_{1}+y_{1}$ and $t_{2}=x_{2}+y_{2}$, so that the evolution equation is of the desired fractal form. Of course, $\kappa=\lambda$ just corresponds to a regular weighted energy fraction with weights $w_{i}+v_{i}$, so this is not a new fractal observable.

Appendix C Software Implementation

The software to perform the RG evolution in this paper is available from the authors upon request. In this paper, we discuss some of the specifics of its implementation. A public version of the code is planned for a release some time in the future.

C.1 Running Coupling

Because we only perform leading-order evolution, the running of $\alpha_{s}$ is strictly speaking only required at leading-logarithmic accuracy. In our implementation, though, the running of the strong coupling is included using the $\beta$ function at $\mathcal{O}(\alpha_{s}^{3})$,

[TABLE]

[TABLE]

The running coupling at the scale $\mu$ is given by solving eq. (82) iteratively to order $\mathcal{O}(\alpha_{s}^{3})$,

[TABLE]

where $L=\ln\frac{\mu^{2}}{\Lambda_{\rm QCD}^{2}}$. Using the PDG value $\alpha_{s}(M_{Z})=0.1181$ gives the boundary condition $\Lambda_{\rm QCD}=0.2275$ GeV. The group theory factors for QCD are $C_{F}=\frac{4}{3}$, $T_{F}=\frac{1}{2}$, and $C_{A}=3$. For applications to the LHC running at 13 TeV, the number of quark flavors is $n_{f}=5$.

C.2 Discretization

The evolution equation in eq. (18) can be solved by binning the values of the GFFs in the $x$ variable. If the GFF domain is partitioned into $N$ bins, eq. (18) becomes a set of $(2n_{f}+1)N$ coupled ordinary differential equations. The evolution equation for the binned GFF for bin $n$, $\widetilde{\mathcal{F}}_{i}(n,\mu)$, is given by181818This equation is written for $N$ equal-width bins for simplicity of notation. The generalization to unequal bins is straightforward, and the software implementation is set up to handle variable bin widths if desired.

[TABLE]

where $x_{n_{1}}$ and $x_{n_{2}}$ are the positions of the midpoints of the $n_{1}$-th and $n_{2}$-th bins. Note that eq. (85) is written in terms of $\ln\mu$ instead of $\mu$, since this is how the evolution was implemented numerically to make the step size and numerical errors more consistent. In principle, the $\delta$ function could be used to carry out the $z$ integral exactly. In practice, it is easier to discretize the $z$ integral and use the $\delta$ function to choose the $x$-bin corresponding to each triplet $(z,x_{1},x_{2})$. This is because inverting $\hat{x}$ to solve for $z$ analytically for general $x_{1}$ and $x_{2}$ is not possible. Doing so in advance separately for each value of $x$, $x_{1}$ and $x_{2}$ can be prohibitively memory intensive for large numbers of bins.

The splitting functions are approximated by the analytic value of their integral over the width of the bin. For our analysis, we need the following splitting functions:

[TABLE]

where $P_{q\rightarrow gq}(z)$ is the splitting function for a quark radiating a gluon with momentum fraction $z$, the integration constant for integrals of the plus distributions are fixed by

[TABLE]

and $\beta_{0}$ is given in eq. (83).191919The $1/z_{+}$ and $\delta(z)$ terms in $P_{q\rightarrow gq}(z)$ and $P_{g\rightarrow gg}(z)$ are necessary because the evolution in eq. (3) requires distributions that are also regulated at $z=0$. When performing the integration, terms with a plus-function regulator must be handled correctly for the endpoint bins. If the regulated functions have the following primitives

[TABLE]

then their integrals over the $n$-th bin are implemented by

[TABLE]

In our implementation, the integration range $z\in[0,1]$ is divided into $n_{\text{rough}}$ bins, and the first and last bin are then further subdivided by a factor of $n_{\text{fine}}$. The user can specify these two parameters. For the results presented in this paper, the values used were $n_{\text{rough}}=1000$ and $n_{\text{fine}}=100$. The finer division of the endpoint bins is necessary to accurately capture the singular behavior of the splitting functions near $z=0$ and $z=1$. For many GFFs, this is not necessary, but consider the weighted energy fractions, whose recursion relation satisfies

[TABLE]

For $\kappa<1$, there are poles in the derivative of $\hat{x}$ at $z=0$ and $z=1$, resulting in a noticeable dependence on $n_{\text{fine}}$. This is shown in fig. 21 for the case of $\kappa=0.5$, with all particle weights one. Once we increase $n_{\rm fine}=100\rightarrow 1000$, the maximum change in the value of the evolved GFFs in a single $x$-bin is less than 0.06$\%$.

C.3 Runge-Kutta Algorithm

After the discretization in eq. (85), the RG evolution is performed with an embedded fifth-order Runge-Kutta method adapted from LABEL:numrecipes. This method requires six evaluations of the right side of eq. (18), which on the $k$th step can be combined to give a fifth-order estimate $y_{k+1}$ of the desired function after a step of size $h_{k}$. These computations can be recombined with different coefficients to give a fourth-order Runge-Kutta estimate $y_{k+1}^{*}$. The difference between these two methods then gives an estimate of the local truncation error. The error estimated this way applies to the fourth-order value $y_{k+1}^{*}$, but we take the (more accurate) fifth-order value. This ensures that our solution is actually slightly more accurate than our error indicates. Estimating the error on this fifth-order solution would require calculating a still-higher order step.

Once a step $h_{k}$ is taken, with an error $\mathcal{E}_{k}$, we would like to choose an appropriate trial value for our next step. This fourth-order error estimate scales as $\mathcal{O}(h^{5})$, so we choose the next step, $h_{k+1}$, to be

[TABLE]

Here, $\mathcal{E}_{k+1}$ is the projected error in the $(k+1)$th step, and $S$ is a safety factor taken to be $0.9$. This formula allows the step size to grow if the error is much smaller than our tolerance. If the error is larger than the tolerance, the step fails, and is retried with a smaller step.

It is important that the algorithm be able to dynamically change step size in order to evolve a solution efficiently while keeping errors within desired limits. At low scales, the strong coupling grows large, and the solution changes rapidly. Numerical precision therefore requires small step sizes in this region. At high scales, asymptotic freedom ensures that the solutions change slowly, so much larger step sizes result in the same level of accuracy. This procedure requires a prescription for the maximal acceptable error. For a system of $M\equiv(2n_{f}+1)n$ coupled ODEs, there is a separate $\mathcal{E}_{k}^{m}$ for each $m\in M$. The step is considered a failure unless every equation is within its error tolerance. The error $\mathcal{E}^{m}_{k}$ for the $m$th equation on the $k$th step is required to satisfy

[TABLE]

The value $\epsilon$ is an overall upper limit which was set to $10^{-9}$ for the GFF evolution. The last numerical term in the denominator is required to avoid artificially large errors when the domain of the GFFs input into the program exceeds the actual support of the GFF. As an additional constraint, our algorithm sets a maximum step size of ${\rm d}\ln\mu\leq 0.4$. Note that the same step size is used for every equation in the system.

Appendix D Numerical Stability

All of the RG results in this paper are based on the numerical solution of eq. (18) for upwards evolution in the scale $\mu$. The reason is because downward evolution is numerically unstable, in the sense that small irregularities in the initial conditions amplify into large fluctuations, especially for the gluon GFFs. This behavior is illustrated in fig. 22, where gluon and quark-singlet GFFs are evolved downward from 4 TeV to 100 GeV.

Heuristically, if evolution upwards in scale is analogous to convolution of the GFFs, evolution downwards is akin to deconvolution, a problem known to be ill-posed. To verify that the instability is inherent to the differential equation, and not merely a numerical artifact, we checked that the envelope shown in fig. 22 is not affected by choosing a smaller step size or more stringent error bound in the Runge-Kutta algorithm. To get a sensible result, one could use a numerical regularization method such as Tikhonov regularization tikhonov1977solutions , though we do not do so here. Note that in general, if the evolution in one direction is stable, such that small fluctuations get washed out, the evolution is expected to be unstable in the reverse direction.

Appendix E Moment Space Details

In this appendix, we give details of the moment space analysis from sec. 4.5, as well as perform similar analyses for the non-associative observables from sec. 5. The moments of the GFFs are defined by

[TABLE]

where the zeroth moment is just the normalization,

[TABLE]

This convention follows the standard nomenclature of probability theory. Applying $\int_{-\infty}^{+\infty}\text{d}x\,x^{N}$ to both sides of the evolution equation in eq. (18) gives the moment space evolution equation,

[TABLE]

In order to proceed further, we need the specific form of the recursion relation, $\hat{x}$. We now discuss the details for each of the sets of observables studied in this paper.

E.1 Weighted Energy Fractions

Inserting the weighted energy fraction recursion relation eq. (25) into eq. (95) leads to

[TABLE]

assuming that $N$ is integer and using the binomial theorem. As in eq. (36), the moments of the splitting functions are defined as

[TABLE]

with the convention that $\overline{P}_{i\rightarrow j,k}(N)\equiv\overline{P}_{i\rightarrow j,k}(N,0)$. For any real $N>0$, they can be expressed in terms of the digamma function $\psi_{0}(N)$ and the Euler-Mascheroni constant $\gamma_{E}$,

[TABLE]

Alternatively, one can use the harmonic number function, $H_{N}=\gamma_{E}+\psi_{0}(N+1)$. These expressions for all positive real numbers are necessary to evaluate the moment space evolution equation in eq. (37) for non-integer $\kappa$. Note that $N$ is shifted up by one from the expression usually seen in the literature, because our convention for moments in eq. (93) is shifted by one as well compared to Mellin moments.

E.2 Node Products

We now insert the recursion relation for the node products from eq. (51) into eq. (95). This leads to evolution equations with additional terms compared to those for the weighted energy fractions. These terms have splitting kernels of the form

[TABLE]

for $a>0$ and $b,c\geq 0$. These integrals are convergent, so no plus function regulators are required. They can also be performed analytically for general $a$, $b$, and $c$. Explicitly, the first moments of the quark-singlet and gluon GFFs evolve as

[TABLE]

The additional constant terms are defined as

[TABLE]

which can be evaluated in terms of $\Gamma$ functions. The additional terms drop out of the equation for the first moments of the non-singlet GFFs, so these still evolve according to eq. (38). The third term in eq. (51) leads to several more terms in the evolution equations for higher moments.

In fig. 23, we plot the $\mu$ evolution of the gluon and quark-singlet GFF moments for node products with $\kappa=\{1,4\}$ and $p=\{-1,0,1\}$. The first and second moments were computed at the scale $\mu=100$ GeV from the GFFs in fig. 10, averaged over the different parton showers and $R$ values (as described in sec. 4.2). These average moments were evolved to the scale $\mu=10^{7}$ GeV using eq. (108) and the corresponding second moment equation. For comparison, the first and second moments of the GFFs extracted from the parton shower average at the scale $\mu=$ 4 TeV are shown as dots and diamonds, respectively.

E.3 Full-Tree Observables

For full-tree observables with recursion relation given in eq. (54), the moment space evolution equations are of the same general form as for the weighted energy fractions,

[TABLE]

but with different splitting kernels,

[TABLE]

To our knowledge, these integrals do not have a closed form solution for general values of the parameters $\kappa$ and $\xi$, but it is straightforward to evaluate them numerically. If $M=0$ or $M=N$, these integrals are sensitive to the plus-prescription in the splitting functions. Explicitly, for the first moment in the quark-singlet basis,

[TABLE]

In fig. 24, we show the evolution of the first two moments of the GFFs for $\kappa=2$, $\xi=\{-2,2\}$, and $p=\{-1,0,1\}$. In this case, the evolution agrees well with the value extracted from the parton shower average at $\mu=4$ TeV.