Transverse-momentum spectra of electroweak bosons near threshold at NNLO

We obtain the next-to-next-to-leading order corrections to transverse-momentum spectra of W, Z and Higgs bosons near the partonic threshold. In the threshold limit, the electroweak boson recoils against a low-mass jet and all radiation is either soft, or collinear to the jet or the beam directions. We extract the virtual corrections from known results for the relevant two-loop four-point amplitudes and combine them with the soft and collinear two-loop functions as defined in Soft-Collinear Effective Theory. We have implemented these results in a public code PeTeR and present numerical results for the threshold resummed cross section of W and Z bosons at next-to-next-to-next-to-leading logarithmic accuracy, matched to next-to-leading fixed-order perturbation theory. The two-loop corrections lead to a moderate increase in the cross section and reduce the scale uncertainty by about a factor of two. The corrections are significantly larger for Higgs production.


Introduction
It is often believed that the virtual corrections are the most complicated piece of higherorder perturbative computations, but in general this is not true. Beyond next-to-leading order (NLO), the real-emission corrections typically present the main difficulty. For example, because of the complicated structure of the real emissions, the NNLO corrections to W , Z and photon production in association with an energetic jet are currently still unknown, despite the fact that the two-loop virtual corrections to these processes have been available for more than ten years [1,2]. Over the past years, a number of techniques have been developed to isolate the soft and collinear singularities at NNLO so that the remaining phase-space integrations can be performed numerically, see e.g. [3][4][5]. Using these methods, the NNLO corrections for dijet production [6] and for Higgs-boson production in JHEP02(2014)004 association with a jet [7] have recently been obtained for the dominant, purely gluonic, partonic channel. While these computations are very challenging, we can expect that NNLO results for all electroweak bosons will become available in the not too distant future.
Soft and collinear singularities complicate the numerical evaluation of phase-space integrals, but the amplitudes themselves greatly simplify in the singular regions. In the present paper we will use this simplification to obtain an approximate NNLO result for Higgs-boson, W , Z and photon production in association with an energetic jet. To do so, we will consider the production of an electroweak boson V near the partonic threshold. Near threshold, the electroweak boson recoils against a low-mass hadronic jet, and the real radiation is either soft, or collinear to the jet or the incoming partons. As a consequence, the partonic cross section in the channel a + b → V + j c factorizeŝ ab is the Born cross section and the partonic Mandelstam variables areŝ = (p a + p b ) 2 ,t = (p a − q) 2 andû = (p b − q) 2 , with q the vector-boson momentum, and q 2 = M 2 V . The jet functions J q and J g describe the collinear radiation initiated by an energetic quark or gluon, respectively. The two-loop quark jet function was computed in [8], and the twoloop gluon result was obtained in [9]. Last year, also the two-loop soft function S ab was computed [10]. Here, we will determine the final NNLO ingredient, the two-loop hard functionĤ ab , for all partonic channels for both vector-boson and Higgs-boson production. The hard function contains the virtual corrections which were computed in [1,2,[11][12][13]. We will convert these results into an infrared-subtraction scheme that is compatible with the jet and soft function calculations. Our paper completes the construction of the twoloop cross section near the partonic threshold. Since the threshold terms usually amount to the bulk of the cross section, we expect that our result is a good approximation to the full NNLO result.
The derivation of the factorization formula (1.1) in Soft-Collinear Effective Theory (SCET) [14][15][16] was given in [17] for photon production and generalized to the W and Z case in [18,19]. An interesting new feature, which first arises at NNLO is the presence of a purely gluonic channel gg → V g in Z-boson and photon production. This was already pointed out in [17], but here we explicitly give the relevant SCET operators and determine their Wilson coefficients using the results [12] for the corresponding loop amplitudes. Numerically, however, we find that this channel only gives a negligible contribution. The factorization formula (1.1) can also be used to resum the threshold terms to all orders. For W and Z production at large transverse momentum, this resummation was performed at next-toleading logarithmic (NLL) accuracy in [20][21][22] and to NNLL in [17,18,23,24]. With all two-loop ingredients in place, we can extent the resummation to N 3 LL accuracy since the necessary anomalous dimensions are known. We have implemented the N 3 LL resummation, as well as the NNLO fixed-order expansion of the threshold cross section into a public code PeTeR ("large-p T Resummation") [25]. Numerically, we find that the resummation effects are small and the N 3 LL resummed results are close to the NNLO threshold results.

JHEP02(2014)004
Our paper is organized as follows. In section 2 we show how to extract the Wilson coefficients of the SCET operators describing V + jet production from the results for fourpoint amplitudes in the literature. The construction of the associated hard functions is then detailed in section 3. We first treat the vector-boson case, where we discuss the qq → V g and qg → V q channels in section 3.2 and the gg → V g process in section 3.3. After this, we construct the hard functions for Higgs production in section 3.4. With all the ingredients in place, we then compute the two-loop threshold cross section in section 4 and study the numerical impact of the corrections in section 5.

On-shell matching and renormalization
The hard functions in SCET are obtained by performing a matching calculation, i.e. by computing the same quantity in QCD and in the effective theory, and then fixing the Wilson coefficients of the SCET operators in such a way that the QCD result is reproduced. The Wilson coefficients are independent of the process and the external states used to perform the matching. By far the simplest possibility is to use on-shell amplitudes, in our case qq → V g, qg → V q and gg → V g, since the loop integrals in the effective theory are scaleless for on-shell momenta and vanish in dimensional regularization. The relevant twoloop QCD amplitudes were obtained in [2,[11][12][13]. We now explain how the SCET operators are constructed and how their bare Wilson coefficients are obtained from QCD results for the on-shell amplitudes, and then perform the renormalization of these coefficients.
The SCET operators mediating the production of an electroweak boson at large transverse momentum p T involve collinear fields associated with the two beam directions and the direction of the associated jet. At leading power in the effective theory, they involve a single collinear field for each of the three directions. To construct the collinear Lagrangians, one introduces a light-cone reference vector for each direction. The vectors n 1 and n 2 point in the beam directions, while n J is along the jet direction. Each reference vector n i has a conjugate light-cone vectorn i , with n i ·n i = 2. Quarks collinear to the direction i are described by a field χ i , which fulfills the condition n i / χ i = 0, so that this field is effectively a two-component field. Also, at leading power, only the components of the gluon field A ν,a i⊥ transverse to its light-cone direction can contribute. Because of these conditions, there is a one-to-one correspondence between helicity states and associated operators. Helicity amplitudes are therefore particularly well suited to extract the SCET matching coefficients, as stressed in [26][27][28].
We will now explain the relation between helicity amplitudes and Wilson coefficients in detail, using the example of the purely gluonic channel which arises for Z and γ production at one-loop and contributes at NNLO to the cross section. The presence of SCET operators mediating the gg → V g process was pointed out in [17], but in contrast to the operators for qq → V g and qg → V q, the purely gluonic operators were not explicitly given since they are not needed at NNLL accuracy. According to our considerations from above, the leading-power operators for gg → V g (and for gg → Hg) have the form

JHEP02(2014)004
The operators for the quark channels such as qq → V g have the same structure, but involve a quark field, an anti-quark field and a gluon. In SCET the collinear operators are smeared over the directions associated with the large external momenta, and the associated hadronic vector current J µ (x) mediating gg → V g consists of a convolution of the Wilson coefficients with the smeared operators Because we have left the color and Lorentz structure of the fields in the operators open, the Wilson coefficients have Lorentz and color indices, too. These are contracted with the operators to ensure that J µ (x) transforms as a singlet under color and a vector (or axialvector) under the Lorentz group. The color structure of the Wilson coefficients can be either symmetric or antisymmetric. In the first case, it is proportional to d abc , in the second case proportional to the structure constants f abc . Working with open Lorentz and color indices is convenient because the coefficients C abc µνρσ are directly related to helicity amplitudes in color space. To see this, we perform the matching using on-shell gg → V g amplitudes. Since the loop integrals in the effective theory are all scaleless for on-shell external momenta, the effective theory amplitudes reduce to tree-level matrix elements multiplied by the Wilson coefficients. Furthermore, because the different collinear sectors no longer interact after soft-collinear decoupling, the matrix element factorizes into individual collinear matrix elements, which in a given sector have the form for an incoming transverse gluon field. Performing the integrations over the variables t i , one then finds that the Fourier transforms of the Wilson coefficients C abc µνρσ (t 1 , t 2 , t J ), contracted with the external polarization vectors, are equal to the helicity amplitudes.
The vanishing of the loop corrections in the effective theory implies that in the relevant integrals, the infrared (IR) and ultraviolet (UV) singularities exactly cancel each other. Since the IR singularities of QCD and SCET are the same, this further implies that the UV singularities of SCET Wilson coefficients are identical to the IR singularities of QCD amplitudes. As a consequence the IR singularities of n-point amplitudes in d = 4 − 2ǫ dimensions can be renormalized multiplicatively [26,27] |M ren ({p}, µ) = lim where the renormalization factor Z is a matrix in color space. It is spin-independent, but depends logarithmically on the external momenta {p} ≡ p 1 , . . . , p n . The renormalized amplitude |M ren ({p}, µ) is equal to the renormalized Wilson coefficient of the leading-power SCET operator with the same quantum numbers as the external states in the amplitude. The inverse of the Z-factor can be written in terms of the anomalous dimension matrix [26] Γ({p}, µ) =

JHEP02(2014)004
where the sum runs over unordered tuples (i, j) of distinct partons, T i is the color generator associated with the i-th parton in the scattering amplitude, which acts on the color index of that parton, and s ij ≡ 2σ ij p i · p j + i0, where the sign factor σ ij = +1 if the momenta p i and p j are both incoming or outgoing, and σ ij = −1 otherwise. The product T i · T j ≡ T a i T a j is summed over a. Generators associated with different particles trivially commute, T i · T j = T j · T i for i = j, and T 2 i = C i is given in terms of the quadratic Casimir operator of the corresponding color representation, i.e. C q = Cq = C F and C g = C A . For more details on the color-space formalism, see [29,30].
In [2,[11][12][13] the results were given in terms of finite helicity amplitudes obtained after removing the IR singularities using Catani's subtraction formula [31]. In the following, we will relate these expressions to the renormalized SCET Wilson coefficients. The entire procedure can be viewed as a scheme change from Catani's subtraction scheme to a standard MS subtraction of the singularities.

Conversion to MS scheme
We first reconstruct the IR-divergent part of the two-loop amplitudes and will then perform the renormalization. We write the expansion of the UV-renormalized, on-shell n-parton scattering amplitude with IR singularities regularized in d = 4 − 2ǫ dimensions as where α s ≡ α s (µ) is the renormalized coupling constant. Note that the superscript (i) refers in this section to an expansion in units of α s /2π, which is the notation adopted in the literature on two-loop four-point functions [2,[11][12][13]. In the SCET literature, the perturbative expansion is usually written in α s /4π. Throughout this section, we will expand in α s /2π to be compatible with the literature on the amplitudes, but we will switch to the standard SCET notation when we present our result for the cross section in section 4 and in the appendices. The helicity amplitudes in [2,[11][12][13] were constructed using Catani's IR-subtraction formula [31], which states that the product is free of IR poles through O(α 2 s ). The amplitudes are however different from the MS-renormalized amplitudes |M ren ({p}, µ) in (2.4), because the subtraction operators I (n) (ǫ) ≡ I (n) (ǫ, {p}, µ) contain terms of arbitrarily high orders in ǫ. The explicit form of the I (n) (ǫ) can be found in appendix A. The above relation can be inverted to reconstruct the expansion coefficients of the IR-divergent amplitude |M(ǫ, {p}) as

JHEP02(2014)004
The SCET Wilson coefficient is now obtained by multiplying the IR-divergent amplitude with the inverse of the Z-factor. With a slight abuse of notation, we write the expansion of the inverse Z-factor in the form The explicit form of the coefficients Z (n) (ǫ) is given in appendix A. The above relations can be used to express the MS-renormalized amplitude |M ren ({p}, µ) in terms of the IR-finite amplitude |M fin ({p}, µ) given in [2,[11][12][13]. At one-loop order, the conversion relation reads where C 0 is the finite term of Catani's one-loop subtraction operator I (1) (ǫ), At two-loop order, the conversion relation takes the form 12) and the corresponding expression for C 1 and the two-loop anomalous dimensions are summarized in appendix A. In the appendix, we also give an explicit formula for the commutator [Γ 0 , C 1 ] in terms of three-particle correlations. The above relations are valid for general n-parton scattering amplitudes. For n = 3 colored partons, which is the relevant case here, one can use color conservation to express the dipoles in terms of the Casimir operators associated with the three external legs, The color structure then becomes trivial, and the one-loop conversion factor (2.11) simplifies to (2.14)

JHEP02(2014)004
The two-loop relation (2.12) contains in addition the structure as well as a commutator [Γ 0 , C 1 ], which vanishes in the three-parton case because of the trivial color structure.

Hard functions from helicity amplitudes
The two-loop four-point helicity amplitudes with one external electroweak boson and three colored particles were computed in [2,[11][12][13]. Having discussed how the amplitudes obtained in these papers can be converted to the MS scheme, we will now show how the hard functions can be assembled from the squared amplitudes. To obtain the result one needs to analytically continue the amplitudes to crossed channels and use parity and charge conjugation symmetry to obtain all helicity configurations from a minimal set. The amplitudes and their analytic continuation to different channels are appended in electronic form to the arXiv submissions of the papers [2,[11][12][13], since the expressions are too lengthy to be given explicitly in a paper.

Kinematics and analytic continuation
We first specify the kinematics and the analytic continuation of the amplitudes, which are common to all channels. For concreteness, we consider the specific process where the vector boson can be off the mass shell. All parton momenta are outgoing and we have The kinematic region which describes the decay of the electroweak boson to three partons is called region (1) in [32] and corresponds to the inside of the Mandelstam triangle shown in figure 1. The amplitudes for vector-boson production can be obtained from the result for V → qqg using crossing symmetry and analytic continuation. The kinematic regions relevant for the considered processes in this paper are (2), (3) and (4) in figure 1. In the crossed channels, the incoming momenta will enter with a minus sign in the definitions (3.2). The helicity amplitudes for a given process are written in terms of spinor products multiplied by coefficient functions which depend on the invariants s ij . In the following, we will denote these coefficient functions by the Greek letters α n , β n , γ n , δ n , where the subscript indicates the kinematic region, and use the letter Ω to denote a generic coefficient  (2), (3) and (4) can be obtained by analytic continuation of the result in region (1). The numbering of regions is the same as in [32], where these four regions are denoted by (1a function. It is convenient to express the functions in terms of dimensionless invariants. A suitable choice for kinematic region (1) is 1 region (1): The coefficient functions are given in [2,[11][12][13] in terms of two-dimensional harmonic polylogarithms (TDHPLs) [33,34] in the variables u 1 and v 1 defined in (3.3). In region (1) these variables fulfill the conditions 4) and the TDHPLs are analytic and real for these values of the arguments. Since some of the invariants s ij will become negative under crossing, the condition (3.4) is violated in the regions (2), (3), (4) and one has to analytically continue the amplitudes and the associated TDHPLs. A systematic algorithm for this continuation was given in [32]. In each region, one defines variables which fulfill conditions analogous to (3.4) and rewrites the amplitudes in terms of TDHPLs in these variables. They are region (2): To compute the hard functions, we then first take the IR-finite amplitude coefficients in region (i), Ω fin i (u i , v i ) for i = 2, 3, 4 included in the arXiv submissions of [2,[11][12][13]. Using the conversion formulae derived in section 2, these results can then be converted into the MS-renormalized coefficicients Ω ren i (u i , v i ). Note that also the conversion terms (2.10) and (2.12) need to be continued appropriately, but the continuation is trivial, since these terms only contain logarithms of the invariants s ij ≡ s ij + i0. To evaluate the TDHPLs numerically, we use the programs [35,36].

The qq → V g and qg → V q channels
The amplitudes for qq → V g and qg → V q can be obtained from the V → qqg amplitude using crossing symmetry and analyticity. Note that particles turn into anti-particles under the crossing operation p i → −p i . To obtain the amplitude with a W + in the final state, one therefore needs to cross the W − decay amplitude. The hadronic current S µ mediating the V → qqg decay was given in [2] in the spinor-helicity formalism [37][38][39]. Here, we will use the same spinor notation as [40] (see e.g. [41,42] for overviews of the various types of spinor notations in the literature). The spinor-helicity formalism is a four-dimensional method. The reason it can be applied here is that the Z-factor is independent of the helicities of the partons. The infrared singularities are thus an overall factor in helicity space, which can be removed after which the formalism can be applied to infrared-finite amplitudes. In [2], the fixed-helicity current S µ (q+, g+,q−) is written in the form [21] ( where R V f 1 f 2 is the right-handed coupling of the vector boson to the quarks. The electroweak couplings are given explicitly in table 1. The coefficients α 1 , β 1 , γ 1 , δ 1 are expanded as where i, j and a are the color indices of the quark, anti-quark and gluon. We note that the coefficient δ 1 (u 1 , v 1 ) is not independent, as it is linked to the other functions by current conservation. To obtain the SCET hard functions, we need to evaluate the amplitudes with JHEP02(2014)004 the renormalized coefficients Ω ren in the MS scheme. The conversion from the IR-finite amplitudes of [2,11] to the MS scheme was discussed in detail in the previous section. Since the IR singularities are independent of the spin of the particles, the coefficients are converted using the same expressions (2.10) and (2.12). We will suppress the superscript Ω ren on the amplitude coefficients in the following, but it is understood that the renormalized quantities are used to obtain the hard functions. The remaining helicity configurations follow from (3.6) by using parity conservation and charge conjugation symmetry of the strong interaction. Parity yields the relation and charge conjugation To obtain the hard functions in a given channel, we need to square the helicity amplitudes. Summing the contribution of different helicities, we then get the hadronic tensor By contracting this with the appropriate lepton tensor, one could obtain the vector-boson production cross section with an arbitrary set of cuts on the leptons arising in the vectorboson decay. Here, we will only be interested in the total vector-boson production rate and we thus only need the hard function The contribution of the helicity configuration S µ (q+, g+,q−) to the hard function in the process V → qqg, for example, is obtained from the representation (3.6) and has the form The subscript on the quantity h 1 denotes the region in which the amplitudes are evaluated. In region (n), we obtain h n (s 12 , s 23 , s 13 ) = 1 s 13 s 23 q 2 |α n (u n , v n )| 2 s 12 2s 12 q 2 + s 13 s 23 It is understood that the variables u n ≡ u n (s 12 , s 23 , s 13 ) and v n ≡ v n (s 12 , s 23 , s 13 ) in the amplitude coefficients are expressed in the invariants s ij via the relations (3.5) valid in region (n).

JHEP02(2014)004
The full hard function is obtained by summing over all helicities. Since the other helicity configurations can be obtained by applying parity and charge conjugation to S µ (q+, g+,q−), see (3.8) and (3.9), the result can be expressed using the same function h 1 . For V → qqg, one then finds (3.14) The contribution with s 13 ↔ s 23 arises from charge conjugation. For the channel qq → V g, one obtains the same expression, but the amplitudes now have to be continued to region (2) and evaluated with the variables u 2 and v 2 , Finally, for region (4), qg → V q, we need The extra minus sign in front of (3.16) compensates the minus sign which arises when crossing a fermion from the final to the initial state. The other difference to the previous two cases is that charge conjugation maps region (4) onto region (3), while (1) and (2)

The gg → V g channel
In section 2 the SCET operators for the partonic channel gg → V g were given and we now extract their Wilson coefficients from the helicity amplitudes provided in [12]. Since we only need the leading-order amplitude, which is free of infrared divergences, we can directly use the result presented in [12]. In analogy to [2], the helicity amplitudes are first given in region (1), corresponding to V → ggg and then analytically continued into the other regions. However, instead of the hadronic current, [12] provides the result after contraction with the lepton tensor which for a right-handed lepton has the form This contraction then gives a set of helicity amplitudes. Denoting the vector part of the coupling of the boson V to the quarks inside the loop by where the coefficients α i n and β i n can be found in the electronic appendix to the arXiv submission of [12]. 2 Their expansion is written in the form The color factor of the hard function, obtained from the squared amplitudes, is (d abc ) 2 = 40/3. Let us note that there is also an axial-vector contribution for the case V = Z at leading order [43,44], which is not given in [12] and will not be included in the following. All other helicity amplitudes follow from (3.20), (3.21) by permutation of the external legs, parity conjugation, and the relation A [12]. To obtain the hard function for vector-boson production, we first square the amplitudes and then integrate over the angle of the leptons in the rest frame of the vector boson to remove the lepton tensor. The integral of the lepton tensor over the direction of the leptons takes the form The same result is obtained for the left-handed current L µ L (p − 5 , p + 6 ) = L µ R (p + 6 , p − 5 ). Due to current conservation, the contraction of q µ with the hadron tensor vanishes so that, up to a prefactor, the integration over the direction is the same as contracting with g µν . The contribution of the right-handed lepton to the hard function in region (1) for the helicity configuration (+ − −), for example, is obtained from 24) and the full hard function is obtained by summing over the helicities. To perform the angle integrations, we first rewrite the product of spinor products in terms of traces, using 2 In [12] the coefficients α a n , α b n and α c n are denoted by α1, α2 and α3, and analogously for the βcoefficients. As before, we use the subscript to denote the kinematic region. In addition [12] does not include the prefactor e Q g V in the definition of the amplitudes AR.

JHEP02(2014)004
identities such as Since the result after the angle average only depends on the three gluon momenta, the γ 5 term in the trace does not contribute. We then obtain the following results for the two helicity configurations h |α a n (u n , v n )| 2 q 2 (s 23 q 2 + 2s 12 s 13 ) and h (+++) n (s 12 , s 23 , s 13 ) = (e Q g V ) 2 8q 2 s 12 |β a n (u n , v n )| 2 s 13 (s 12 + s 13 ) 2 s 23 + |β b n (u n , v n )| 2 s 23 (s 12 + s 23 ) 2 s 13 + |β c n (u n , v n )| 2 2 s 12 q 2 + s 13 s 23 − 2Re β a n (u n , v n )β b * n (u n , v n ) s 12 q 2 − s 13 s 23 + 2Re [β a n (u n , v n )β c * n (u n , v n )] s 13 (s 12 + s 13 ) As before, the subscript on h (i) n indicates that this is the result in the kinematic region (n). Opposite helicity configurations give identical contributions, and the remaining configurations are obtained by permuting the gluon momenta. Using these relations, the full hard function in region (1) is obtained as (3.28) The factor 2 accounts for the opposite helicity contributions.

JHEP02(2014)004
We are interested in the result in the kinematic region (2), which corresponds to the process g(p 1 ) + g(p 2 ) → g(p 3 ) + V (q). Proceeding exactly as discussed in the previous subsection, one first analytically continues the amplitudes α i 1 and β i 1 to region (2) and expresses them in the kinematic variables u 2 and v 2 relevant in this region. The last term in equation (3.28) requires to exchange p 1 and p 3 , and one therefore also needs the amplitude in region (4), which describes the process g(p 2 ) + g(p 3 ) → g(p 1 ) + V (q). The result for the full hard function in region (2) then reads (3.29) The analytically continued amplitudes are included in the arXiv submission of [12].
To get the lowest-order cross section in the gg → V g channel, we need to average over the spins and colors of the incoming gluons. This giveŝ In the resummed cross section, the hard function will be multiplied by a convolution of the gluon jet function J g with the soft function S gg , see (1.1). To write the factorization theorem in the form (1.1), we have introduced a hard functionĤ gg→V g , which is normalized to one at lowest order. However, the leading-order cross sectionσ

Hard functions for Higgs production
The main focus of our paper is on vector-boson production, but for completeness we now also provide the hard functions for Higgs production in association with a jet, since the construction is completely analogous to the vector-boson case. The factorization theorem (1.1) is valid also for Higgs production and it involves the same jet and soft functions as in the vector-boson case. With the hard functions given here, the resummation can be extended to N 3 LL accuracy also for Higgs production.
Since the infrared singularities of the amplitudes are independent of the spin, the conversion from Catani's subtraction scheme to MS is obtained with exactly the same formulae (2.10) and (2.12) as in the vector-boson case. The two-loop helicity amplitudes for Higgs production in the heavy top-quark limit were given in [13]. For H → g(p 1 ) + g(p 2 ) + g(p 3 ), one has The dimensionless variables u n and v n relevant in the different kinematic regions were defined in (3.3) and (3.5) and the expansion of the coefficient functions α n , β n is now written as

JHEP02(2014)004
where C t is the Wilson coefficient of the effective Lagrangian which mediates Higgs production in the large m t -limit and v is the vacuum expectation value of the Higgs field. The two-loop value of C t is listed in appendix B. When squaring the amplitudes, the color factor is (f abc The factor of two accounts for the equal opposite-helicity contributions. This result looks different than the hard function for gg → V g in (3.29) because the corresponding hard function was obtained from the helicity configuration (+ − −) instead of (+ + −). The one-loop hard function H gg→Hg was also given in [45]; we agree with this result. For the quark channel, H → q(p 1 ) +q(p 2 ) + g(p 3 ), there is a single independent amplitude, The expansion of the coefficient γ n is written in the form (3.32), but with color structure (t a ) ij instead of f abc . This yields a factor (t a ) ij (t a ) ji = C F d F (where d F = N c is the dimension of the fundamental representation) in the squared amplitude, which has the form where the factor 2 accounts for the identical, parity-opposite contributions, and the crossed terms arise from charge conjugation. The coefficients α n , β n and γ n can be obtained in electronic form from the source files of the arXiv version of [13].

JHEP02(2014)004 4 Two-loop cross section near threshold
The singular threshold cross section can be obtained by performing the fixed-order expansion of the resummed partonic cross sections, whose explicit expressions are given in [17,19]. In these papers, the resummation is achieved by solving renormalization group (RG) equations for the hard, jet and soft functions and evolving them to the factorization scale µ f , where they can then be combined with the PDFs. This allows one to evaluate each contribution at its characteristic scale. For example, to avoid large logarithms in the hard function, the starting scale µ h of the RG evolution is chosen to be µ h ≈ p T . For the jet and soft functions lower starting scales µ j and µ s are appropriate, as discussed in detail in [17,19]. The simplest way of going back to the fixed-order expressions is to switch off the resummation by taking the limit where the scales µ h , µ j , µ s and µ f all become equal. In Laplace space, where the cross section factors into a product of the hard, jet and soft functions, taking the limit is completely trivial. All the RG-evolution factors switch off and the fixed-order result is simply the product of the Laplace-transformed functions. In momentum space, the limit is a bit more delicate since the cross section becomes distribution valued in this limit. Starting with (20) in [19] and taking the limit in which all scales coincide whenever it is trivial, we are left with d 2σ sing abc dy dp 2 for the channel a + b → V + j c . To factor out the tree-level cross sectionσ (0) ab , we have normalized the hard functionsĤ ab to one at the lowest order and have indicated this by putting a hat on the normalized functions. The jet functionj c and the soft functions ab appearing here are the Laplace-transformed functions. Their explicit form can be found in appendix B. At the n-th order in perturbation theory, the functions are polynomials of order 2n in logarithms of the Laplace variables. In the above representation, these logarithms are replaced by derivatives acting on an m X -dependent kernel. This type of solutions for the RG equations of the jet and soft functions was introduced in [46]. To take the limit η → 0, we first need to expand this kernel in a series of distributions The ⋆-distributions appearing on the right-hand side are generalizations of the usual plusdistributions to dimensionful variables. Their explicit form can be obtained by rewriting the integration over the invariant mass m X in the form The expansion of the cross section is now straightforward. To present the result, we write the perturbative expansion of the normalized hard function in the form and also introduce expansion coefficients p (n) i which capture the contribution of the product of the jet and soft function at the n-th order in perturbation theory. The coefficients p (n) 0 multiply δ(m 2 X ) and the higher coefficients 0 < i ≤ 2n the ⋆-distributions arising in the expansion (4.2). The coefficients h (n) and p (n) i depend on the partonic channel, but in the following we suppress the channel indices a, b and c for better readability. To two-loop order, the cross section then has the structure d 2σsing dy dp 2 The explicit form of the one-loop coefficients in the above formula is The lengthy two-loop coefficients p (2) i are listed in appendix C. The Casimir operators relevant for the different channels are

JHEP02(2014)004
The nonlogarithmic one-loop coefficients of the gluon [17] and quark [47,48] jet functions read while the coefficients for the soft function read [17] c Sqq 1 The two-loop coefficient for the quark jet function has been calculated in [8], for the gluon jet function in [9], and the two-loop coefficients for the soft functions have been calculated in [10]; they are listed in appendix B. To show the relative size of the corrections, we have normalized the hard functions for the qq and qg channels to one at the lowest order, as indicated by the hat on the normalized functions. With α s (µ) ≈ 0.09, the corrections are moderate, of the order of a few per cent. The contribution proportional to N v,a V arises from diagrams where the vector boson couples to an internal quark loop instead of the external quarks. By charge conservation, such contributions do not arise for W -bosons, so that N v,a W ± = 0. For photons, we have

Numerical studies
where the sum runs over the quark flavors in the loop and the denominator arises because we have factored out the charge in the definition of the hard function. For Z-bosons, there are contributions from both the vector and the axial part of the coupling. The axialvector part at one-loop order can be found in [49]. In references [2,11] the two-loop vector JHEP02(2014)004 part was computed, but the two-loop axial part is at present still unknown. The relevant couplings are Note that the contributions proportional to N v,a V are numerically very small. The normalization of the constant N v Z differs from N F,Z in [2,11] because we consider the squared amplitude instead of the amplitude itself.

(5.4)
These numbers do not include the small perturbative corrections to the Wilson coefficient [C t (m t , µ)] 2 in (3.33). We observe that the higher-order terms are dramatically larger than in the vector-boson case, in line with the findings of [7]. Larger corrections are expected since the higher-order contributions to gluonic quantities are enhanced by factors of C A /C F and also because the leading cross section is O(α 3 s ) for the Higgs case instead of O(α s ) as in vector-boson production. However, for a meaningful assessment of the size of the higherorder corrections, one will need to evaluate the full cross section. Also, to make reliable predictions, one should check how large the corrections to the heavy top-quark limit are. The above values correspond to p T ≈ 0.5 TeV for which the effective theory treatment is no longer appropriate. We will study the Higgs case in more detail in the future and will restrict ourselves to vector-boson production in the following.
In our previous work [18,19], we have used a Mathematica code to compute the cross sections. With the large size of the expressions for the two-loop hard functions, this code becomes prohibitively slow and we have now developed a C ++ code PeTeR [25] to compute the cross section, which will be made public in the future. The code computes the resummed cross section near the partonic threshold as well as its fixed-order expansion. In addition, it also computes the full NLO fixed-order cross section. In a future paper, we will present a detailed phenomenological study of vector-boson production, including the twoloop corrections as well as electroweak Sudakov effects, which were recently treated using the same threshold-resummation framework [50]. For the moment, we focus on the size of the two-loop QCD corrections and check how much they change the cross section. To do so, we use the same input parameters as in our previous paper [19], namely the NNLO MSTW 2008 PDF set and its associated α s (M Z ) = 0.1171 [51]    We treat all partons as massless except for the top quark, which is integrated out from the theory. The hard function for Z-boson production receives tiny contributions from the axial-vector coupling, see (5.1). At one-loop order they are due to triangle diagrams. A similar contribution is present for the gg channel [43,44], which is not included so far but might be of a similar order of magnitude as the NLO triangle contribution. For simplicity, and because the two-loop axial corrections are not known, we set N v V = 0. Numerically the two-loop N v V terms are negligibly small. A list of values for the integrated cross section σ(p T > 200 GeV) is shown in table 2 for different LHC center-of-mass energies. The table presents three different approximations i) the fixed-order threshold cross section ii) the resummed results, and iii) the results obtained after matching to the known NLO fixed-order result. The entries LO, NLO sing. , NNLO sing. show the perturbative expansion of the threshold cross section, which consists of the singular distributions defined in (4.2). Since the LO partonic cross section is proportional to δ(m 2 X ) it is purely singular. Beyond leading order, the cross section also has regular pieces not associated with soft and collinear radiation. As the table shows, the regular pieces obtained from the difference NLO−NLO sing are of moderate size. For example, for Z-production at √ s = 8 TeV, the singular pieces amount to about 70% of the NLO correction. The fact that the singular pieces amount to the bulk of the cross section is true in many other cases as well, and we therefore expect that the singular pieces will provide a good approximation to the full NNLO correction. The column NNLO sing. +NLO shows the result obtained if both the full NLO result and the singular pieces at NNLO are included. For the factorization and renormalization scales, we use This value is close to p T and was adopted as the default scale µ h for the hard function after a numerical study in [19]. The scale uncertainty is obtained by varying the scale µ by a factor two around the default value. Figure 2 shows the resulting uncertainty bands for Z and for W + or W − production at NLO and NNLO sing +NLO. The results are normalized to the NLO result at the default scale choice. We find that including the two-loop singular terms corresponds to a shift of about +5% of the cross section and decreases the scale uncertainty by a factor of two, compared to NLO. The factorization formula (1.1) can be used to resum the singular pieces to all orders using RG evolution in SCET. For W and Z production, this was done in [18,19]. To perform the resummation, we evolve the hard, jet and soft functions from their characteristic scales to the factorization scale. We adopt here the same default scales as in the papers [18,19]. Resummation to N n+1 LL accuracy requires the hard, jet and soft functions at N n LO. Comparing the resummed N n+1 LL and the threshold fixed-order results N n LO sing. , we find that they are numerically very similar. Resummation is thus not a large effect, since the characteristic scales for the jet and soft functions are not much below the hard scale. Their numerical values depend on the fall-off of the PDFs towards larger x, which enhances the threshold region. After a numerical study, following [52], the values and µ s = µ 2 j /µ h were adopted in [19]. Since the numerical values of the jet and soft scales are not much lower than the hard scale, the logarithms which are resummed are of moderate size.
The highest order result obtained in [18,19] was denoted by N 3 LL p where the label "p" (for partial) indicated that the two-loop hard function was missing. With this ingredient in place, our results now have full N 3 LL accuracy. 4 However, since all the logarithmic pieces of the hard function follow from RG-invariance, they were already included in N 3 LL p of [18,19]. The logarithms were introduced such that their contribution vanishes at µ =  Figure 3. Left: relative change of the NNLO sing cross section compared with an approximation used in [18,19], in which the two-loop constant was chosen to vanish at µ = p T . Right: relative contribution of triangle diagrams at NLO (red line) and of the gluon-gluon channel at NNLO (times 5, blue dashed line).
i.e. the two-loop constant was defined as the value of the hard function at µ = p T . To see how much the results in [18,19] change due to presence of the two-loop constant, we plot its value compared to the cross section without the two-loop constant. The size of the effect depends on the transverse momentum. As shown in the left panel of figure 3, it reduces the cross section by about two per cent at p T ≈ 100 GeV and enhances it by a similar amount at large p T ≈ 2 TeV. The change is within the scale uncertainties of the N 3 LL p results and smaller than the total two-loop effect, which is of order +5%, see figure 2.
There are two more effects, which we briefly address. The first is the contribution of the gluon-gluon channel in Z production, which first arises at NNLO and is shown by the blue dashed line in the right panel of figure 3. The channel gives a very small positive contribution to the cross section. It peaks around p T = 50 GeV and is smaller than 1‰. The second effect, shown by the red line in the right panel of the figure is the triangle contribution at NLO which arises due to the axial coupling of the Z boson. This contribution is mostly negative and smaller than 3‰. As we stressed earlier, there are also axial contributions at NNLO, in particular also in the gluon-gluon channel, which were obtained in [43,44] but are not included here. According to these papers, the axial corrections to the gluon-gluon channel are bigger than the vector contributions and could be comparable in size to the (very small) axial contributions at NLO.

Conclusion
Radiative corrections to hard-scattering processes simplify considerably near the partonic threshold. In this work, we have used these simplifications to obtain the NNLO corrections to transverse-momentum spectra of photons, W , Z and Higgs bosons. Our results are valid at large transverse momentum p T of the electroweak boson, where the invariant mass of the recoiling jet is small compared to p T . As the threshold terms often capture the bulk of the radiative corrections, we expect that our results are a good approximation to the exact NNLO results. In addition, the threshold terms can serve as a check on the full NNLO results, once they become available.

JHEP02(2014)004
The starting point of our analysis is the threshold factorization of the cross section into hard, jet and soft functions. Building on earlier work, in which we computed the two-loop collinear and soft functions, we computed the last missing NNLO ingredient, the hard functions, in this work. To this end, we converted known results for on-shell V +jet amplitudes into MS-subtracted hard functions defined in Soft-Collinear Effective Theory. The conversion procedure presented in sections 2 and 3 is completely general, and applies similarly to other processes. Our calculation also provides the last missing ingredient to resum the threshold terms to N 3 LL accuracy.
We have implemented the NNLO threshold corrections and the N 3 LL resummed results into a C ++ code PeTeR [25], which will be made public in the future. For W and Z production, we find that the NNLO threshold corrections are moderate. They enhance the cross section by about 5%, and they reduce the scale uncertainty by about a factor of two. In addition, we have also given resummed results at N 3 LL accuracy, matched to NLO fixed-order results. Numerically, we find that the resummation effects, i.e. terms beyond NNLO, are not very important. Our final results for the integrated cross sections with p T > 200 GeV are given in the last two lines of table 2. For Higgs production, the corrections to the hard functions are much larger than in the vector-boson case and resummation will likely be more important. We will present numerical results for the Higgs cross section in the future.
For an accurate description of LHC data for vector-boson production at high-p T , one also needs to implement electroweak corrections which are large and negative. These Sudakov-type corrections have recently been studied using the same threshold-resummation formalism. In a next step, we will perform a detailed phenomenological analysis of vectorboson production, including both electroweak and QCD corrections. The hard functions determined in the present paper are relevant not only for hadronically inclusive boson production, but are also needed for resummations of more exclusive one-jet observables such as jet-mass spectra or jet-veto cross sections.

JHEP02(2014)004
We define the QCD β-function and its expansion as so that the lowest two coefficients have the explicit form In the following, we will expand all anomalous dimensions in units of α s /4π, and we denote the expansion coefficients in the form Up to two-loop order, the cusp anomalous dimension γ cusp is given by and the collinear anomalous dimensions γ q and γ g are The renormalization factor Z is obtained by solving its RG equation, which is driven by the anomalous dimension matrix Γ in (2.5). The two-loop expression has the form ln Z = α s 4π Expanding the inverse Z-factor in units of α s /2π, one obtains

JHEP02(2014)004
The one-loop subtraction operator appearing in Catani's formula for the IR divergences is Apart from the pole terms, we will need the explicit expressions for the first two coefficients and The two-loop subtraction operator is defined as where the last term has not been specified in [31], but is was stated that it only contains single poles. Using the expression for the Z-factor (A.7) one can derive this term. Explicitly, we find where the two sums run over all unordered triplets of distinct parton indices. The terms in the first two lines are equal to 1 8 Γ 0 , C 0 . This commutator can be simplified by noting that the contributions which involve four different partons vanish because the color generators associated with different partons commute. This expression for H (2) R.S. (ǫ) was derived in [26,27], but the term in the second line, which involves the collinear anomalous dimensions γ i 0 was missed. This extra contribution was discussed in appendix D of [53], where it was shown that it can only contribute for amplitudes with more than four external particles.
The two-loop conversion relation in (2.12) involves a commutator of the one-loop anomalous dimension Γ 0 with the O(ǫ) term in the expansion of I (1) (ǫ). This commu-

JHEP02(2014)004
tator can be simplified to B Hard, jet and soft functions at two-loop order

B.1 Hard function
The helicity amplitudes in [2,[11][12][13] were only given for the choice µ 2 = q 2 . The full µ dependence can be reconstructed by solving the associated RG equation [26], driven by the anomalous dimension Γ in (2.5). One finds [17] H û,t, The logarithms for the different channels are The anomalous dimensions can be extracted from the general result [27], explicitly, The results for the other crossed channels like gq can be obtained as usual by replacinĝ t ↔û.

JHEP02(2014)004
For Higgs production, the last term of γ Hqq (α s ) and γ Hqg (α s ) must be changed to − 3β(αs) 2αs because the leading-order cross sections start at O(α 3 s ). Also, the above result is relevant for the full hard function, which includes the factor [C t (m 2 t , µ 2 )] 2 , from the Wilson coefficient of the operator (3.33) which mediates Higgs production in the large-m t limit. The NNLO value of this coefficient is [54,55] C t (m 2 t , µ 2 ) = 1 + If the scale-dependent factor [C t (m 2 t , µ 2 )] 2 is divided out, the anomalous dimension of the hard function changes by ∆γ H ab (α s ) = −2γ t (α s ), where The anomalous dimension γ t (α s ) is related to the QCD β-function [56,57] since the operator is proportional to the Yang-Mills Lagrangian.

B.2 Jet function
The expression for the Laplace-transformed jet functionj c (L, µ), with c = q or c = g, reads The Casimir operators C S ab for the different channels are (B.13)

JHEP02(2014)004 C Coefficients of the two-loop threshold cross section
Here we list the expansion coefficients that appear in the two-loop threshold cross section (4.6). The one-loop coefficients p