Power Corrections in the N-jettiness Subtraction Scheme

We discuss the leading-logarithmic power corrections in the $N$-jettiness subtraction scheme for higher-order perturbative QCD calculations. We compute the next-to-leading order power corrections for an arbitrary $N$-jet process, and we explicitly calculate the power correction through next-to-next-to-leading order for color-singlet production for both $q\bar{q}$ and $gg$ initiated processes. Our results are compact and simple to implement numerically. Including the leading power correction in the $N$-jettiness subtraction scheme substantially improves its numerical efficiency. We discuss what features of our techniques extend to processes containing final-state jets.


Introduction
Accurate predictions for high-energy scattering processes at hadron colliders rely upon calculations to higher orders in the perturbative expansion of QCD. Fully differential predictions are needed, in order to correctly model the final-state cuts imposed in all experimental analyses and directly compare theory with data. The calculation of higher-order corrections is complicated by the fact that the real-emission and virtual corrections which contribute to the cross section exhibit infrared singularities that cancel only after they are combined. At next-to-leading order (NLO) in the strong coupling constant several well-established techniques have been successfully applied for many years to accomplish this task [1][2][3]. In the past several years, new schemes [4][5][6][7][8][9][10][11] have been proposed that enable calculations through the next-to-next-to-leading order (NNLO) in perturbative QCD, and permit precision comparisons of theoretical predictions with data from the Large Hadron Collider (LHC) In this work we focus on the N -jettiness subtraction scheme for higher-order calculations [10,11]. This conceptually appealing idea uses the N -jettiness event shape variable τ N [12] as a resolution parameter to isolate and cancel the double-unresolved singular limits, where two partons become soft and/or collinear, that complicate the calculation of NNLO cross sections. The idea of using a physical observable to isolate singularities and construct subtraction terms at higher orders stems from the idea of q T -subtraction [6], which uses the transverse momentum of a color-singlet object to accomplish NNLO computations for color-singlet processes. This was further generalized to calculate top-quark decays [13] and tt production in leptonic collisions [14]. The N -jettiness subtraction scheme further extends this idea to handle processes with arbitrary final-state jets. When τ N is large, an N + 1 jet configuration is guaranteed, leading to a NLO contribution to the N -jet process that can be obtained with standard techniques. When τ N is small, the NNLO results can be obtained using an effective theory approach [15][16][17][18][19]. The N -jettiness subtraction scheme has proven quite successful in enabling NNLO phenomenology. It has led to some of the first calculations for vector boson production in association with a jet [10,[20][21][22][23] and Higgs production in association with a jet [24] at the LHC through NNLO. It has also led to first predictions for inclusive jet production at NNLO in electron-nucleon collisions [25], and has reproduced known results for color-singlet production through NNLO [11,[26][27][28] Color-singlet production through NNLO using the N -jettiness subtraction scheme has been publicly released in the numerical code MCFM v8.0 [28].
The N -jettiness subtraction scheme relies upon the introduction of a cutoff τ cut N that separates the N + 1 jet configuration from the doubly-unresolved limit. The below-cut region is expanded in τ cut N /Q, where Q denotes the hard momentum transfer in the process, in order to allow for an effective field theory calculation. The cutoff must be chosen small so that the power corrections in τ cut N /Q are negligible. However, the below-cut and above-cut contributions separately depend on logarithms of τ cut N /Q that only cancel after the two regions are combined. Since these regions live in different phase spaces and are numerically integrated separately, these logarithms introduce numerical noise that challenge the efficiency of the method. Although the numerics can already be controlled sufficiently for phenomenological applications, it is desirable for computational efficiency to reduce the sensitivity of the method to the power corrections. An explicit calculation of at least the leading power correction would allow N -jettiness subtraction to be used with larger τ cut N , reducing the computational cost of the approach.
In this work we discuss the analytic calculation of the dominant power corrections through NNLO 1 . In Section 2 we show in detail the derivation of the leading-logarithmic power correction at NLO for an arbitrary N -jet process. We summarize which features of this result generalize to the NNLO level. In Section 3 we explicitly calculate the leadinglogarithmic power correction at NNLO for color-singlet production mediated by both qq and gg initial states. In Section 4 we study the numerical impact of the power corrections on the N -jettiness subtraction scheme. Our main results for the power corrections at NNLO are summarized in Section 3, in the form of simple analytic expressions amenable to numerical implementation.

Derivation of the NLO power correction
In this section we briefly review the leading-power factorization theorem and present a detailed derivation of the leading-logarithmic power correction at NLO for an N -jet process. We emphasize the features of the calculation which extend to the NNLO level.

Review of the leading-power factorization theorem
In the N -jettiness subtraction scheme [10,11] for a generic collider process involving N jets in the final state, the N -jettiness event-shape observable [12] serves as the resolution parameter between the N + 1 jet configuration and the doublyunresolved limit. Here, the p k denote the four-momenta of final-state QCD partons. The n a,b are light-like vectors for the initial beam directions and the n i are light-like vectors denoting the directions of the final-state jets in the problem. The n i are determined by pre-clustering final-state radiation using a standard jet algorithm. The Q i are variables characterizing the hardness of the beam jets and final-state jets. The minimum in Eq. (2.1) defines the contribution of p k to τ N according to which direction p k is closest. The smallτ N cross section is derived using the all-order leading-power factorization theorem for the cross section [30], obtained using the Soft-Collinear-Effective Theory (SCET) [15][16][17][18][19]. We schematically write the differential cross section in the small-τ N limit as where the operator definition of each components can be found in Ref [30]. The beam function B [31,32], the jet function J i [33,34] and the soft function S N for jets [35] and for the massive case [36] are all known to the required NNLO level. The results of Eq. (2.2) expanded to fixed order in the strong coupling constant can also be obtained using the method of regions [37]. This entails expanding the full QCD matrix elements and the phase space consistently in τ N assuming either a soft-momentum scaling p s ∼ Qτ N or a collinear momentum scaling p c ∼ Q 1, √ τ N , τ N 2 , which exhaust the possible leading singular regions. In writing the collinear momentum scaling we have adopted the usual Sudakov decomposition of p c . In the method of regions approach, the collinear and soft behaviors are disentangled and the fixed order results of Eq. (2.2) can be recovered as the sum of soft and collinear contributions. To avoid double counting between the collinear and soft regions, a zero-bin subtraction [38] is usually required in order to reproduce the leading singular results of QCD. The leading-power factorization theorem is exact only when τ N → 0. Therefore in practical applications of N -jettiness subtraction a very small τ cut N is introduced, and Eq. (2.2) is used below τ cut N . The below-cut cross section receives power corrections with the leading behavior α n s log 2n−1 (τ N ). Including the power corrections would allow larger τ cut N to be used, and could in principle improve the numerical performance of the N -jettiness subtraction scheme.

Power corrections at NLO
In this section, we calculate the leading power correction at NLO to the N -jettiness factorization theorem in Eq. (2.2). We illustrate our derivation focusing on the τ 1 measurement in deep inelastic scattering (DIS). This general case contains both an initial-state hadron beam and a final-state jet, allowing us to demonstrate all possible technical details. We define N -jettiness using the hardness measures Q −1 a = x √ s/s ll and Q −1 1 = 2E J /s ll , where E J is the jet energy, although the derivation holds more generally. Although we show in detail the derivation of the leading-logarithmic terms (log τ 1 at NLO) in the leading-logarithmic power correction (LL P ), the remaining terms which scale as τ 1 can also be evaluated in a straightforward manner. We will present the NLO power correction for DIS, jet production in e + e − collision and Higgs production in gg fusion in the end of this section, which represent the possible cases for the N -jettiness NLO power correction. The NLO power corrections for a generic N -jet process can be obtained as a linear combination of these three cases with color factors assigned properly, plus corrections derived from the N -jet Born matrix element.

General features for LL P at NLO
Before showing the explicit calculation, we note that on general grounds, the NLO cross section for measuring τ N takes the form where z parameterizes the energy of the final-state radiation and τ N controls the collinear singularity. As z → 1, the emission becomes soft and forces τ N → 0, which justifies the upper limit in the z-integral.
Here the form of f (τ N ) depends on the kinematic details. f (τ N ) starts at O(τ N ) and vanishes as τ N → 0. The numerator N is a regular function in both limits z → 1 and τ N → 1, and includes information from both the matrix element and possible phase space cuts. If we expand N in terms of τ N and 1 − z we can see that N 0 will give rise to the contribution covered by the leading-power factorization theorem in Eq. (2.2). N 1 will generate the log(τ N ) power correction, while the remaining terms will contribute to terms linear in τ N or higher orders in powers of τ N . Therefore determining N 1 in which the emission becomes soft is our primary goal. We begin by listing general features of the LL P that we find at NLO that hold true at NNLO as well.
• The results are free of divergences. All -poles cancel among themselves in the power corrections. This is a clear requirement of the LL P ; the differential cross section in τ N is a physical observable free of divergences, and it can be expanded in the small-τ N limit to obtain the LL P .
• The LL P comes solely from the soft limit. The configurations which contribute to the leading-power expression in Eq. (2.2) but not to its leading logarithms, do not give rise to LL P . An example of a leading-power configuration that does not contribute to the LL P is in the gl → qql channel, when one of the q is grouped with the beam to contribute to τ .
• Soft quarks contribute to the LL P . Their contributions can be determined unambiguously from the leading power splitting kernels.
• Power divergences occur in the power corrections due to expanding terms such as The power divergence can be eliminated by rescaling The rescaling leads to the appearance of derivatives of the parton distribution functions (PDFs) ∂ x f (x) in the power correction.
We have used both a rigorous QCD calculation and the method of regions to obtain our results. We have checked that these two methods lead to identical expressions. In the following discussion we present the calculation using the method of regions, since we will later extend this approach to the NNLO level. A similar procedure has been applied in Ref. [39,40] with focus on the threshold power corrections in Drell-Yan production. We discuss the derivation using the DIS process as an example. In DIS, following the N -jettiness definition Eq. (2.1), the calculation can be organized by the beam and the jet contributions depending on whether the radiation is grouped with the beam or with the jet direction to contribute to the jettiness. We note that both the beam and jet contributions are well-defined, IR-safe physical observables. The global N -jettiness is the sum of these two contributions. In the method of regions, both beam and jet contributions receive their dominant part from the configurations in which the radiation momenta become collinear (scales as collinear) or soft (scales as soft) which will be defined later when we discuss the beam and jet contributions. In each contribution, beam or jet, the final results will be the sum of the collinear and the soft sectors, with proper zero-bin (overlap between collinear and soft) subtracted out to avoid double counting.
to the cross section have already been included Eq. (2.2). The phase space for this process can be written as Here, the measurement function Θ enforces that k is grouped with the hadron beam to contribute to τ 1 which defines the beam contribution. P a and p b are the four-momenta for the incoming hadron and lepton respectively. The Bjorken-x relates the hadron momentum P a to the partonic momentum as p a = x P a . l denotes the outgoing lepton 4-momentum, while q and k are the momenta for final-state partons. We always assume that k is potentially unresolvable while q is always hard. The last δ-function involving τ 1 defines the jettiness . We note that this definition makes τ 1 dimensionless. The superscript on the differential phase space denotes that this is an NLO expression.
To proceed, we parameterize the momentum k using the light-cone coordinates defined by p µ a and q µ : Note that in the beam contribution p a · k ∼ O(τ 1 Q 2 ) and k 2 ⊥ ∼ q · k p a · k as required by the jettiness definition and the on-shell condition for k. Here q · k can either be large (of order Q 2 ) or small (of order Q 2 τ ), which defines the collinear scaling and soft scaling of the momentum, respectively. However for the collinear scaling, when we perform the phase space integration, the momentum component q · k will unavoidably reach a region in which q · k ∼ Q 2 τ . In this region the momentum scaling overlaps with the soft scaling which defines the zero-bin. We therefore need to subtract out the zero-bin in the collinear sector to avoid double counting.
In the following, we detail the evaluation of the collinear scaling contribution q ·k ∼ Q 2 . The soft ones can be obtained similarly by assuming q · k ∼ Q 2 τ 1 . Writing k in terms of the light-cone decomposition, the δ-function for energy-momentum conservation in Eq. (2.5) can expressed as where we have introduced a variable z = s ll 2pa·q , with s ll = 2p b · l. The collinear scaling q · k ∼ Q 2 means 1 − z ∼ O(1). Here we have dropped in the δ-function the k ⊥ dependence. This is allowed to the logarithmic accuracy in which we are interested, since any term linear in k ⊥ will vanish after being averaged over the solid angle of k µ ⊥ , while any k 2 ⊥ will scale as τ 1 (1 − z ) which will not contribute to the logarithms in the power correction. 3 From the τ 1 definition in Eq. (2.5) and the momentum conservation expression in Eq. (2.7), it is straightforward to find To avoid the possible occurrence of power divergences in deriving the power correction as a consequence of expanding (1 − z + z τ 1 ) −1 in τ 1 , we have rescaled the variable z using This fixes the ambiguity of order τ 1 in defining O(1) variables. We note that the upper bound of the z-integration is 1 − τ 1 . However, we can safely integrate z all the way to 1.
The range [1 − τ 1 , 1] is where q · k ∼ Q 2 τ 1 and the zero-bin subtraction should be performed as we discussed before. After correctly implementing the zero-bin subtraction, this range receives no LL P . The phase space in Eq. (2.5) is factorized into a Born piece which only involves leading power momenta which scale as O(Q), and a radiative phase space: (2.10) The Born phase space is The radiative phase space is needed to obtain the Born phase space. Assuming the collinear scaling q · k ∼ Q 2 , p a · k ∼ Q 2 τ 1 and z ∼ 1, we find the phase space for the collinear sector in the method of regions approach where Ω 2−2 is the 2 − 2 dimensional solid angle from k ⊥ .
The phase space for the soft sector and the zero bin can be derived straightforwardly from Eq. (2.13) by assuming q · k ∼ Q 2 τ 1 and dropping the q · k dependence in the fist δ-function. The explicit form will not be shown here. We also note that the measurement function Θ can be set to 1 in the collinear sector up to the accuracy we are working with.
The contribution from expanding the PDFs in the beam sector is which receives no power correction. There is also contribution from the flux which again receives no power correction in τ 1 . Before we turn to the matrix element evaluation, we notice that we can make the simplification where we have omitted terms of order O(τ 2 ). We now consider the squared matrix elements for this process assuming the same collinear scaling as in the phase space. There are two partonic channels to consider, ql → l q g and gl → l qq. We first note that in the beam sector the gl channel does not contribute to the LL P , as discussed in our presentation of the general features of the NLO LL P . For the ql channel, using the simplification of scalar products above we find that the matrix element for the ql channel with a soft gluon reduces to , (2.18) and for ql channel with an unresolved quark we have the matrix element M (0) is the LO matrix element for ql → q l with leading-power kinematics. For the unresolved quark, the LL P comes from the configuration in which a quark is emitted from the final state q but grouped with the beam to contribute to τ 1 . Although this "anticollinear-grouping" configuration does not contribute to the leading-power singular terms in Eq. (2.2), it does have an effect in the power correction. We further note that to get Eq. (2.19), instead of expanding the full NLO matrix element for real emission, Eq. (2. 19) can be determined directly from s −1 gq P qg (x) by relating s qg ∼ (1 − z) and 1 − x ∼ τ 1 . Combining the phase space in Eq. (2.14) and the matrix elements in Eq. (2.18) and Eq. (2.19), we find the power correction from the collinear sector to be where from the effective theory point of view, √ τ s ll fixes the collinear scale, as expected.
Following similar procedure by assuming q · k ∼ τ Q 2 , we find the contribution from the soft scaling is Once we combine these, we arrive at the following logarithmic power correction from the beam contribution: The ratio in the logarithm reflects the scale hierarchy between the collinear and soft sectors. In the power corrections, no remaining singularities in should arise. All the poles must cancel amongst soft and collinear regions since the beam contribution itself is a well defined physical observable. We indeed find that they do.

Jet contribution
We next study the jet contribution following the same steps as for the beam region. Momentum conservation for xP a + p b → l + q + k can be simplified to Unlike the beam contribution, we have defined z = 2pa·q s ll . The one-jettiness definition becomes Deviations from this result due to pre-clustering the partons with different infra-red safe jet algorithms will not contribute to the LL P . We can parameterize where as required by the jettiness measurement, k · q ∼ Q 2 τ 1 while k · p a can either be of order Q 2 or Q 2 τ 1 , which defines the collinear scaling and soft scaling, respectively. In the following, we will focus on the collinear sector in which k · p a ∼ Q 2 . This also implies that 1 − z ∼ O(1). Here we have rescaled z using to avoid power divergence from expanding the matrix element in τ 1 . Similar to the beam contribution, the phase space is factorized into a Born part with and a radiation piece The first factor is again the Jacobian J = (1 + τ 1 )z d−3 from the variable change {x, q} → {x B , q B } needed to reach the Born phase space. Assuming that k µ follows the collinear scaling p a · k ∼ Q 2 , q · k ∼ Q 2 τ 1 and z ∼ 1, we find the phase space for the collinear sector in the method of regions: (2.30) Deriving the phase space for the soft and the zero-bin subtraction follows similar steps but with the additional assumptions that p a · k ∼ Q 2 τ 1 . The power correction due to expanding the PDFs around x B is given by while the flux contributes as Having derived the power correction coming from the phase space, we move onto the matrix elements. Noting the simplification in which all terms linear in k ⊥ or proportional to τ 1 (1 − z) have been dropped, we find that the matrix element for the ql channel with an unresolved gluon reduces to , (2.34) and the matrix element for the gl channel with an unresolved quark simplifies to This second matrix element comes from a soft quark emitted from the initial state p a but grouped with q to contribute to τ 1 . Again M 0 is the LO matrix element for ql → q l with leading-power kinematics.
The contribution from soft sector can be obtained in a similar manner with the assumption that p a · k ∼ Q 2 τ 1 . Putting together all components, we find the power correction from jet contribution is Again, the -poles are absent in the final result, since the jet contribution is itself a physical observable. Here the PDFs are evaluated at x = x B . In the gl channel, we have normalized the result to the ql channel color and spin average by multiplying by a factor of N C N 2 C −1 . dσ 0 is the Born-level matrix element and phase space with the PDF removed.

NLO power correction for DIS
Combining the contributions from the beam and jet contributions found in the previous sections, we obtain the full power correction the DIS τ 1 distribution as where Q i is the electric charge of the ith quark, L = − log(τ 1 ) and the PDFs are evaluated at x = x B . We note that the power correction for DIS comes from the Jacobian J in Eq. (2.14), from expanding PDFs in Eq. (2.31), and from the soft quark configuration.

NLO power correction for hadronic production in e + e − collisions
The power correction for N -jettiness (or equivalently thrust) for hadronic production in e + e − collisions is found to be dσ (1) e + e − →hadrons = dσ 0 where L = − log(τ ). The derivation follows closely the one presented for the DIS beam contribution. The result reproduces the known logarithmic power correction found from both fixed-order calculations [41] and within the framework of SCET [42], which demonstrates the validity of the method of regions approach in studying the power corrections. Here the power correction comes from Jacobians due to rescaling the Born momentum q 1 or q 2 to q 1,B or q 2,B in the process l + l − → q(q 1 )q(q 2 )g(k), and from the soft-quark matrix element.

NLO power correction for ggH and Drell-Yan
Following a similar procedure as in the previous sections, we can evaluate the power corrections for ggH. We define τ 0 using the hardness measures Q a = Q b = 1 in Eq. (2.1), implying that τ 0 has units of energy. For the gg channel, we find and L gg = L gg (x 1 , x 2 ) is the gluon-gluon luminosity. We note that dσ 0 is the Born-level differential cross section for gg → H with the PDFs removed, and Y is the rapidity of the Higgs. For the qg + gq channel, we have Here, L follows the definition above and L g 1 q 2 is the gluon-quark luminosity. For Drell-Yan production of lepton pairs through a vector boson V , the q iqj channel power correction gives while the qg + gq channels contribute where Q j is the change carried by the final-state soft quark j. V ji = δ ji for Z production and is the CKM matrix for W production. Since we are not measuring the final state flavors a sum over j arises. dσ 0 is again the Born-level cross section for qq → V with the PDFs removed. One should also replace m H → m V in the definition of the logarithm L for the Drell-Yan case. In both ggH and Drell-Yan cases, the net NLO power correction comes from expanding the PDFs and from the soft quark matrix element.

Derivation of the NNLO power correction
We next present the derivation of the power corrections at NNLO. At O(α 2 s ), we have to deal with both real-virtual (RV) and real-real (RR) contributions. The double-virtual correction has been entirely included in the leading power factorization theorem in Eq. (2.2). For RV, the phase space integration is identical to the NLO phase space, For RR, the calculation is more involved. We begin by stating some of the features of the LL P that we observe at the NNLO level.
• All -poles cancel between RV and RR in the LL P .
• Soft limits lead to LL P at NNLO just like at NLO. The LL P soft currents for the qg channels for both Higgs and Drell-Yan production are deducible from the leadingpower splitting kernels.
• Explicit calculations show that the LL P comes from the strongly-ordered limit of the matrix elements at NNLO. For instance, in the case of a qg final state, E g E q . We also notice from our explicit calculations that within the strongly-ordered limit, the final results for the LL P can be obtained using an independent-emission approximation in the phase space integrals.
• Configurations (the Abelian piece) which contribute to the leading logarithms in the leading-power result also contribute to the LL P . Their contribution to the LL P at NNLO can be written as a convolution in τ N between an NLO leading-power leading logarithmic contribution and an NLO LL P .
Although we can not prove or disprove on general grounds these observations seen in our explicit calculations, we conjecture that all of the above features observed for zero-jettiness in color-singlet production generalize to a generic N -jet case. At NNLO, the full calculation is more lengthy than at NLO. Here we sketch the procedure of obtaining the final result, making sure to discuss all relevant features. We consider the following examples from the Drell-Yan process to illustrate our method: the RV corrections coming from the qg → V q channel, and the RR corrections to the qg → V qg channel.

Sketch of the RV calculation
We start with RV. The construction of the phase space follows identically the procedure at NLO. The only complication is that the matrix elements no longer have a Taylor series in τ 0 . They possess fractional powers coming from the loop integrations that must be isolated in order to obtain the correct logarithms, since the logarithms come from the expansion of τ −n 0 factors hitting 1/ poles, where n is an integer. This can be done because the softquark limit contributing to the LL P comes completely from the "anti-collinear" grouping of the quark, just as at NLO and as discussed below Eq. (2.19). The factorization of the RV matrix elements in this collinear limit are well known [43], and contain the necessary decomposition into fractional powers. We then use the factorization of the matrix element in the collinear limit to write the matrix element as a sum of and In both M RV,C−loop , we have kept only terms that will contribute to the LL P . We have assumed that the quark is emitted from leg P b for P a P b → V + X and is grouped with P a in order to contribute to the LL P . Here m is the virtuality of the vector boson V and we use the notation that k − = n b · k and τ = n a · k = k + . We note that the collinear factorization of the one-loop RV amplitude schematically contains two pieces: a tree-level splitting amplitude times a one-loop amplitude for the process qq → V , and a oneloop splitting amplitude times a tree amplitude for qq → V . The M RV,C−loop comes from the second structure. Now following the same procedure for the NLO calculations, we can get the RV contribution to the LL P straightforwardly. The RV contributions are given by the sum of and I (2) Here we have normalized the result to the qq channel color and spin average by multiplying with a factor N C N 2 C −1 , which turns an overall factor C F to T R . We have suppressed the dependence on the luminosity and Born cross section for simplicity. The full RV contribution to the LL P cross section for final state qg is thus given by dσ (2) RV,qg+gq→V = dσ 0 I (3.5)

Sketch of the RR calculation
For the RR contribution we find that the matrix elements leading to the LL P can be written as a sum of Eq. (A.1) to Eqs. (A.5), which can be derived using the soft limit of the s −2 qgg P qgg splitting kernel [44]. To evaluate the cross sections, both a rigorous QCD phase space integration and the method of regions are used to derive the final results. For RR we need to consider collinear-collinear, soft-soft and collinear-soft scalings in the method of regions, together with suitable subtraction of zero-bins. We find that the final results can be obtained using an independent emission approximation in the phase space integrals following a strongly ordered limit in which E g E q . The relevant strongly-ordered soft current can be found in Eq. (A.6). We organize our calculation by partitions in which both qg are grouped with n a,b to contribute to τ , and qg grouped separately with n a or n b to contribute. We sum all the partitions to find the full contribution. Following these steps leads us to the final results. As an example of the integration of one of the contributions in the Appendix over the relevant phase space Φ[k 1 , k 2 ] of k 1 and k 2 , we find for Eq. (A.2) for q(x a P a )g( where the ellipsis denotes omitted terms which will not contribute to the LL P . We have normalized the result to the qq channel color and spin average which turns a factor of C F to T R , and have suppressed the dependence on the Born cross section and the luminos- RR,qg+gq→V = dσ 0 I Here, the I ij are presented for completeness in Eq. (A.10).

NNLO power corrections for ggH and Drell-Yan
We now summarize our final expressions for the LL P power corrections at NNLO. For the qq channel in Drell-Yan we have the power correction which is accurate to the LL P . We note that the results can be obtained by a convolution in τ 0 between an NLO leading-logarithmic leading-power term (given in Eq. (B.1) for completeness) with an NLO leading-logarithmic term at subleading power, given in Eq. (2.41). We note that these expressions contain sub-leading logarithms associated with the scale dependence that come "for free" as a result of our derivation. We keep these in our final results, although a strict expansion keeping only log 3 (τ 0 ) structures is possible. For the qg + gq channel, the RR contribution to the LL P can be derived using the soft currents in Eq. (A.1) to (A.5). Though each term gives a somewhat lengthy expression, when summed up the final results for the power correction are simple and compact. We have for the full result, after combining RV and RR: Here we have again showed the explicit dependence on the quark charge Q j , and the CKM matrix V ji in the case of W -boson production (this can be set to a Kronecker delta in the case of Z-boson or γ * production). We note that the results for the RR contribution can be obtained using Eq. (A.6), which is the strongly ordered limit of Eq. (A.1) to (A.5) in which E g E q . In the strongly-ordered limit, the matrix element factorizes into a product of a sub-leading soft quark current and a known leading-power soft-gluon current involving three eikonal directions. We note that dσ 0 is the Born-level differential cross section for the qq → V process with the PDFs extracted.
For gluon-fusion Higgs production, the LL P contributions can be obtained in the same way as for Drell-Yan. For the Abelian case, it is found that the contributions can again be obtained by convoluting the NLO leading-logarithmic terms at leading power with the NLO LL P . For the qg final state, the RV can be extracted from, for instance, the 1-loop correction to s −1 qq P g→qq [43] with the proper crossings. The RR can be derived by suitably changing the color factors in RR for Drell-Yan. After performing the phase space integrations, we find that all the -poles cancel as required, and the LL P for the ggH are given by: for the gg channel at the LL P accuracy and dσ (2) for the gq + qg channel. Here we have normalized to the gg channel color and spin average.

Numerical results
We now study the numerical consequences of the derived power corrections in the Njettiness subtraction formalism. We focus on the Drell-Yan and gluon-fusion Higgs production channels at the LHC. The NNLO corrections to the Drell-Yan and ggH processes have been implemented in MCFM v8.0 using N -jettiness subtraction [28] and a thorough study of the impact due to the missing power corrections by comparing with the known exact results [45,46] has also been presented therein. In this section we follow exactly the same settings as used in Ref. [28] for H/Z/W production at a 13 TeV LHC with NNLO the MSTW2008 PDF set [47]. To study the impact of the power corrections calculated in the previous section, we generate the O(α 2 s ) NNLO coefficient using the N -jettiness subtraction scheme with and without the power corrections.
We begin with a numerical validation of the calculated power corrections using Wboson and Z-boson production at a 13 TeV LHC with scale choice µ R = 2m V and µ F = m V /2. This scale choice is made to increase the size of the NNLO coefficient. Due to the simple structure of the Drell-Yan cross section, the power corrections can be fitted to high accuracy by generating NNLO results with different τ cut N values, as has been first performed in Ref. [28]. In Fig. 1, we compare the calculated LL P (red solid line with dots) at O(α 2 s ) with the fitted results (green dashed line) for both Z-boson and W + -boson production. We consider inclusive Z-boson production, and impose the following final-state cuts in the case of W + production: y(lep) < 2.5, E T > 30 GeV. (4.1) From Fig. 1, we see that the calculated LL P and the fitted power corrections agree very well, and converge to each other in the small-τ 0 region. The discrepancy seen in the larger τ 0 region is due to the sub-leading logarithms which are included in the fit but not in the LL P . Now we turn to the predictions for the NNLO coefficients for inclusive Z/W production in which the impact of the power corrections are found to be relatively large [28]. In Fig. 2, we compare the N -jettiness subtraction scheme with and without the power corrections for both Z-boson production in the upper panel and W + -boson production in the lower panel. We have plotted the difference between the O(α 2 s ) coefficient computed using Njettiness and the exact O(α 2 s ) coefficient [45] and have normalized this to the known result. The vertical axis in fig. 2 characterizes the deviation from the exact NNLO correction 4 . The blue solid line with squared dots represents the version without the power corrections generated using the current MCFM v8.0 [28]. The red line with round dots shows the results when adding on the analytically-calculated power corrections, while the green dashed line is the result obtained using the exact O(α 2 s ) coefficient subtracting out the power corrections fitted numerically. From Fig. 2, we can see dramatic improvements in the convergence of N -jettiness subtraction when the LL P is added. Without the power corrections, τ cut 0 should be set to below 10 −3 GeV to reproduce the exact NNLO coefficients. The cut can be relaxed by a factor of 10 when the power corrections are included. As a consequence of the larger allowed τ 0 cut, the computing time and the numerical efficiency are greatly improved.
In Fig. 3, we show the comparison between N -jettiness subtraction with and without power corrections for Z-boson and W + -boson production with cuts on the lepton rapidity, and on the missing energy in the case of the W -boson. With the presence of the cuts, the convergence is already better than in the inclusive case. When the LL P power corrections are included, we see again a substantial improvement in the convergence of the N -jettiness subtraction scheme. At least a factor of 10 in increasing τ cut 0 is observed, which leads to more efficient realization of the NNLO calculation.
Finally, in Fig. 4, we consider gluon-fusion Higgs production at the LHC with scale choice µ R = µ F = m H . Again we show the O(α 2 s ) coefficients both with power corrections included (blue solid line with squared dots) and without them(red solid line with round dots). We have normalized the results to the known exact NNLO coefficient [46]. In the ggH case, the coefficients predicted without the power corrections already converge faster than the Drell-Yan case. Once the power corrections are included, similar improvement as in the Drell-Yan case is also observed for gluon-fusion Higgs production.

Summary
In this manuscript we have studied the leading-logarithmic power corrections in the Njettiness subtraction scheme. We have derived in detail the NLO power corrections for an arbitrary N -jet process, and have presented the important features of the derivation of the leading-logarithmic power corrections at NNLO for color-singlet production from both qq and gg initiated processes. The final expressions for the NNLO power corrections can be found in Section 3. We have found that for color-singlet production in hadronic collisions, the LL P at NNLO comes completely from the strongly-ordered soft limits. To get the LL P at NNLO, we only need the information from the NLO LL P , the leading-power collinear   Figure 3. The difference between the NNLO coefficients for Z-boson and W + -boson production with lepton and missing energy cuts at the LHC obtained using N -jettiness subtraction with and without power corrections, normalized to the exact NNLO coefficients. splitting kernels and the LO matrix elements. More interestingly the Abelian piece of the NNLO LL P is given by the convolution of the NLO leading logarithms in the leading power and the NLO LL P . We note that a similar convolution structure also exists in the threshold case when the Drell-Yan threshold results are carefully studied [45]. We conjecture that these features hold for the LL P for the production of an arbitrary number of jets.
The final results for the LL P , which are of the form α s log(τ ) at NLO and α 2 s log 3 (τ ) at NNLO, are compact and can be easily added to existing numerical implementations of the N -jettiness subtraction scheme. Once the LL P are included, the resolution parameter τ cut N can be relaxed, substantially improving the numerical convergence of the approach. We have demonstrated these improvements by studying both vector boson and Higgs production at the LHC in a variety of settings. In all cases the incorporation of the LL P allows the value of τ cut N to be increased by nearly an order of magnitude while maintaining the same level of agreement with known results for color-singlet production at NNLO.
In the future, it will be important to complete the derivation of the NNLO power corrections for an arbitrary number of jets following the approach outlined here, as well as obtain the power corrections beyond LL P to further improve the efficiency of the N -jettiness subtraction scheme. Other interesting directions included predicting and resumming power corrections within the framework of SCET [48,49].

Acknowledgments
We thank A. Isgrò for helpful conversations. R. B. is supported by the DOE contract DE-AC02-06CH11357. F. P. is supported by the DOE grants DE-FG02-91ER40684 and DE-AC02-06CH11357. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357. This research was also supported in part by the NSF under Grant No. NSF PHY11-25915 to the Kavli Institute of Theoretical Physics in Santa Barbara, which we thank for hospitality during the completion of this manuscript.
Note added: While we were finalizing this manuscript, Ref. [50] appeared, which calculates the NNLO leading-logarthmic power correction for color-singlet production through the qq partonic process. Although their approach is different from ours, their results for Drell-Yan production fully agree with the LL P expressions presented here.

A The soft current for the qg channel
We list here the soft current for the qg final state in Drell-Yan used to derive the LL P at NNLO. The q 1 (k 1 )g 2 (k 2 ) soft current can be derived by studying the leading-power NNLO splitting function [44] and is given by the sum of the following five terms: Hereŝ = m 2 ,t i = p a · k i andû i = p b · k i . The leading singular contribution in the strongly-ordered limit E g E q is given by We note that S s.o. factorizes into the product of a leading-power soft gluon current that knows about three eikonal directions, and a sub-leading soft quark current.
The final results obtained after integrating Eq. (A.1) through Eq. (A.5) over k 1 and k 2 can be identified with the terms in Eq. (A.6) in a one-to-one manner. Performing the phase space integration over S (2) 1 , S 2 to S (2) 5 , gathering all the pieces and splitting into different color factors according to Eq. (A.6), we find in the qg final state for Drell-Yan, the LL P receives contributions from a 2C F − C A term: Contributions from C A -correlated terms (from matrix element contributions proportional to (k 1 · k 2 ) −1 ) are given by + . . . , (A.9) while the contributions from C A non-correlated terms is The ellipsis in these equations denote terms that do not contribute to the final LL P and that have therefore been omitted. We have normalized the results to the qq → V color and spin average by multiplying with N C N 2 C −1 , which turns an overall factor of C F to T R . We have obtained I

B Leading-logarithmic terms at leading power
The leading-logarithmic coefficients at leading power for the Drell-Yan process at the NLO are given by where the dependence on the Born cross section has been suppressed. We note that the convolution gives where we have only kept the leading-logarithmic contribution, as denoted by the ellipsis.