UvA-DARE (Digital Academic Repository) One-loop jet functions by geometric subtraction

,

If these restrictions are tight, they lead to large logarithms in the corresponding cross section.For example, for Higgs plus one jet production with a veto on additional jets with transverse momentum above p veto T , the cross section takes the following form where σ 0 is the leading-order cross section, and the coefficients c n,m are independent of p veto T .For a tight veto p veto T m H ∼ p jet T , the expansion in α s deteriorates due to the large logarithms and resummation is crucial to improve convergence and reduce the theory uncertainty.Resummation captures the dominant effect of higher-order corrections, effectively treating ln(m H /p veto T ) ∼ 1/α s .Large logarithms arise because the cross section involves multiple scales that are widely separated.Resummation of these logarithms can be achieved by factorizing the cross section into components that each involve a single scale, using diagrammatic methods in QCD, see e.g.[1][2][3][4][5][6][7][8], or Soft-Collinear Effective Theory (SCET) [9][10][11][12][13].For exclusive Higgs plus one jet production, discussed in eq.(1.1), this takes on the following (schematic) form σ ∼ HSBBJ [14,15].The hard function H describes hard scattering, the soft function S encodes the effect of soft radiation, and the beam functions B and jet function J account for initial-and final-state collinear radiation.The structure of this factorization not only depends on the process, but also on the observable and can involve convolutions between ingredients (though it is simply a product in the above example).Because each ingredient in the factorization involves a single scale, the large logarithms can be resummed by evaluating each ingredient at its natural scale and using the renormalization group to evolve them to a common scale.Alternatively, an automated approach to resummation was pursued in refs.[16,17].
In this paper we focus on calculating one-loop jet functions, which enter in resummed cross sections starting at next-to-leading logarithmic (NLL ) accuracy.Resummation at NLL includes the two-loop cusp anomalous dimension and one-loop (non-cusp) anomalous dimensions.Jet functions have been calculated for a wide range of observables, including the invariant mass [18][19][20][21][22][23], the family of e + e − event shapes called angularities with respect to the thrust axis [24][25][26] or Winner-Take-All axis [27,28], Sterman-Weinberg jets [29,30], the cone and the k T family of jet algorithms for exclusive [30,31] and inclusive [32,33] jet production.Jet functions have also been considered for a range of jet substructure observables, such as the jet shape [34][35][36].In our calculations we treat quarks as massless and restrict to infrared-safe observables.An example of a massive quark (initiated) jet function is given in refs.[37,38], and an example of an infrared-unsafe jet observable is the electric charge of the jet [39,40].
We briefly comment on the other ingredients in the factorization: A general approach to calculating soft functions has been developed in refs.[41][42][43][44].In particular, the Soft-SERVE package [44] provides two-loop soft functions for processes with two collinear directions (i.e. two jets in e + e − or 0 jets in pp collisions), and an extension to N jets is in progress [45].Hard functions can be obtained from the IR finite part of helicity amplitudes, as long as the color of the initial (final) particles is not averaged (summed) over, see e.g.ref. [46].
The difficulty in calculating jet functions lies in the phase-space integration, which depends on the observable.When feasible, an analytic approach is superior.However, there are observables for which even the one-loop jet function is highly nontrivial, such as jet broadening [25] and the jet shape [36], for which fully analytic results are difficult to obtain or have not been obtained yet.The numerical approach we develop here offers a promising alternative, addressing the collinear and soft divergences in a general way, thereby automating the calculation of one-loop jet functions for a broad range of observables.At minimum, our work provides a valuable cross check for analytic calculations.
The poles in the dimensional regulator are obtained analytically, possibly up to an integral over the azimuthal angle, and depend on the collinear and soft behavior of the observable.This soft behavior is described by a power law, and therefore simply characterized by the exponent and coefficient.Extracting these parameters may require solving non-trivial algebraic equations, and we develop a procedure to simplify this step.The full details/complications of the measurement only enter in the finite term, which can be integrated numerically.We have implemented our approach in a Mathematica package, Geometric One-loop Jet functions (GOJet), which accompanies this paper.GOJet can handle a large class of infrared-safe observables, including all the observables listed above.
Using GOJet we provide explicit examples of the method for the angularities with respect to the Winner-Take-All axis, the cone and k T -clustering jet algorithms and the jet shape.Furthermore we calculate for the first time the one-loop jet function for angularities with respect to the thrust axis including recoil.We cross check our result against existing results in the literature for the specific case of jet broadening [25] and for the case of no recoil [24,47].
The remainder of the paper is structured as follows: In section 2 we discuss how we use geometric subtraction to calculate jet functions, including a simple example.The GOJet package, which provides a Mathematica implementation, is discussed in section 3.In section 4, we use our package to calculate several one-loop jet functions, and we conclude in section 5.

