Jet Veto Clustering Logarithms Beyond Leading Order

Many experimental analyses separate events into exclusive jet bins, using a jet algorithm to cluster the final state and then veto on jets. Jet clustering induces logarithmic dependence on the jet radius R in the cross section for exclusive jet bins, a dependence that is poorly controlled due to the non-global nature of the clustering. At jet radii of experimental interest, the leading order (LO) clustering effects are numerically significant, but the higher order effects are currently unknown. We rectify this situation by calculating the most important part of the next-to-leading order (NLO) clustering logarithms of R for any 0-jet process, which enter as $O(\alpha_s^3)$ corrections to the cross section. The calculation blends subtraction methods for NLO calculations with factorization properties of QCD and soft-collinear effective theory (SCET). We compare the size of the known LO and new NLO clustering logarithms and find that the impact of the NLO terms on the 0-jet cross section in Higgs production is small. This brings clustering effects under better control and may be used to improve uncertainty estimates on cross sections with a jet veto.


Introduction
An important ingredient in precision measurements of Higgs boson properties is the Standard Model prediction for cross sections used in the experimental analyses. These analyses often involve cuts on the final state hadronic activity, vetoing candidate jets with transverse momentum p T above a veto scale p cut T . Exclusive jet cross sections, such as the H + 0-jet cross section in which there are no jets with p T > p cut T , play a key role in channels where the Higgs cannot be reconstructed and jet binning is an effective discriminant between various Standard Model backgrounds. The uncertainty in theoretical predictions of these exclusive jet cross sections can be substantial, and accurately predicting these cross sections and providing robust estimates of their uncertainty is the focus of considerable effort [1][2][3][4][5][6][7][8][9][10][11][12]. A good example is the H → W W * → + − νν channel, where a fixed order analysis of the exclusive 0-jet and 1-jet cross sections have theoretical uncertainties of O(17%) and O(30%) respectively and dominate the systematic uncertainties in the measurement [13,14].
There are two challenges in understanding cross sections with a jet veto. The first is the fact that the jet veto scale, which is typically between 25 and 40 GeV, is well below the hard scales in the process, such as the Higgs mass m H . This leads to large logarithmic corrections of the ratio p cut T /m H . Fixed order perturbation theory supplemented with resummation of these jet veto logarithms is a powerful tool to make predictions of exclusive jet cross sections and the associated uncertainties. Recent work in this vein has substantially lowered the H → W W * 0-jet and 1-jet uncertainties [6,[8][9][10][11].
The second challenge is the effect of jet clustering. Clustering from the jet algorithm leads to a complicated dependence on the jet radius R, a dependence that also takes part in and complicates the resummation of jet veto logarithms. These effects start at O(α 2 s ) 1 , the first order in which there can be two final state partons and nontrivial clustering can take place. At O(α 2 s ) the clustering effects have been fully determined and have a substantial effect on the cross section [6,10]. No higher order terms are known, meaning the contribution of jet clustering to the cross section is effectively only known at lowest order. It is challenging to make reliable uncertainty estimates of the higher order clustering effects from only the leading order term. Calculating clustering contributions at the next order is the focus of this work.
We take the 0-jet cross section in Higgs production via gluon fusion as a test process, although results here can be reapplied to other color-singlet processes through a simple exchange of color factors. The clustering effects have two types of terms: those proportional to logarithms of m H /p cut T ("single logarithmic terms") and terms independent of p cut T ("finite terms"). At O(α 2 s ) all of the clustering terms have been determined, and both single logarithmic and finite terms have contributions proportional to ln R [3,5].
In fact, the general form of the single logarithmic terms is known [5]: at each order, there is a new contribution multiplying the resummed cross section of the form U (n) clus (R, p cut T ) = exp where σ 0 (p cut T ) = σ sing 0 (p cut T ) + σ ns 0 (p cut T ) , σ sing 0 (p cut T ) ∝ U The exponentiation in eq. (1.1) comes from the resummation of the veto logarithm, and C n (R) is a function that can be decomposed in terms of logarithms of R as C n (R) = C (n−1) n ln n−1 R 2 + . . . + C (1) n ln R 2 + C (0) n (R) , where C (0) n (R) is finite as R → 0. Furthermore, the coefficients of the logarithms of R can be connected to finite terms with logarithms of R [10]. For a qq initiated color-singlet process, one factor of C A is replaced with a C F in eq. (1.1).
Standard jet radii used in experimental analyses are R = 0.4, 0.5. Formally, if one considers R ∼ p cut T /m H then the logarithms of R are as important as the veto logarithms and resummation should be performed. This implies that the leading coefficient C (n−1) n at every order is part of a next-to-leading logarithmic (NLL) series. It is not known if any of these coefficients are related, but if they are not then they have to be explicitly calculated at each order, rendering resummation 2 of the entire series impossible. Note that the exponentiation of each term is known, as in eq. (1.1).
At O(α 2 s ), the leading clustering term has the coefficient [3] C For p cut T ∈ [25,30] GeV and R ∈ [0.4, 0.5], this term makes a contribution to U (2) clus which increases the cross section by an amount between 9% and 15%, which is on par with, or larger than, the theoretical uncertainties on the most recent predictions! Clearly the higher order terms are important to not only provide precise predictions, but to understand the uncertainty associated to these clustering effects. 2 Here we want to distinguish between exponentiation and resummation. Exponentiation involves capturing all higher-order terms (usually at a certain logarithmic order) that are generated by a given term, as in eq. (1.1). Resummation is more complex and involves capturing all terms at a logarithmic order. In the case here we can exponentiate a known C (n−1) n to capture the NLL terms originating from it, but we cannot resum the NLL clustering logarithms unless all C (n−1) n are known.

Overview
In this work we calculate the coefficient C (2) 3 , which is the coefficient of the ln 2 R ln m H /p cut T term in the O(α 3 s ) correction to the cross section. At first glance, the calculation of C (2) 3 is very challenging as it naively requires a three loop calculation. However, since the clustering effects start at O(α 2 s ), the full coefficient C 3 (R) is intrinsically an NLO quantity. There are three major simplifications that blend nicely in the calculation of C (2) 3 . They are: • Collinear factorization formulas that simplify the matrix elements relevant for the calculation of C 3 . One can exploit the fact that the clustering logarithms arise from the region of phase space where final state particles are close together. This means the leading clustering term at each order (given by the coefficient C (n−1) n ) is determined by matrix elements with a collimated final state, allowing one to take advantage of powerful collinear factorization formulas [15,16]. This allows us to attack C (2) 3 rather than the entire C 3 (R).
• Factorization properties of soft-collinear effective theory (SCET) [17][18][19][20][21], which allow for a direct calculation of C 3 via the soft function. Factorization properties from SCET allow the fixed-order cross section to be factorized into beam and soft functions, where beam functions describe radiation along the beam directions (see Ref. [22]) and the soft function describes global radiation across the whole event. This factorization allows one to cleanly separate the dynamics that gives rise to the ln m H /p cut T multiplying C can be framed as an NLO calculation in the soft function (one where the tree level contribution starts at O(α 2 s ) in the soft function).
• Subtractions for NLO calculations, which separately regulate the infrared divergent regions of phase space in the real and virtual contributions and allow one to perform the calculation numerically. Subtractions are an efficient method to handle the divergences present in real and virtual matrix elements. These divergences cancel in the total result, but since the phase space of the real and virtual contributions are distinct, a numerical implementation of the calculation that implicitly sums these canceling divergences is challenging. Subtractions separately render the real and virtual contributions finite by reproducing the divergences in the real matrix element in a way that can be analytically integrated to provide a cancellation of the divergences in the virtual matrix element.
The first two simplifications work in tandem to produce an NLO calculation that determines C (2) 3 . However, the phase space restrictions required for the calculation make the singularity structure quite complex, and the matrix elements are still lengthy. The calculation is made tractable by using subtractions. Using the simplifications in the matrix element from QCD and SCET techniques, the number of singular regions of phase space is substantially reduced from a complete H + 2-jet NLO calculation. When introducing the subtractions, we will highlight the connection between the matrix elements used for subtractions and the matrix elements of SCET.
In sec. 2 we show how the calculation of C 3 (R) can be reduced to a soft function calculation, and in sec. 3 we determine the matrix elements needed to calculate C (2) 3 using splitting functions. In sec. 4 we set up the calculation of C 3 as an NLO calculation using FKS subtractions. In sec. 5 we determine the relevant terms in the virtual matrix elements needed for the calculation, as well as their combination with the integrated subtractions and the contribution to the result. In sec. 6 we perform the calculation of C (2) 3 , analyzing the impact on the H + 0-jet cross section and its uncertainties and discussing the all-orders series of clustering logarithms. Finally, in sec. 7 we give our outlook and conclude.

Clustering Logarithms from the Soft Function
The veto on jets found by a clustering algorithm can be viewed as a local veto, because the jet algorithm acts on local clusters of radiation. This should be contrasted with a global jet veto, an event-wide measure that restricts jet activity that does not depend on a jet algorithm. For example, the H + 0-jet cross section can be defined by a veto on beam thrust, and precision calculations can be carried out to high accuracy without the complications that arise in a clustering algorithm [1]. Relevant to our case, the E T of the event is an effective global veto that is independent of the jet algorithm, where At O(α s ) there is only one particle in the final state and requiring E T < p cut T is equivalent to the jet veto. The difference is useful to understand clustering effects. In the first term the jet veto is performed (where p jet T is the leading jet p T ), and in the second a global veto on E T is performed. The global veto term is independent of the jet algorithm and any clustering effects, and may be used to understand the structure of the veto logarithms. In this case it is helpful to isolate clustering effects, and hence the R-dependence, into the ∆σ clustering correction 3 .
Since the O(α s ) clustering correction vanishes, ∆σ (1) (p cut T ) = 0, the known O(α 2 s ) clustering effects are a LO quantity. This implies that ∆σ (3) (p cut T ), which contains the complete C 3 (R), could be determined from existing H + 2-jet NLO codes. However, ∆σ also contains higher powers of veto logarithms arising from exponentiation of lower order terms, and extracting C 3 (R) requires an overwhelming computational investment (that is not very feasible) as it requires working in a small R and p cut T regime to accurately extract the logarithmic dependence in C 3 (R). Therefore we opt to focus solely on the leading ln 2 R terms and calculate them through more direct methods.
In the small R limit, the H + 0-jet cross section can be factorized in SCET into hard (H), beam (B a,b ), and soft (S) functions [4,5]: The bare soft and beam functions individually contain rapidity divergences that are not regulated by dimensional regularization. These rapidity divergences can be regulated in different ways, and in this work we use the rapidity renormalization group [23,24]. This method regulates the rapidity divergences with an explicit factor in matrix elements, introducing a scale ν that functions much like µ in dimensional regularization.
In SCET with the rapidity regulator the clustering logarithms from C n (R) are divided between the beam and soft functions. At O(α n s ), the contributions from C n (R) are [5] total ∝ α s 4π Therefore C n (R) can be calculated from the soft function, simplifying the computation.

The Soft Function
The soft function in SCET is a forward scattering matrix element of soft Wilson lines with a veto measurement on the final state, The Y n are soft Wilson lines, and n = (1, 0, 0, 1),n = (1, 0, 0, −1) are light-like vectors oriented along the beam direction. The measurement function for the veto cross section is where the jets are formed by clustering the soft particles in the final state. For the clustering correction (relative to the E T global veto), the measurement function is which defines a soft function correction ∆S. At O(α s ), ∆M vanishes, and at O(α 2 s ) and O(α 3 s ) in the small R limit 4 it is ∆M (2) The p T i are the transverse momenta of the different particles, and at O(α 3 s ) the outcomes of the jet algorithm are classified by constraints of the form θ(2-jet; {1+2, 3}), which for example requires the algorithm to yield two jets, one with particles 1 and 2 and the other with particle 3. These are jet algorithm dependent, and here we use the k T -type clustering algorithms (k T , Cambridge/Aachen, and anti-k T [25][26][27][28][29]), which allow us to study the algorithm dependence of the coefficient we calculate. The O(α 2 s ) phase space constraint depends only on the separation ∆R 12 = ∆y 2 12 + ∆φ 2 12 between the final state particles, and is common between k T -type clustering algorithms. The constraint θ(∆R 12 > R) requires the two particles to be in different jets. A study of the singular limits in these measurement functions shows that ∆M (2) vanishes in the soft or collinear limits, and ∆M (3) is nonvanishing only if a single parton becomes soft, a pair of partons become collinear, or in the combined soft-collinear limit. Hence the calculation is LO at O(α 2 s ) and NLO at O(α 3 s ). From the expressions in eq. (2.8), one sees that the measurement function vanishes if either of the following is true: These serve as useful bounds on the phase space. We note that this discussion is independent of the soft function, and also applies when considering the full QCD calculation. The soft function is invariant under boosts along and rotations around the beam direction, suggesting a certain set of phase space variables. Each on-shell (final state) phase space integral can be expressed in terms of the transverse momentum p T , the rapidity y, and the azimuthal angle φ.
Since the matrix elements and measurement are independent of the total rapidity (y t ) and azimuthal angle (φ t ), a useful combination of the above variables is the following, Note that it is also possible to define a dimensionful variable for the total p T and rescale each transverse momentum by it; since this is the only dimensionful phase space variable its dependence in the matrix element is fixed. We will find use for these variables, and define For the phase space integration we find it is more useful to rescale by the fixed p cut T and keep the x i . The phase space in terms of our chosen variables is For those not familiar with soft functions it is worthwhile to note that an arbitrary amount of momentum can go into the final state, as evidenced by the fact that there is no momentum conservation relation present in the phase space.

Projecting Out the Veto Logarithm
Consider the soft function contribution at a given order, with the d-dimensional, -loop matrix element denoted T ( ) a 1 ...an (Φ n ), where all spin and color have been summed over and the a i serve as final state parton labels. We only include the single c-web terms in the matrix element, as these are precisely the terms contributing to C n (R) [5,30]. Note that the number of parton labels a i denote how many particles are in the final state (the remaining are loop momenta which are integrated over). The bare soft function clustering correction contribution from these matrix elements is where Sym a 1 ...a k is the symmetry factor for the particular channel. The only non-unit symmetry factors relevant to this calculation are Sym gg = 1/2! and Sym ggg = 1/3!. For simplicity, we first consider more detailed properties of the fully real matrix element, T n . We can make the dependence of the matrix element on the dimensional regularization scale µ, the coupling, and the parent phase space explicit: The measurement function only depends on Φ n . One notices that neither the matrix element nor the measurement function depends on y t and φ t variables in the parent phase space, and it seems we can integrate over them with impunity. However, since the y t integral spans an infinite range, we get an unregulated divergence. This is a rapidity divergence, and there are canceling rapidity divergences in the soft and collinear sectors that must each be regulated. The rapidity renormalization group provides a regularization scheme to regulate rapidity divergences and renormalize them much in the same way that dimensional regularization regulates virtuality divergences [23,24]. Resummation of rapidity logarithms using the rapidity renormalization group follows familiar steps to conventional resummation. For a single c-web, the regulator factor in the soft function is where p 3g is the component of the group momentum for the c-web along the Wilson line direction (the beam direction). Since we only deal with a single c-web, the group momentum is the sum of all momenta in Φ n . For our purposes, we only need the leading singularity: This leading divergence produces a finite term proportional to ln ν/p cut T , which is the logarithm that appears in eq. (2.4). It is useful to note that the O(η 0 ) term contains kinematic factors that regulate the combined collinear limit of the c-web (which implies it will only contribute to subleading clustering logarithms).
Integrating over Φ t , the result is Note that there is no additional divergence in the rapidity regulator parameter η besides the overall 1/η in eq. (2.20). If we are only concerned with the term containing the logarithm of p cut T , we can expand in η and obtain the coefficient of ln ν/p cut T , which is directly related to C n (R) [see eq. (2.4)]. There are remaining divergences in , though. The measurement function cuts off the ultraviolet region of phase space, meaning they are infrared singularities which cancel between real and virtual contributions.

The Real Matrix Elements in Terms of Splitting Functions
The leading clustering logarithm dependence [ln 2 R at O(α 3 s )] arises from the collimated limit of the final state. This implies that if we want to calculate C (2) 3 alone, we can use matrix elements in this limit. At tree level, these matrix elements can be built from lower order amplitudes by exploiting collinear factorization of tree-level amplitudes [15,16,[31][32][33]. This factorization manifests at the lowest level in terms of the factorization of color-ordered amplitudes for individual helicity configurations. At a less exclusive level, this factorization can be phrased in terms of splitting functions. In fig. 1 we give a schematic picture of the matrix elements.  We note that to calculate the single ln R terms at O(α 3 s ), whose coefficient is C 3 , the matrix elements outside of the fully collimated limit are needed. In general, the number of logarithms of R is given by the number of collinear propagators in the matrix element [5], meaning that the ln R terms arise from only one pair of partons becoming collinear (instead of all of them). Since the O(α 3 s ) matrix elements used here are in the triple collinear limit, they will not capture the complete C (1) 3 . We adopt the basic notation in [15,16]. If A c 1 ,...;s 1 ,...
is the tree-level amplitude to produce momenta p i with parton flavors a i , colors c i , and spins s i , then the matrix element squared, summed over color and spin, is If a set of momenta {p 1 , . . . , p n } become collinear, then the matrix element factorizes, The momentum p is the collinear limit of p 1 + · · · + p n with flavor a. The spin-polarization tensor T ss is obtained by summing over all other spins, The momenta p i can be decomposed in the collinear limit as wheren p = (1, −p) is a null vector and the component of each momentum transverse to the collinear direction is Note that this decomposition satisfies i z i = 1 and i k µ ⊥i = 0. In our case, we are interested in the entire final state being collimated. This implies that the parent matrix element is that for single gluon emission (which then undergoes a collinear splitting), whose spin-polarization tensor is For a qq initiated process, the color factor C A is traded for C F . Note that (−g µν )E µν = 1.
We will contract this spin-polarization tensor with the O(α 2 s ) g → gg, g → qq and the O(α 3 s ) g → ggg, g → gqq splitting functions.

Splitting Functions
The 1 → 2 and 1 → 3 polarization-dependent splitting functions are given in Refs. [15,16]. One can also perform the average over polarizations [33]. For gluons, The polarization tensor D µν is valid for axial gauge withn p · A = 0. In principle, the polarization dependent terms can make a nonzero contribution to C 3 . However, this is not the case; we will show that the polarization dependent terms depend on orientation angles of the collinear system relative to the rest of the event that C (2) 3 is insensitive to, and average to zero upon integration over the total phase space 5 . This is reasonable since the measurement functions that determine C 3 depend only on the relative positions of the collinear partons. To see the polarization dependence explicitly, we will find the difference to the polarization averaged splitting functions, defininĝ In app. A we give the relevantP ss a 1 ...an and P a 1 ...an , which originally appeared in Ref. [15,16], and study ∆P a 1 ...an here.

1 → 2 Splitting Functions and the Born Matrix Elements
The 1 → 2 gluon-initiated splitting functions are given in eq. (A.1). The polarizationdependent terms are If p 1 and p 2 are the collinear momenta, then z = x 1 /(x 1 + x 2 ) and ∆y, ∆φ are the rapidity and azimuthal angle separations. The separations ∆y and ∆φ make up the angle ∆R, and we can define an angle θ which expresses how the collinear system is oriented relative to the transverse plane, The measurement function and the polarization-averaged matrix elements are independent of θ, which means that its only dependence is in the polarization dependent term ∆P. In the collinear limit, the integration over the full phase space will average over this angle and the polarization dependent terms will vanish: This justifies our use of the polarization-averaged O(α 2 s ) splitting functions. Hence the total Born matrix elements are (in 4 dimensions) If we include the symmetry factor of Sym gg = 1/2! in the gg matrix element, these agree with the known O(α 2 s ) soft function matrix elements when taken into the collinear limit [34]. These can be used to calculate the O(α 2 s ) clustering logarithm coefficient C 2 , as shown in app. C.

1 → 3 Splitting Functions and the Real Matrix Elements
The story above repeats itself with the O(α 3 s ) splitting functions. The g → ggg and g → gqq splitting functions [given along with the polarization averages in eq. (A.2)] have polarization dependent terms that depend on the orientation of the collinear system and average to zero upon integration over the phase space. These polarization dependent terms, ∆P, have no contribution from terms inP µν proportional to g µν , since g µν (E µν − D µν ) = 0. All other terms have 3 basic structures, and we list their projection with E µν − D µν : z 1 s 23 cos 2θ 23 + z 2 s 13 cos 2θ 13 − z 3 s 12 cos 2θ 12 − 1 2 z 1 z 2 s 12 cos 2θ 12 + s 13 cos 2θ 13 + s 23 cos 2θ 23 , Other choices of the indices on k ⊥ have the same form. Each term contributing to ∆P is proportional to one of these structures, and each structure is in turn proportional to cos 2θ ij for some i, j. This implies that, like at O(α s ) 2 , the polarization dependent terms cancel in the contribution to C 3 . The matrix elements for real emission are thus given in terms of the polarization averaged splitting functions, (3.14) Note that there is no kinematic hierarchy used within the splitting functions (such as a strongly ordered limit), and that the complete 1 → 3 splitting functions are used in the matrix elements.

Calculating Clustering Logarithms with Subtractions
In the NLO calculation that determines C 3 , there are canceling divergences in the real and virtual matrix elements, and the phase space cuts implemented by the measurement function are complex. This situation is well-suited to use a subtraction to separately regulate the real and virtual contributions. The universal soft and collinear factorization properties we have exploited to determine the matrix elements are precisely those which allow for subtractions that can regulate the real and virtual divergences simultaneously.
The basic form of a NLO calculation is where V 2 and B 3 are the virtual and real matrix elements. The soft, collinear, and collinearsoft divergences in the real emission cancel with the virtual contribution. If we can define a set of subtraction terms S i (Φ 3 ) that match the singularities of the real matrix element, then they can regulate both divergences at once: where Note that the subtraction terms always come with the Born measurement function; this is because any infrared safe measurement cannot resolve a soft or collinear splitting in the singular limit (and hence depends on one fewer degree of freedom). These Born events may differ for different singular limits. Typically, the projection from Φ 3 onto a Φ i 2 occurs via a map, using a phase space factorization formula to write Φ 3 in terms of Φ i 2 and radiation variables. In our case, however, the parent gluon in the splitting can have arbitrary momentum, meaning that the phase space naturally factorizes into radiation variables for each emission (basically making the map trivial).
The subtractions we use are essentially a variation of FKS subtractions, which separately handle soft, collinear, and collinear-soft singularities [35,36]. In FKS subtractions, these singularities are isolated into regions where at most one collinear and one soft singularity are present [35]. This is done by introducing the factor S ij is nonvanishing only if partons i and j become collinear or i is soft (or the simultaneous soft-collinear limit). The S ij satisfy i,j S ij = 1 and the following limits:  This machinery is useful once we pair it with the subtraction terms. In our case, we will be able to write the singular limits of the matrix element in terms of factors that multiply the Born matrix element. These subtraction factors are: S i kl : soft radiation of particle i between particles k and l , C ij : collinear splitting into particles i and j (i < j) , CS ij : soft-collinear splitting into particles i and j, with i soft .
And the real matrix element becomes, under specific singular limits, (4.7) The above assumes that a 3 = g, so that soft and collinear-soft limits are actually singular. We can take the singular limits of the real matrix elements analytically and define the subtractions in terms of plus distributions. We will explicitly use the g → ggg matrix elements as an example in the rest of this section, and the general case easily follows. Integrating the subtractions is performed in app. B.
Each subtraction will regulate a particular divergence in the real matrix element, and the integrated subtraction compensates for the terms that are introduced. It is often very useful to introduce auxiliary cuts on the singular variables in the subtractions and integrated subtractions, which serve as a self-consistency check on the calculation (the dependence on these cuts should cancel in the total cross section) and can help probe the singularity structure (in our case, the dependence on R). In the usual FKS subtractions, these cuts are on the soft gluon energy in soft and collinear-soft subtractions and the splitting angle in the collinear and collinear-soft subtractions. For example, instead of integrating over all soft gluon energies E g , one integrates over the range 0 < E g < E c . The dependence on the artificial parameter E c must cancel between the subtraction and integrated subtraction. In our case, we will find the following cuts useful: x g < x c : for soft subtractions, with x g the soft gluon p T fraction , ∆R < R c : for collinear subtractions, with ∆R the opening angle of the splitting . (4.8) The collinear-soft subtraction will have both of these cuts, and the integrated subtractions will depend on x c and R c in the appropriate places.
The matrix elements for the subtractions are closely related to matrix elements in SCET. This is not surprising, as both are built from the singular limits of QCD: we show that the soft, collinear, and collinear-soft FKS subtractions are given by soft, naive jet, and zero-bin jet function matrix elements. We will explore these connections further in this section.

Soft Subtractions
In the limit that gluon 3 becomes soft, the g → ggg matrix element becomes This expression can be understood as a sum of eikonal factors multiplied by the Born matrix element for gluons 1 and 2. The first term is the eikonal factor for soft gluon exchange between gluons 1 and 2 (proportional to the Born matrix element), which defines the subtraction term, . (4.10) For g → ggg, the color operators T 1 and T 2 obey the relations T 2 where T t is the color operator for the initial gluon that splits. It is also useful to define the color operator T r ≡ −T t for the rest of the event, which also obeys T 2 r = C A . The subtraction term in eq. (4.10) can be recognized as the O(α s ) soft function matrix element for soft gluon exchange between soft Wilson lines 1 and 2. This is natural, as both the subtraction and the soft function matrix elements exploit the eikonal factorization properties of QCD amplitudes.
The second and third terms in eq. (4.9) come from soft gluon exchange between either gluon 1 or gluon 2 and the rest of the event. Because the g → ggg system is color-connected to the rest of the event, other partons can radiate soft gluons into the collinear system. If there is another parton i that exchanges soft gluon 3 with gluon 1, then the usual eikonal factor is (4πα . (4.11) However, since gluon 3 must be collinear with gluon 1, θ 13 θ i3 and and thus the eikonal factor reduces to (using the collinear limit) . (4.13) The kinematic factor is i-independent, and so we can sum over all colors. This is an example of coherent soft gluon emission by the rest of the event, and is well-studied in related contexts [37]. This yields , (4.14) where −2T 1 · T r = C A for g → ggg, which matches the term in eq. (4.9). Note that this subtraction term is actually in a collinear-soft limit. This is arising in the soft subtraction because we are demanding the final state partons are collimated, necessitating the additional expansion. Thus, these subtractions will reappear in the collinear-soft case (where they remove double counting of divergences with the collinear subtraction), and in that case the color connections are different (so we will retain distinct labels for S 3 1r and CS). The matrix element in eq. (4.14) can be recognized as the zero-bin matrix element in the jet function, which coincides with the soft function matrix element taken into the collinear limit (where the soft gluon is collinear to one of the Wilson lines). This combined collinear-soft limit has also been studied previously in SCET [38], and can be understood in terms of an additional factorization in precisely the collinear limit that we are studying. That is, soft gluons in the splitting function also have collinear scaling, and Ref. [38] terms them csoft gluons. Since the emission of these csoft gluons is controlled by the collinear system, as well as coherent soft radiation from the rest of the event, their emission factorizes at the level of operators in SCET using fundamentally the same collinear factorization properties used to write the overall matrix elements. The total soft gluon emission matrix element, in eq. (4.9), is given precisely by the csoft function as formulated in Ref. [38]. In the csoft function the "rest of the event" is represented by a soft Wilson line (V ) that radiates soft gluons from the anti-collinear direction.

Collinear Subtractions
In the collinear limit of gluons 1 and 3, the matrix element takes the form where the azimuthal terms are dependent on the polarization of the 1-3 collinear system relative to gluon 2 [39], and vanish in the total contribution to C 3 . The polarizationindependent terms are given by tree-level collinear factorization applied to the 1-3 leg of the underlying Born matrix element T (0) gg (k 1 + k 3 , k 2 ), and the subtraction factor is More generally, we need to choose the appropriate splitting function for the process, This subtraction term may be recognized as the O(α s ) jet function matrix elements (before a zero-bin subtraction [40]). As with the soft subtraction, this is expected as the jet function is based on the same collinear factorization properties (see Ref. [41]) as the collinear subtraction terms.

Collinear-Soft Subtractions
The collinear-soft limit can be taken either from the soft or collinear limits above, and there are particular soft subtractions (those where the soft gluon is exchanged with the "rest of the event") whose kinematic dependence matches the collinear-soft subtraction (see sec. 4.1). The collinear-soft limit accounts for the double counting of divergences between collinear and soft, and is Above P (0) gg is the soft limit of the g → gg splitting function. The subtraction factor is, in terms of splitting functions, and more generally the g → gg splitting function is replaced by the right one for the process. It can also be written in the form of eq. (4.14) (with a different color factor), . (4.20) As discussed in sec. 4.1, this subtraction term is given by the zero-bin matrix element in the jet function, or equivalently by the O(α s ) soft function taken into the limit where the soft gluon is collinear to one of the Wilson lines.

The Regulated Real Emission
Having detailed the subtractions, we can now use them to regulate the real matrix elements. The sum over all regions is Sym a 1 a 2 a 3 T (0) a 1 a 2 a 3 = (S 12 + S 21 + S 13 + S 31 + S 23 + S 32 ) Sym a 1 a 2 a 3 T (0) a 1 a 2 a 3 . For the g → ggg case, the exchange symmetry means that we can drop the symmetry factor by choosing one region, This equivalence is obviously not point-by-point, but is true integrated over phase space. This limits the number of subtractions we need, since S 31 T ggg is singular only in a fraction of the limits that T (0) ggg is. Using the limits in eq. (4.6), the ggg channel is regulated in the combination 23) and the contribution to the soft function we call R ggg . Note that the soft and collinear limits of the S ij , S soft 31 and S coll 31 , average to 1/2 in integrating over the Born phase space. For the gqq case all the subtractions must be included (although only the gluon can become soft or collinear-soft). For this channel we switch to the labels g, q, andq. In each sector we take the limits for each subtraction, and the regulated combination is We find that R ggg and R gqq are indeed finite when integrated over phase space. These regulated pieces will be calculated numerically, and the calculation discussed in sec. 6.

Virtual Matrix Elements and Integrated Subtractions
So far, we have made use of several factorization properties to define the real matrix elements in terms of splitting functions. Extracting the corresponding virtual contribution is more involved. We are not calculating a "complete" cross section, in the sense that we are neglecting subleading terms in ln R, and hence only including matrix elements singular in the appropriate limit. Our approach is to focus only on the terms which give a contribution to C (2) 3 , which are directly related to one-loop splitting amplitudes.
In this section we will show that the only terms in the virtual matrix elements that can make a contribution to C (2) 3 are particular terms that come from the one-loop splitting amplitudes. Pairing their dependence with the integrated subtractions we can determine the contribution to C (2) 3 analytically.

The Virtual Matrix Elements
The collinear factorization property in eq. (3.2) extends beyond tree level, although more naturally in terms of amplitudes. We use the notation in Ref. [15,16]. If |M (0) n+1 (p 1 , p 2 , . . .) is the n + 1-particle tree-level amplitude which is a vector in color and spin space, then in the collinear limit of particles p 1 and p 2 the factorization is where P is the collinear limit of p 1 + p 2 . The splitting matrix Sp (0) is a matrix in color and spin space, and is related to the tree-level splitting amplitude Split tree via a simple color matrix. At one loop, the virtual corrections satisfy a similar property (see, e.g., Ref. [42]): That is, the loop corrections to the unfactorized amplitude divide into loop corrections to the splitting matrix and loop corrections to the underlying n-particle amplitude. This factorization is divided into two terms, and the crucial property in our case will be that the s 12 dependence is isolated into the tree-level and 1-loop splitting matrices. Furthermore, the 1-loop amplitudes for the n-parton configuration (where the collinear partons have been clustered) satisfy the decomposition n is a color operator containing the IR poles and |M (1) fin n is finite as → 0. The same decomposition holds for the 1-loop splitting amplitude,  The last two lines are finite as → 0 and the only singular dependence on s 12 comes from the 1/s 12 present in the splitting matrices. It will become evident below that we can neglect these terms. The first line is a reflection of the fact that the IR poles are proportional to the Born amplitude. Indeed, the matrix element squared is (1) The integrated subtraction is also proportional to the Born matrix element, and the IR singularities cancel at the level of the prefactor to the Born matrix element. That is, since the virtuals and integrated subtractions share the Born matrix element as a prefactor to the poles, including subleading terms in (that vanish when → 0), we do not have to worry about 1/ 2 , 1/ poles multiplying subleading terms in in the Born matrix element to generate finite terms. This is well known, but it is a crucial property that we must exploit in determining the contribution of the virtual matrix elements to C 3 . Let us consider what kinematic dependence in the virtual and integrated subtractions can give rise to ln 2 R terms. Since both contributions are proportional to the Born matrix element, the relevant question is what kinematic dependence is required in the prefactor. It is clear that the logarithms of R arise from the angular phase space integrals; in the Born contribution the relevant integrals are The 1/∆R 2 factor arises from the propagator of the collinear parton that splits, which gives a factor of 1/s. To get another logarithm of R, we must have a factor of ln ∆R, as In the virtual matrix elements, the ∆R 2 dependence arises only from the invariant mass s 12 of the collinear pair of partons. Therefore, only finite terms in the virtual matrix elements that contain a logarithm of s 12 can contribute to C 3 . This means we only need to track this dependence in the virtuals, and do not need the complete one-loop matrix elements. Since the splitting matrix Sp (1) , and in particle I C , is the only part of the virtuals that depend on s 12 , we can focus only on this term.
The divergent terms in the (UV-renormalized) 1-loop splitting matrix are [43] I (1) . The constants C i are the Casimirs for each parton in the splitting, with C g = C A and C q = C F , and the constants γ i are For the gg channel, 1 = 2 = 12 = g and so the relevant factor in the virtual matrix elements is, including the symmetry factor Sym gg = 1/2! for the gg channel, g − ln 2 ∆R 12 + ln ∆R 12 γ g − 2 ln x 1 x 2 x 1 + x 2 + ∆R 12 -independent, finite . gg . The integrated subtractions must cancel the ∆R 12 dependence in the divergent terms as well as the finite ln 2 ∆R 12 term (which would generate a ln 3 R). The qq channel has 1, 2 = q,q and 12 = g, meaning the relevant factor in the virtual matrix elements is (5.13)

Adding in the Integrated Subtractions
Similar to the virtual corrections, the terms in the integrated subtractions that can contribute to C (2)  3 are those with an additional ln ∆R 12 multiplying the Born matrix element. The integration subtractions are found in app. B, and it can be seen that only the soft subtractions with soft gluon exchange between particles 1 and 2 depend on ∆R 12 . Since soft gluon emissions from the gg and qq Born configurations contribute to the ggg and gqq channels respectively, this means the integrated soft subtractions in the ggg channel should be paired with the gg virtual matrix elements, and the gqq integrated subtractions paired with the qq virtual matrix elements 6 .
Using eqs. (4.23) and (4.24) to determine the relevant ∆R 12 dependent terms, we find 1 ln ∆R 12 + ln 2 ∆R 12 + 2 ln x c ln ∆R 12 , (5.14) which multiplies T (0) gg , and qq . Adding these integrated subtractions to the virtual matrix elements in eqs. (5.12) and (5.13), we obtain the regulated virtual matrix elements V S gg and V S qq that contribute to C 3 . For the gg channel, we have and for the qq channel, The x c dependence, which is a parameter of the subtraction formalism, must cancel between the regulated gg virtual and ggg real contributions and the qq virtual and gqq real contributions, and will provide a check on the calculation.
These matrix elements can be analytically integrated, which we do in the next section when we describe the total calculation of C 3 .

Calculation of C (2) 3
To calculate C 3 , we must take the regulated real and virtual matrix elements determined in the previous sections and integrate them over the phase space against the measurement functions. 6 If we were accounting for all divergences, we would have to take into account the fact that the gqq channel can have divergences originating from the gg channel via a g → qq splitting, although these do not contribute to C 3 .

The Regulated Virtual Contributions
The regulated virtual matrix elements in eqs. (5.16) and (5.17) are integrated over the 2 particle phase space as shown in eq. (2.20), with Sym a 1 a 2 T (0) a 1 a 2 replaced by V S a 1 a 2 . The evaluation of the virtual contributions follows closely the calculation of C (1) 2 , which is given in app. C. The only difference is the angular integral over ∆R 12 which returns a double logarithm of R and is given in eq. (5.8), and integrals over x 1 , x 2 . There are three integrals needed for each channel: The values of these integrals are, for the gg channel: T F n f (−305 + 12π 2 + 174 ln 2 + 72 ln 2 2) .
2 and k (3) is proportional to a clustering logarithm constant s (1) 2 in Ref. [10]. Using these results, the virtual contributions to the soft function are The contributions to C (2) 3 can be extracted from these results by pulling out a factor of (α s C A /π) 3 ln(ν/p cut T ) ln 2 R 2 . We note that these contributions are the same for all k T -type jet algorithms, as all algorithms in this group have the same phase space constraints with two particles in the final state. The calculation of the regulated real contributions, R ggg and R gqq , is performed nu-merically using the integration routine VEGAS in the CUBA library [44]. The calculation requires approximately 2 · 10 9 events for a given set of parameters (R, R c , x c ) sampled in the 7-dimensional phase space to achieve uncertainties on C 3 below 5%. The calculation of R ggg and R gqq uses FKS-type subtractions. For the ggg channel symmetry implies there is only one unique sector, while for the gqq channel there are two: qq and gq (which equals gq). In the gqq channel we also find it useful to divide the result into color factors (C 2 A n f and C F C A n f ), meaning there are four separate contributions to the gqq result. For each sector, the contribution to C We use the anti-k T algorithm as our primary jet algorithm, but in fig. 3 we investigate the jet algorithm dependence of the calculation. Since the virtual contributions are the same for all k T -type algorithms, the algorithm dependence of C (2) 3 is probed through the regulated real contributions. Using R ggg as an example, we find that the differences between algorithms are within the uncertainties of the calculation. This suggests that C (2) 3 is the same for the k T -type algorithms, and hence that the coefficient that we will extract is (at least somewhat) universal.  . Jet algorithm dependence in the calculation, shown for the ggg channel as an example. The fractional uncertainty for the anti-k T algorithm is shown in black, along with the percent difference in R ggg to the anti-k T result for the C/A (red, dashed) and k T (blue, dotted) algorithms. The fact that the difference between algorithms lies within the uncertainty of the anti-k T result suggests that the value of C is the same for the k T -type algorithms.
The coefficient a 2 in eq. (6.5) contributes to C (2) 3 , while the other coefficients are not used 7 . The entire fit procedure is performed over a grid of R c = {0.2, 0.5, 1.0} and x c = {0.1, 0.3, 1.0}. The variation in R c checks that the C (2) 3 contribution is independent of its value, while the variation in x c must cancel with the integrated subtractions. Indeed, we find this to be the case, and the final C (2) 3 value (adding the real and virtual results together) for the 9 fits over different R c and x c values are very consistent with each other 8 .

Results and Impact on the H + 0-jet Cross Section
With the results for the regulated real and virtual contributions in hand, we can combine them to determine C (2) 3 . We add the analytically known virtual terms from eq. (6.4) to the fits to the real terms and obtain 9 independent evaluations of C 3 is taken as the statistical average of these 9 determinations. In fig. 4, we show the value of C for the ggg and gqq channels for each fit as well as the averaged result.
We can evaluate the ggg real + gg virtual and gqq + qq virtual channels separately, and further divide the gqq + qq results into C F C A n f and C 2 A n f color channels. The results from each channel, and the total contribution, are 3,gqq,C F = 0.1405 ± 0.0011 , C 3,gqq,C A = −0.6854 ± 0.0091 , C 3,ggg + C 3,gqq,C F + C     Although we have found that at O(α 3 s ) the leading ln 2 R clustering logarithms are small, it would be interesting to determine the complete O(α 3 s ) clustering corrections through a H + 2-jet NLO calculation. Such a calculation is most powerful if it not only determines the numerical size of the correction for phenomenological ranges of parameters but also extracts the O(α 3 s ) rapidity anomalous dimension contributions. This would be an important part of extending the resummation for H + 0 jets to N 3 LL, as the anomalous dimensions connected to the global veto contributions may be simpler to determine (due to the lack of a clustering algorithm).

Conclusions
In this work we have calculated the leading O(α 3 s ) clustering logarithms for the jet vetoed cross section. Clustering effects from the jet algorithm start at O(α 2 s ), and those terms are numerically important. We have calculated the most important terms in the O(α 3 s ) correction, those with the form α 3 s ln 2 R ln m H /p cut T , and find that they are numerically less important than the O(α 2 s ) terms. This brings the higher order clustering effects under better control, and is the first step in determining the complete clustering effects that take part in the N 3 LL resummation of the H + 0-jet cross section.
In addition to the improved description of jet vetoed cross sections, this work combines several techniques to perform the calculation. What is a naively N 3 LO calculation is reduced to an NLO calculation through the combined factorization properties of QCD and SCET. Subtraction techniques for NLO calculations are then used to perform the calculation. Our success suggests that these techniques may be useful in other SCET calculations, and from a pedagogical perspective it would be interesting to explore how SCET interfaces with fixedorder subtractions beyond NLO. In sec. 4 we have seen an example of how soft, collinear, and collinear-soft subtractions at NLO are given by SCET matrix elements. This is expected, since subtractions and SCET matrix elements are based on the same singular limit of QCD amplitudes. At NNLO, the factorization properties and organization methods of SCET may help organize the larger set of subtraction terms. g → qq : The 1 → 3 gluon-initiated splitting functions and their polarization averages are g → ggg : We rescale the other subtractions terms similarly and denote them with hat, C and CS. eq. (B.1) implies that the integrated soft subtraction is, once we add in the phase space cut on the auxiliary x c , For the case of the 1r or 2r soft gluon exchange, where the form matches the collinear-soft subtractions, we also add a cut on R c . Similarly, the integrated collinear-soft subtraction is (B.4) For the collinear subtraction, the Born matrix element is expressed in terms of momenta k 1 + k 3 (in the collinear limit) and k 2 . This means the p T fraction variables must be written in terms of x 1 +x 3 and x 2 , while the angular variables are suitable as-is (the collinear splitting is parameterized by y 13 and φ 13 , the Born by y 12 and φ 12 ). We perform the change of variables The variable z parameterizes the collinear splitting. In terms of these variables, the collinear limit of the real matrix element is Adding the cut on R c , the integrated collinear subtraction is (B.7) Note that the Born matrix element is expressed in terms of the momenta k 1 + k 3 and k 2 . However, since we have written the phase space integration in terms of these momenta, they are really just dummy variables at this point, and this is why the collinear subtraction factors out cleanly.

B.1 Soft Subtractions
There are two types of soft integrated subtractions: those which exchange the soft gluon between two collimated partons (e.g., 1 and 2) and those which exchange the soft gluon between one collimated parton and the rest of the event in a color coherent way (e.g., 1 and r). In the first case, the integrated subtraction is The p T fraction x 3 can be integrated over directly, and the angular integrals are straightforward by using ∆R 2 23 = (y 13 − y 12 ) 2 + (φ 13 − φ 12 ) 2 , (B.9) and integrating over y 13 first (which is independent of ). The result is This result is basically identical to the integrated FKS soft subtraction taken in the limit ∆R 12 1. A similar soft function was calculated in Ref. [45]. The second type of soft subtraction matches the kinematic dependence of the collinearsoft subtraction. For the collinear-soft subtraction we will impose cuts on x c and R c , while for the soft subtraction we will only impose a cut on x c . The integrated soft subtraction is (B.11) The collinear-soft subtraction has an additional constraint θ(∆R 13 < R c ). The angular integrals are again straightforward, but there is a constant term depending on R c whose functional form is difficult to obtain. These integrals to O( ) are (with and without the R c constraint): and is plotted in fig. 5. This soft subtraction is therefore

B.2 Collinear Subtractions
The integrated collinear subtraction can be written from eqs. The angular integrals have already been carried out in the soft subtractions above, and the integral over z is straightforward. The result is In this case we have not included a symmetry factor for the g → gg splitting. Otherwise, the result is very similar to the FKS integrated collinear subtraction, save for the different constant. For q → qg splittings, the result is and for g → qq we have The above coefficients are the standard ones, The superscripts (g) and (q) in the integrated subtractions denote the C A and n f parts respectively.

B.3 Collinear-Soft Subtractions
The collinear-soft subtraction CS 31 , in eq. (4.20), is the same as the soft subtraction S 3 1r in eq. (4.14) with a different color factor. Thus the integrated collinear-soft subtraction is given by the soft results above with the R c constraint and the appropriate color factor,

C Calculation of the O(α 2 s ) Clustering Logarithm
In this appendix we give the calculation of the O(α 2 s ) clustering logarithm, which has been computed previously but whose elements are recycled in computing the sum of the virtual and integrated counterterm contributions to C The 2 is a constant that was previously computed in Ref. [10], and is equal to (C.5)