General Method
In section 2.1 we will discuss geometric subtraction and how we apply it to calculate oneloop jet functions.Technical aspects related to the treatment of Heaviside theta functions in our calculation and infrared safety are discussed in sections 2.2 and 2.3, respectively.We illustrate our method by calculating the jet function for the e + e − angularity event shapes in section 2.4, with further examples in section 4.

Subtraction scheme
The jet function depends on the flavor i = q, g of the initiating parton and the jet observable, and has a perturbative expansion in α s At tree level the jet consists of a single quark or gluon, and in general J (0) i = 1 in the appropriate units. 1 The one-loop contribution is given by the collinear limit of two final-state partons Here s denotes the invariant mass of the two partons, and z and 1 − z the momentum fractions of the partons.The squared matrix element is contained in Q i (s, z, φ), with P i (z) the (sum of) splitting function(s).The calculation is performed in d = 4 − 2 dimensions and the MS-renormalization scheme with renormalization scale µ is employed.
For certain observables an additional rapidity regulator η and corresponding rapidity scale ν are required [48][49][50][51][52][53], which is included in eq.(2.2) for generality.This arises when the collinear and soft functions have the same invariant mass scale µ, with transverse momentum measurements being the typical example.For the extension of eq. ( 2.2) to a two-loop example, see ref. [54].
The measurement in a jet function can often be written as To avoid distributions, we require the user to rewrite the measurement as a Heaviside theta function by integrating, i.e.Θ[O − f (s, z, φ)], where we are now cumulative in O. 2 We therefore assume that the measurement M obs (s, z, φ) is a Heaviside theta function, which cuts out a certain region of the collinear phase space, as illustrated in figure 1 (suppressing φ dependence).An advantage of cumulative distributions is that they involve logarithms rather than plus distributions: In section 2.2, a technical point related to rewriting measurement delta functions in terms of theta functions will be discussed.There are also measurements that are naturally theta functions.For example, the k T -family of jet algorithms requires both particles to be clustered into a jet with radius parameter R, , where p T is the transverse momentum of the jet.In principle these phase-space constraints M obs can depend on the azimuthal angle φ as well, but since there is no singularity associated with the φ integration, we will only include φ when needed.
The jet function in eq. ( 2.2) has divergences as s → 0 (collinear divergence), and z → 0 and z → 1 (soft divergences), which occur at the phase-space boundaries in figure 1. Infrared-safe observables must always either include or exclude the entire collinear divergence (the red line in figure 1), as will be discussed more in section 2. view of collinear subtraction, one can consider the jet function (as long as it contains the collinear divergence) as a collinear counterterm.Different observables can then be viewed as different schemes, differing in the extent that soft and soft-collinear divergences are included in the observable.For instance, region 1 of the general observable illustrated in figure 1 only contains the collinear and part of the soft-collinear singularities.By contrast, region 2 only contains part of the soft and none of the collinear divergence.Region 3 does not contain any soft or collinear divergent parts of phase space and does therefore not have to be regulated.Another possibility would be to consider an observable which corresponds to the complement of region 1, which naively causes problems because it develops a logarithmic singularity for s → ∞.However, its one-loop jet function is given by minus the jet function for region 1, because the integral over the full collinear phase space results in a scaleless integral.
To define a general subtraction scheme for calculating jet functions for infrared-safe observables, we follow the approach of geometric subtraction [55].We would like to define a finite part of the jet function as follows: where we introduced the dimensionless slicing parameters A and B, that remove the soft and collinear divergence, and which we subsequently want to take to zero.The central idea of geometric subtraction rests on the identity: where we exploited that a is small on the second line to replace f (x) by f (0) in the second term.However, the expression on the second line is now regulated for any 0 < a ≤ 1, leading to a duality between slicing and subtraction schemes.To obtain the full jet function from the above finite part, counterterms need to be added to reinstate the part of the integral that is removed by the cuts.The counterterms generated in this way are added back in integrated form, regulated dimensionally and if needed also with a rapidity regulator, and may give a finite contribution to the jet function.While a subtlety arises in general when different limits do not commute, here we do not face this problem as the collinear and soft singularities are factorized.For the small A limit in eq.(2.4) we can then straightforwardly apply eq.(2.5).However for the parameter B nothing is gained from this procedure, because the jet function is already in the limit of small s and the counterterm generated is the original integral itself.
To obtain a simpler counterterm in the s < B region, we can however use a simpler observable, which we choose to be the jet mass, as a collinear counterterm.(This was also used in the geometric subtraction scheme [55].)Since the region of the s-z plane corresponding to the jet mass is box-shaped, we will refer to this collinear counterterm as the box.A subtlety now appears due to the difference of soft and soft-collinear divergences included in the box counterterm and the given observable M obs , which as discussed above may not be the same.To deal with this problem we introduce separate soft counterterms for both the box counterterm and the M obs term in the region s < Bµ 2 , as discussed in detail below.
These considerations lead us to the following final decomposition of the jet-function into finite and divergent parts: where the arguments s, z, φ are suppressed and A, B are positive real numbers with A ≤ 1.The first term in G i,obs,1 corresponds to the finite part defined in eq.(2.4), and the other terms correspond to integrated counterterms.It is straightforward to check that the sum of G 1 , G 2 and G 3 is equal to the original one-loop jet function.
The advantage of the above decomposition is that G 3 is observable independent, G 2 only depends on the soft limit of the observable (which can be encoded by a few parameters at one-loop order, see eq. (2.8)) and G 1 is finite.In eq.(2.6), Q 0 and Q 1 denote the soft z → 0 and z → 1 limit of Q. Explicitly,  A graphical representation of our subtraction scheme in eq.(2.6).We have only included the soft counterterms for z → 1 for legibility.Shown are the restriction on the measurement from the observable (blue line), the soft limit of the observable (red line), the box (green line), the cut on z arising from A (pink line).Blue plus (minus) areas correspond to positive (negative) contributions of the full integrand Q i M obs , while red plus (minus) areas correspond to positive (negative) contributions of Similarly, M obs,0 and M obs,1 denote the soft z → 0 and z → 1 limit of the measurement M obs .The soft limit can contain multiple boundary conditions on the phase space, which we account for by writing M obs,0 and M obs,1 as a sum of Heaviside theta functions that constrain the integration over s as a function of z.Moreover, they will follow a power-law behavior parametrized by where the sum on r is over different regions (see figure 1), and the parameters c i , α i depend on the observable, and can depend on φ as well. 3We also allow for a constraint Φ on the azimuthal angle, as will be discussed in section 2.3.Depending on the observable, each soft boundary condition will therefore follow one out of three distinct behaviors shown in figure 1: the upper boundary of R 1 corresponds to α < 0, the lower boundary of R 2 to α = 0, the upper boundary to α > 0 and R 3 does not extend into the soft region.Finding c 0,1 and α 0,1 can be nontrivial, and we will discuss a strategy to do so for an involved example in section 4.2.
We will now discuss the decomposition in eq.(2.6) in more detail, using the graphical representation in figure 2 for the k T algorithm.In order to get a finite G 1 in figure 2a, we subtracted the collinear singularity and the soft singularities.The collinear singularity is removed by the box, replacing M obs by M obs − 1 when s ≤ Bµ 2 , such that M obs (s = 0, z, φ) − 1 = 0.The soft singularities get accounted for by subtracting the z → 0 and/or z → 1 limits of the integrand.Indeed, one can see that in figure 2a the blue plusses and red minuses cancel as z → 1.The resulting integral G 1 is now finite.For general observables, G 1 in eq.(2.6) may be hard to calculate analytically, and one has to resort to numerical integration techniques.In the examples in section section 4, we will use the Cuba implementation of Vegas [56] to perform the integrations.Convergence problems in the numerical integration may arise due to the mismatch of the observable and its soft approximation, which generally can lead to integrable singularities.If these problems are severe it can help to find an explicit remapping of the counterterm, which decreases the mismatch between the observable and its soft limit.We present a method for how this can be achieved with a worked through example in app.B.
Let us now discuss the integrated counterterms.Due to their simplicity, the counterterms can be calculated analytically, which we discuss for a single region r in the sum in eq.(2.8).Let us first focus on the soft counterterms, which are contained in G 2 shown in figure 2b.The soft limits of the integrand Q i M obs are given by Q i,0 M obs,0 and Q i,1 M obs,1 , see eqs.(2.7) and (2.8).The constants c i and α i are user input in our code, see section 3.For values α = 1, no rapidity regulator is needed and η can be set to 0, leading to the following soft counterterm . (2.9) For α = 1 one needs a rapidity regulator and the corresponding expression is given in app. A. The box counterterm G 3 in figure 2c is given by The integral over φ has been carried out for Θ(Φ) = Θ(φ + − φ)Θ(φ − φ − ) leading to the function The evaluation of this integral and its expansion to order 2 is presented in app.C. The chosen subtraction bears fruit in the simplicity of the integrated counterterms.The corresponding Laurent series in can be expressed solely in terms of the Riemann zeta function at integer values, given that only pure Gamma functions appear.From an analytic point of view, the potentially more complicated pieces are instead captured in the finite part, which depends on the details of the observable and can be calculated numerically to arbitrary high order in .Notice that the soft counterterm G i,2 can give rise to more complicated integrals if the coefficients c ± i depend on the azimuthal angle φ.One may be able to carry out this integral analytically in certain cases, but this can certainly not be done in general.This is not a problem, because one can expand in and η before integrating over φ.

Delta and theta functions
In our subtraction scheme we assume that the observables restrict the integration to certain regions of phase space via Heaviside theta functions.However, many observables O are naturally expressed in terms of Dirac delta functions, requiring one to rewrite it using where f is a function of the kinematics of the collinear splitting, and possibly external parameters.The sign ± should be chosen such that the theta function does not vanish at tree-level, which ensures that the poles are included in the one-loop jet function.For example, if O ≥ 0 and at tree-level O = 0, one needs to choose the plus sign in eq.(2.12).
In perturbative QCD one often works with the following convention for the Dirac delta function, This differs from the definition given in standard math literature where the lower boundary b must be strictly less than zero.If the delta function that encodes the measurement satisfies eq.(2.13), this has implications for the definition of the Heaviside function on the right-hand side of eq.(2.12).In particular, one must demand then that Θ(0) = 0. To see this, consider a function g(x) with 0 ≤ x ≤ 1.From we conclude that Θ(0) = 0.While this is not of much concern when a theta function is integrated over, there are situations where it must be taken into account.As an example, the jet shape calculation involves a jet function describing the energy fraction z inside a cone, see section 4.3.Switching to a cumulant variable for z, we need to choose δ(z −. . . ) = −d/dz[−(z − . . .)], because 0 ≤ z ≤ 1 and z = 1 at tree-level.If we now want to calculate the average momentum fraction from the cumulant tree-level result we have to take θ(0) = 0 to find agreement with the direct evaluation using the delta function.

Infrared safety and limitations on the observable
While so far our discussion was mostly based on the s-z plane, there are observables which depend also on the azimuthal angle φ.The integration domain is then parametrized by coordinates (s, z, φ) and IR safety requires the full s = 0 plane to be included or excluded by the observable, i.e. the set of points However, our method allows for a special class of IR-unsafe observables, where only subdomains of the collinear plane with the azimuthal angle bounded between constant values are included/excluded by the observable, i.e.
with 0 ≤ φ − < φ + ≤ π.This is illustrated in figure 3a.An IR-unsafe observable which is not of this form, and currently not supported by GOJet, is illustrated in figure 3b.Here φ ± vary as functions of z across the collinear plane in such a way that not the full s = 0 plane is included in the integration domain.For s > 0 the bounds on φ can depend on z.
GOJet can also handle IR-unsafe observables that include just z = 0 and/or z = 1 of the s = 0 plane, which only require soft counterterms.

Example: Angularities with the Winner-Take-All axis
We will now illustrate our scheme by considering the family of e + e − event shapes called angularities [57] parametrized by b 4 .Here Q is the center-of-mass energy, and the sum runs over all particles i in the final state with energy E i and angle θ i with respect to some axis.The final expression is only valid in the small-angle limit, which is appropriate for the jet function calculation, highlighting that e b probes the angular distribution with exponent 1 + b > 0.
While angles were originally taken with respect to the thrust axis, we will here use the Winner-Take-All axis [27].For the one-loop jet function this axis is simply along the most energetic particle in the jet, so the only non-zero contribution in the sum on i in eq.(2.19) comes from the least energetic particle, with θ i the angle between the two partons in the jet.Noting that s = 2p For angularity exponent b < 1, the observable is unbounded from above, similar to the top curve of region 2 in figure 1.In the notation of eq.(2.8), we see that the soft limit of the observable is characterized by The one-loop contribution to the jet function is obtained by plugging in these these constants in eqs.(2.9) and (2.10) to calculate G 2 , performing the integration over s and z for G 1 , and adding these contributions to the box G 3 .Performing the integration over s analytically and the integration over z numerically for b = 2, we obtain q,e 2 = where we used µ = Q(e c 2 ) 1/3 to calculate the constant contribution and reinstated the logarithmic behaviour afterwards.Our result agrees with the expression in refs.[27,28] up to order 10 −11 . 5For b = 0 the rapidity regulator is required.In that case we find q,e 0 = in agreement with ref. [27].

GOJet Program
The GOJet Mathematica-package automatically performs the subtraction, given the observable and its soft limit (see eq. (2.8)) as input.One can either let Mathematica perform the numerical integration or choose to export the integrand.The latter feature may be useful if NIntegrate either has difficulty converging or is not fast enough.In such cases it can be advantageous to use algorithms such as Vegas, that are faster due to their implementation in C++ or Fortran.A general overview of the various functions included in the package is given in section 3.1.A detailed description of their input is given in section 3.2, with a worked-out example in section 3.3.

Functions
There are a total of 12 different functions, listed in section 3.2, which the user can access.As indicated by their names half of these are for calculating gluon jet functions while the other half are for calculating quark jet functions.Restricting to the former, PolesGluon returns the pole terms in and η for the gluon jet function and GluonJet returns the integrand of the finite terms, by which we here refer to the 0 η 0 -term.In addition, GluonJetN performs the numerical integration over the cube 0 ≤ s, z, φ ≤ 1 of this integrand.This integration domain is the result of mapping s → s/(1 − s) and φ → πφ, which also stabilizes the integration over s.Note that GluonJet also contains the 0 η 0 -pieces of the counterterms G 2 and G 3 , which are already integrated over analytically.For the convenience of the user these pieces are simply added in integrated form since they are not altered by the trivial numerical integration over the unit cube.
Let us now discuss the arguments of the functions in general terms.The first arguments encode the measurement O and its soft limit O 0 and O 1 corresponding to the limits z → 0 and z → 1, respectively.The observable should generally be IR safe, with some exceptions discussed in section 2.3.Furthermore, we require certain restrictions on the form of the soft limits.Specifically, it is not possible to restrict the φ-integration boundaries via O 0 and O 1 , whose format is fixed.It is however possible to apply s, z-independent constraints on the boundaries of the φ-integration through the separate argument Φ, which are the same for the finite part as well as the counterterms.
The next set of arguments specify the regularization and IR scheme: the need of a rapidity regulator or collinear regulator is controlled by the switches rr and box, respectively.The explicit cut for the soft limits and box is specified by A and B (see eq. (2.6)).The independence of the final result on these parameters provides a useful cross-check for the calculation.A specific choice of these parameters can also be used to improve the convergence of the numerical integration.For the gluon jet function, the number of quark flavors is specified through the argument nf.The number of colors has been fixed to three, but the full dependence on the Casimirs can be easily reconstructed from the answer.The final set of arguments enables the user to specify the integration method or output format for the integrand.
Finally, we also allow for more complicated observables, where the phase-space restriction due to the measurement breaks up into more than one region.The corresponding functions have "Regions" appended to their name, and contain additional arguments specifying possible dependence on external parameters in the regions.

Input format
Here we specify the syntax of each of the functions: The variables used to describe the input are: • O: The list of argument(s) of the Heaviside theta function encoding the bounds imposed by the measurement.More specifically, O contains the arguments of the Heaviside theta functions M obs in eq.(2.2).For the case of a single region, the elements of the list correspond to the arguments of Heaviside theta functions, whose product constrain the region.In the case of multiple regions, O is a list of lists.The entries of the outer list correspond to the different regions, each entry is again a list of constraints containing the arguments of the Heaviside theta functions M r obs constraining the particular region.This allows the user to implement arbitrary sums of products of Heaviside theta functions.
• R 1 (R 0 ): List of lists which contain arguments of Heaviside theta functions which depend only on external parameters for each region in the limit z → 1 (z → 0).The length of this list is therefore equal to the number of soft regions that emerge in the soft limit.Regions that do not depend on external parameters need {1} as input in their respective position in the list.The number of soft regions can be less than the number of regions, but should match with the lists for O 0 and O 1 below.In particular, regions may merge or disappear in the soft limit.R 1 (R 0 ) can also be used in cases with just one region where there is dependence on external parameters in the soft limits.
describing the lower and upper boundary of the region in the limit where z → 1 (and similarly for z → 0), see eq. (2.8).If there is no lower boundary, c − 1 is just 0. When considering multiple regions, O 1 (O 0 ) is a list of lists where each region has an upper and a lower boundary of the aforementioned format.
• Φ: List of arguments of the Heaviside theta functions that impose constraints on the azimuthal angle φ, i.e., the input {φ + − φ, φ − φ − } will constrain φ − < φ < φ + .In the case of multiple regions that contain collinear and/or soft divergences we require the range on φ to be the same for all regions.(Arbitrary constraints on φ can of course be encoded in O; but these are not allowed to survive singular limits; that is the they should match the boundaries imposed by Φ in these limits; see section 2.3 for more details.) • rr: Boolean variable specifying whether a rapidity regulator should be included, which we implemented as This corresponds to the more conventional factor (ν/((1 − z)z ω)) η , for the scale choice ν = 1 2 ω.The user can always reconstruct the full dependence on the scale ν a posteriori, given the knowledge of the 1/η pole.
• box: Boolean controlling whether a box is needed to handle the collinear divergence.
It should be included when the region of phase space includes s = 0 and not otherwise (in line with the restrictions outlined in section 2.3).
• A: Real number specifying the region where the soft counterterms are subtracted.Explicitly, the z → 0 (z → 1) counterterms are subtracted in the phase-space region where z < A (1 − z < A), and therefore 0 < A ≤ 1.
• B: Postive real number specifying the size of the box.
• s: Variable used to describe the invariant mass of the parton that initiates the jet.
In the code we have made this variable dimensionless by rescaling with the renormalisation scale µ 2 , i.e., s = s µ 2 .
• z: Variable encoding the momentum fraction z of one of the partons in the collinear splitting.
• φ: Variable corresponding to the azimuthal angle of the collinear splitting.
• nf: Variable specifying the number of (massless) quark flavors.This variable does not need to be set to an integer, but can be left in symbolic form.
• format: String specifying the output form of this function.One can choose between "Mathematica", "Fortran" and "C".
• file: String with the filename to which the integrand will be exported.For an empty string the integrand will be printed to the screen.
• method: This string can specify which method NIntegrate uses in Mathematica, and we refer the reader to the Mathematica documentation for the available options.For an empty string the default method of NIntegrate will be used.

Example: k T clustering algorithms
To illustrate the use of our code we now calculate the jet function for the family of k T clustering algorithms.At one-loop order, where there are at most two particles in the final state, they are clustered into a single jet if the angle between them is less than the jet radius parameter R, which for the case of an e + e − collider corresponds to a single region where E is the jet energy.The z → 0 and z → 1 limits of eq.(3.2) are described by These are no lower constraints, i.e. c − i = 0. Calculating this observable requires a box since the s = 0 line is inside the domain of integration.Since α i = 1, a rapidity regulator is not needed.The constraint in eq.(3.2) due to the measurement does not depend on φ, and so we take Φ = {}.
We now calculate the quark jet function.As eq.(2.8) is a relatively simple expression, for which the jet function can be easily calculated analytically, we will use Mathematica to perform the numerical integration over the subtracted integral by using QuarkJetN with the the 'LocalAdaptive' integration method.In the following we set µ = ER for simplicity.Note how this, since the variable s corresponds to s µ 2 , cancels the factor E 2 R 2 in the obsevable.
In [1] LocalAdaptive"; box = True; rr = False; A=0.6; B=20; From this answer it is straight forward to reconstruct that the full color-dependence of the regulated one-loop quark jet function is given by: The poles match exactly with the result by [31] and the finite term agrees up to order 10 −6 .Similar agreement is found for the gluon jet function: The accompanying Mathematica notebook contains several hands-on examples to further illustrate the use of the different functions.

Applications
To validate the method and corresponding code, the jet functions for several known examples have been checked.Some of these were used throughout the paper to explain our approach, namely the k T family of clustering algorithms (section 3.3), and angularities with respect to the WTA axis (section 2.4).In addition, we provide results in section 4.1 for the cone algorithm and in section 4.3 for the jet shape.The latter is more challenging due to its azimuthal-angular dependence, which arises because the jet axis is along the total jet momentum and thus sensitive to recoil of soft radiation.In section 4.2 we present, for the first time, the one-loop jet functions for angularities with respect to the thrust axis, taking into account recoil.Although for b > 0 this recoil is formally power-suppressed, it can be numerically large [47].

Cone jet
At one-loop order, the condition that both partons are within a cone jet in an e + e − collisions is that their angle with the jet axis is less than R (for pp 6 ).Since the jet axis is along the total jet momentum, one simply needs to consider the angle with the parton that initiates the jet, leading to the following condition As we focus on the finite term in the jet function, we fix µ = ER finding which agrees up to order 10 −6 with ref. [31].

Angularities with recoil
In this section we determine, for the first time, the one-loop angularity jet function that includes the recoil of the thrust axis due soft radiation.While this recoil is power-suppressed for b > 0, ref. [47] noted that it has a numerically large effect and presented a factorization framework to include it.The one-loop jet function we calculate here will start to contribute at NLL accuracy.This should be contrasted with the calculation in section 2.4, where we considered the angularity with respect to the WTA axis.To clearly distinguish these two cases in the notation, we will use τ n instead of e b , where n refers to the thrust axis.The setup underpinning our calculation is illustrated in figure 4.Here θ is the angle between the thrust axis n and the direction n of the initial collinear parton due to the recoil from soft radiation, which is treated as an external parameter in our calculation.The momenta of the two massless partons in the jet are denoted by p 1 and p 2 , where we The setup of our calculation.The recoil is quantified by θ. use (un)primed coordinates to denote light-cone components with respect to the n ( n) direction.Explicitly, and similarly for p 2 .Here we chose n µ = (1, 0, 0, 1) and nµ = (1, 0, 0, −1), z is the momentum fraction of the parton, s the invariant mass of the jet, and Q the center-of-mass energy of the e + e − collision.The expression in the recoiled frame follows from the definition of z and s through p − 1 = zQ and s = (p 1 + p 2 2 , as well as p µ 1⊥ = −p µ 2⊥ and the on-shell condition p 2 1 = p 2 2 = 0. Note that |p i⊥ | 2 = z(1 − z)s.The rotation between the two frames is described by | in the small θ approximation, where φ is the azimuthal angle around the n axis.The large momentum components are the same in both frames, p − i = p − i .The expression for the angularity τ n becomes where b > −1.Using the delta function trick (see section 2.2), we switch to a cumulative measurement, writing the observable as Unfortunately is it not possible to invert eq.(4.6) to obtain an analytic solution for s and subsequently extract the soft limit z → 0. We can, however, use the power-law ansatz in eq.(2.8) to find the soft behavior of the observable.Since the equation is symmetric in z → 1 − z, we focus on finding the soft behavior in the z → 0 limit.Using in eq. ( 4.5) and taking the z → 0 soft limit, we find There is a single solution for s in either of the soft limits and therefore this observable only has an upper boundary over the full range of b, i.e. c − 0 = 0.The leading terms in eq. ( 4.8) are used to solve for α + 0 and c + 0 , and differ for −1 < b < 0, b = 0 and b > 0. We will analyze the last case in some detail and only provide the solutions for the others.
Assuming b > 0, the leading behavior in the z → 0 limit of eq.(4.8) is and from this we infer where Similarly, for b = 0 we obtain For −1 < b < 0 the solution is a bit more difficult and reads  In order to use GOJet, we rescale s and choose an energy scale µ.To be able to smoothly turn off the recoil, we choose µ in terms of the angularity, µ = Q (τ c n ) 1/(1+b) .The only independent variable left is then given by k in eq.(4.11).To be complete we also give the resulting observable input for GOJet: The jet function for θ = 0 (without recoil) was calculated analytically in refs.[24,47] and we obtain the same results as can be seen in figure 5a.The error bars indicate the uncertainty from our numerical integration.Ref. [47] includes a zero-bin subtraction [58] to avoid double counting with the soft function in their factorization, which we do not include.The zero-bin subtraction depends on the details of the factorization theorem (indeed it vanishes in ref. [24]), so we do not offer this as a standard functionality of GOJet.The numerical integration for small values of b is particularly challenging (as can be seen for b = 1  8 ), because the sub-leading terms with respect to the leading soft behavior of the observable in eq.(4.10) are particularly large in this case.A more detailed discussion of this issue and a method to cope with it is presented in app.B. In figure 5b we reproduce the known results for b = 0 (broadening) and general recoil [52].Our new results for general b including the effect of recoil, are shown in figure 6.The error bars are not shown in this plot as they are negligibly small.

Jet shape
As another nontrivial example, we calculate the jet function for the classic jet shape observable, reproducing the one-loop result of ref. [36].The jet shape describes the average energy fraction z r inside a cone of angular size r around the jet axis.As in section 4.2, recoil from soft radiation displaces the jet axis from the initial parton by an angle θ.This breaks the azimuthal symmetry, requiring one to integrate over φ.We have checked that our poles match exactly with the poles in [36] for all values of θ and r.The difference between the finite term is always below 0.5%.This has been illustrated in figure 7a for gluon jets and figure 7b for quark jets.We note that run time is not an issue, as less precision is needed in phenomenological results and the distribution can be interpolated.Our calculation represents the second independent calculation of this observable and thereby delivers a useful cross check of the results of ref. [36].

Conclusions
In this paper we developed an automated approach for calculating one-loop jet functions, and provide an implementation in the accompanying Mathematica package called GO-Jet.We use geometric subtraction [55] to isolate the soft and collinear singularities.The collinear counterterm does not depend on the details of the observable, except that certain observables do not require it.We find that the soft counterterm depends on the behavior of the observable in the soft limits, which can be described by a power law.While the user must provide GOJet with this power law as input, we present a strategy to extract this in a highly nontrivial example.We employed cumulative distributions, such that observables correspond to integrating over certain regions of phase space, and thereby avoiding plus distributions.We have demonstrated our approach by reproducing the known one-loop jet function for a range of observables, and calculating, for the first time, the jet function for angularities including recoil.For broadening (b = 0 in our conventions) the effect of recoil must be kept [25], while for b > 0 it is formally power suppressed but can be numerically large [47].For b close to 0, we encountered numerical convergence issues, due to an integrable divergence.We addressed this problem by substantially improving the counterterm through a remapping.
Our approach focusses on IR-safe observables, and we did not address the IR-unsafe case.Jet functions containing IR divergences are sensitive to nonperturbative physics, and our purely partonic calculation must be supplemented by a (universal) nonperturbative function that subtracts these divergences.A prime example is initial-state jets, which are described by beam functions [59].Beam functions contain infrared divergences, which are removed by matching onto parton distribution functions, leaving finite matching coefficients.
The automated approach and code presented here provides a very useful tool, calculating jet functions at one-loop order.Very few two-loop jet functions are known, and an automated approach would allow many resummation calculations to be extended to NNLL or N 3 LL accuracy.At this order the singular limits become more complicated, the order of subtractions matter, and the parametrization of the observable in these limits will no longer be a simple power law, complicating the counterterms.

B Counterterm Mapping
In this appendix we discuss how to improve the convergence of the soft subtraction through a mapping.For simplicity, we consider only the soft singularity at z = 0, for which the finite term generated by the geometric subtraction is of the form: Here we suppressed the dependence (and integrals) over s and φ, extracting the 1/z singularity from the integrand Q, i.e. f = zQ.While this integrand is by construction integrable, poor numerical convergence may be caused by mismatch of the observable O and its soft limit O 0 .This problem can become particularly severe if O(z) has a fractional power series in z, as we illustrate below.
To improve the convergence of the integral, we apply the following mapping (to the counterterm only): This maps the interval 0 ≤ z ≤ 1 onto itself, as long as z + g(z) > 0, and the subtracted integral will remain the same as long as the function g(z) decreases faster near z = 0 than z itself, i.e., it satisfies lim The resulting curves are plotted in figure 8, highlighting the improvement due to the remapping.A Vegas run using 5 • 10 9 points for the finite part of the quark jet function of this observable yields −48.63(2) without the mapping, while we obtain −48.745(9) after the mapping.The true value is −48.7731,indicating that the remapped counterterm yields a result significantly closer to the true value.In both cases it becomes clear that the offset is not completely covered by the uncertainty.While the remapping may thus improve convergence, it may not completely solve the issue.

C Azimuthal Integral
In this appendix we evaluate the integral One can convert this integral into a Gauss-type hypergeometric integral using the transformation cos φ = 1 − 2x.However this leads to square roots in the denominator which do not naively lead to a polylogarithmic expression.Instead, one can rewrite the integral as a contour integral in the complex plane using the transformation z = e iφ , such that sin φ = z The integrand can be chosen to have branch cuts on the real axis for z < 0 and for z > 1.
For 0 < a, b < π, which is the range of physical interest, no branch cuts are ever crossed.
It is convenient to perform the integral on a contour along the real axis from 0 < z < A with 0 < A < 1, i.e., While the divergence at z = 0 requires careful treatment, this drops out in the difference of the two terms in eq.(C.5).We performed the integral using the Maple package Hyperint [60], finding that the integral can performed order by order in in terms of harmonic polylogarithms.This is to be expected, given that its singularities are located at z = 0, −1, 1.Up to order 2 we can express the result in terms of the classical polylogarithms:

Figure 1 :
Figure1: For a general observable the phase space can be constrained to several regions (blue).The collinear singularity C (red line), soft singularities S 0 and S 1 (purple lines), and soft-collinear singularities (black dots) are indicated.

Figure 2 :
Figure2: A graphical representation of our subtraction scheme in eq.(2.6).We have only included the soft counterterms for z → 1 for legibility.Shown are the restriction on the measurement from the observable (blue line), the soft limit of the observable (red line), the box (green line), the cut on z arising from A (pink line).Blue plus (minus) areas correspond to positive (negative) contributions of the full integrand Q i M obs , while red plus (minus) areas correspond to positive (negative) contributions of Q i,1 M obs,1 .

Figure 3 :
Figure 3: IR unsafe observables that our code (a) can and (b) can't handle.
we obtain the following measurement function for a cut on the angularity e b ≤ e c b ,

Figure 5 :
Figure 5: The offset between our results for (a) different values of b with θ = 0 and (b) different values of the recoil parameter k with b = 0 and the known results from the literature is shown.

Figure 6 :Figure 7 :
Figure 6: The results for the finite part of J 1 q for different values of b as a function of k.