The Higgs Transverse Momentum Distribution at NNLL and its Theoretical Errors

In this letter, we present the NNLL-NNLO transverse momentum Higgs distribution arising from gluon fusion. In the regime $p_\perp\ll m_H$ we include the resummation of the large logs at next to next-to leading order and then match on to the $\alpha_s^2$ fixed order result near $p_\perp \sim m_h$. By utilizing the rapidity renormalization group (RRG) we are able to smoothly match between the resummed, small $p_\perp$ regime and the fixed order regime. We give a detailed discussion of the scale dependence of the result including an analysis of the rapidity scale dependence. Our central value differs from previous results, in the transition region as well as the tail, by an amount which is outside the error band. This difference is due to the fact that the RRG profile allows us to smoothly turn off the resummation.


Introduction
We are transitioning into an era of precisions Higgs physics. Presently there is sufficient data to study various decay modes of the Higgs [1][2][3], and soon there will be sufficient data to study the the Higgs differential cross section, which can be used to better understand the underlying production mechanism and to search for new physics [4][5][6][7][8]. While the new physics would be most prominent at large values of p ⊥ , the predominance of the events will be in the lower range. Furthermore, to isolate the Higgs' decays from backgrounds, events are binned according to their highest-transverse momentum jet. The 0-jet bin, corresponding to no central jets above a certain transverse-momentum threshold, plays a critical role in the current Higgs analysis. 1 Thus there is increased motivation for making precise predictions for the distribution when p ⊥ m h since back to back central jet production is power suppressed [27].
The small p ⊥ region of parameter space is polluted by large logarithms which must be resummed in order to retain systematic control of the theoretical errors. Resummations have been previously been discussed within an SCET [9][10][11][12] framework [29,30,[32][33][34][35][37][38][39] as well as in the CSS resummation formalism [40][41][42][43][44][45]. While formally most of the these JHEP12(2015)097 results agree at a given order in the resummation procedure, the central values as well as the predicted errors between different resummation techniques will differ. 2 Another crucial difference between various theoretical predictions is how the transition between the fixed order result at large p ⊥ and the resummed result is handled. In this paper we present a result for the resummed cross section at NNLL + NNLO within the confines of SCET and the rapidity renormalization group (RRG) [26,27]. Our motivation for this analysis is two-fold. By working within the RRG formalism we are able to consistently turn off the resummation in the tail region. It has been shown in the context of B → X s + γ [60], thrust [61] as well as the Higgs jet veto calculation [56], that if one does not turn off the resummation then one can over estimate the cross section, in the region where fixed order perturbation theory should suffice, by an amount which goes beyond the canonical error band in the fixed order result. This overshoot happens despite the fact that the resummed terms are formally sub-leading in the expansion. The reason for this overshoot has been shown [56,60,61] to be due to the fact that there are cancellations between the singular and non-singular terms in the tail region and that this cancellation will occur only if the proper scale is chosen in the logarithms. This will be born out in our analysis as well. By using the RRG we have a natural way to interpolate between the resummed and fixed order regions. This allows us to present the full spectrum such that in the large p ⊥ region our results match onto the NNLO result, as previously discussed within the context of the jet veto cross section [17,18]. In this sense our paper fills a gap in the literature. As we will discuss below, the RRG plays a important role in the theory error determination. We also discuss in detail how our error estimates and full spectrum differ from those of previous work.

Systematics
In this section, for the sake of completeness, we review the well known systematics of the calculation. We work within the confines of the large top mass approximation. This approximation is known to work extremely well -much better then one would naively expected -especially away from the tail of the spectrum. We will be not be considering p ⊥ large enough for these corrections to be relevant to the error budget. A full discussion of such errors as well as others will be discussed at the end of the paper.
At the scale m t , full QCD is matched onto the large top mass effective theory at order α 2 s , after which we match onto SCET at the scale m h . In SCET we will work to leading order in a systematic expansion in p ⊥ /m h , and the errors due to non-perturbative corrections, which are suppressed by Λ QCD /p ⊥ , are included in the final error analysis.
As is well known, logs in the perturbative expansion -in impact parameter spaceexponentiate and allow us to organize the series as follows

JHEP12(2015)097
The standard terminology is such that keeping only F corresponds to "LL" leading (double) log, while if G is retained we would have NLL etc. As p ⊥ gets larger the EFT breaks down and the fixed order calculation becomes the relevant quantity. We will utilize the order α 2 s (NNLO) result [19,20] 3 in this large p ⊥ region. We do so despite the fact in the small p ⊥ region we will be only keeping terms of order α s . This is not inconsistent as the error in the two disparate regions are distinct and uncorrelated. In the resummed region we will be working at NNLL in the perturbative expansion and at leading order in the power expansion p ⊥ /m h . By utilizing the RRG scale we will smoothly turn off the resummation and match onto the full NNLO fixed order calculation.

Factorization in SCET and anomalous dimension
A factorization theorem for Higgs production at small p ⊥ within SCET and the RRG formalism was developed in [27]. The starting point in the large top mass effective theory is the gluon fusion operator The matching coefficient for this operator (C t ) is known to two loops [24]. At the scale m h , 4 we match onto SCETII. This theory is the version of SCET in which the collinear modes and the soft modes have the same invariant mass, which distinguishes it from SCETI where there exists a hierarchy in masses between these two modes. Due to this equality in virtualities a new set of divergences -rapidity divergences -arise, which are not regulated by dimensional regularization. The rapidity regulator introduces a new scale which acts as a boundary between the collinear and soft modes as discussed in [26,27]. The reader may consult ( [27]) for the details of the formalism. At leading order in λ = p ⊥ /m h the differential cross section for higgs production at low transverse momentum may be written as S is the soft function defined as defined in terms of the light-like Wilson lines As is common with fixed order calculations, the authors of [20] use the term NLO, since the leading order result is a delta function in p ⊥ . 4 We ignore the running between the top mass scale and m h as these logs will be sub-leading in our power counting.

JHEP12(2015)097
f αβ is the transverse momentum distribution function (TMPDF) which is matched onto the PDFs at the scale p ⊥ Λ.
W n is a collinear Wilson line in the fundamental representation defined in x-space by n/n = (1, 0, 0, ±1) are the null vectors which correspond to the directions of the incoming protons. In impact parameter space, we can write the functions S and f αβ in terms of their inverse Fourier transform.
The cross section is given by H(m h ) is the hard coefficient which depends on the scale µ = m h and and we define H(m h ) = 8m 2 h |C s | 2 , where the coefficient C s is known to two loops [25]. The soft function as well as the TMPDFs depend upon both µ, the usual scale introduced within the context of dimensional regularization which factorizes the hard modes from the soft and collinear modes, as well as ν, the rapidity factorization scale which distinguishes between the soft and collinear modes.

Renormalization and resummation
Each of the pieces of the factorization theorem H, S and f , have natural scales with which they are associated. H is independent of ν and its natural µ scale is m h . The soft and collinear functions have natural scales (µ = p ⊥ , ν = p ⊥ ) and (µ = p ⊥ , ν = m h ) respectively. When these objects are evaluated at these scales, they will be devoid of large logarithms. However, given the natural distribution of scales we can see that it is not possible to choose a µ and ν such that all of the individual pieces sit at their respective natural scales. Thus it is expedient to choose a (µ, ν) value such that the maximum number of pieces sit at their natural scale, and then use the RG and RRG to sum logs for the remaining pieces.  The individual functions satisfy the following RGE's.
where κ i = (µ, ν) and F ∈ {S, f αβ }. The anomalous dimensions satisfy the relations All the renormalization group equations are diagonal in impact parameter space and hence are straightforward to solve. The scales for RG and RRG operations commute, and as a consequence of this, the following relations are generated where Γ cusp is the cusp anomalous dimensions of two light-like Wilson lines [49]. Figure 1 shows the path in (µ, ν) space we have chosen to resum the logs. The hard part is run from the scales m h down to the scale 1/b 0 ∼ p ⊥ , while the jet function is run in ν space up to the scale m h , where the scale b 0 is defined as be γ E /2. In performing the hard resummation, π 2 are resummed by analytically continuing the matching scale from space-like to time-like kinematics relevant for Higgs production [50]. It has been shown that the inclusion of the π 2 in the resummation does improve the perturbative convergence of the hard matching expansion [31,53], however this method only sums a subset of π 2 and one can not claim that this method leads to a systematic reduction in theory errors. 5

JHEP12(2015)097
At NNLL, the systematics require that we keep the cusp anomalous dimensions at three loops, the non-cusp at two loops and matching at both the high and low scales at order α s . All of the necessary anomalous dimensions have been calculated previously in the literature. The hard anomalous dimensions can be extracted from [55]. We utilize the one loop matching for the jet and soft function from ( [27]). To run the jet function in ν we need only the two loop cusp piece since log in jet rapidity the anomalous dimension is not large.
Given the consistency conditions (eq. (4.4)), the anomalous dimension for the jet function in ν is of the form γ f ν = Γ cusp log(µb 0 ) + γ f ν,non−cusp . For the resummation path that we have chosen, µ ∼ 1/b 0 , so that the γ f ν ∼ Γ cusp +γ f ν,non−cusp , which implies that for ν running at the order NNLL, we only require the two loop cusp and two loop non-cusp pieces to the logarithmic order we work. The running of the hard function in µ still incorporates the full 3 loop cusp anomalous dimension, so that the the running in both µ and ν is consistent at NNLL. We also the need the two loop "non-cusp" piece in ν which can be obtained from the two loop soft function or ref. [32]. The results for the hard anomalous dimension are tabulated in [56].
The matching of the TMPDF onto the PDF may be written as where the sum is on species of partons, and the parton PDF's are defined as We had adopted the mostly minus metric such that conventions that p α ⊥ p β ⊥ g ⊥αβ = − p 2 ⊥ and make use of the 't Hooft-Veltmann scheme, so the external transverse momenta remains in 2 dimensions, as do the external polarizations on the operator (the free Lorentz indices). The scheme choice is advantageous, as it allows one to renormalize the operator directly in p ⊥ space. However, the higher order anomalous dimensions needed for NNLL resummation have not been evaluated in this scheme, so for the purpose of resummation in this paper we follow the conventional dimensional renormalization(CDR). At tree level the TMDPDF matches on to the PDF as follows: The one loop renormalized TMPDF, necessary to extract the order α s matching coefficients JHEP12(2015)097 ⊥ig/g ), is given by where we have written the expression in terms of plus distribution L n = 1 2πµ 2 , the properties of this distribution are collected in the appendix of [27].
The resummation is carried out in impact parameter space and when Fourier transforming back to p ⊥ space, the integrand diverges along the path of integration due to the Landau pole at 1/b ∼ Λ QCD . However looking at the form of the cross section, the integrand over b contains the Bessel function of first kind with argument bp ⊥ and an exponent which contains powers of ∼ log(m h b). The combination of these two factors quickly damps out the integrand outside the region 1/p ⊥ > b > 1/m h . Thus to avoid hitting the pole, we can simply put a sharp cut off for b at a value sufficiently higher than the largest value of 1/p ⊥ ∼ 0.5 GeV −1 but smaller than 1/Λ QCD . A convenient value is found to be b ∼ 2GeV −1 . It has been numerically ascertained that the variation of this cut-off between the values of 1.5 − 3 GeV −1 produces an error only in the fourth significant digit of the integral.

The fixed order cross section and power corrections
When we match the full theory onto SCET we drop terms which are beyond leading order in λ = p ⊥ /m h . However, a comparison of the full theory fixed order cross section( [20]) with the effective theory calculation reveals that in order to achieve ten percent accuracy, we need to include these higher order terms in λ beyond p ⊥ > 30 GeV. In principle we could achieve this by keeping higher order terms in the factorization theorem. However, this would only be necessary if we wished to resum the logs associated with these power corrections. Given that these power corrections are only relevant in regions of larger p ⊥ , these logs are not numerically large and thus represent a small corrections which need not be resummed. At large p ⊥ we should turn off the resummation and use the fixed order NNLO result.
where dσ res dp 2 ⊥ is the resummed cross-section at NNLL computed at leading order in SCET.
This includes all the logarithmically enhanced terms at NNLL. dσ ( fin) dp 2 ⊥ contains all the rest of the terms not captured by the effective theory at NNLO. Mathematically, we can define dσ (fin) dp 2 where dσ (res) is obtained by expanding out the resummed exponent and truncating the terms to the desired order in α s , while dσ is the complete fixed order cross section at the same order in α s . This is made explicit in the matching done for the particular case of √ s = 13 Tev as shown in figure 2. This clearly shows the relative importance of the finite terms as compared to the logarithmic terms, giving us an idea of the scale at which the resummation needs to be turned off. Thus we subtract out the logs from the NNLO correction so that we do not double count, and turn off the resummation using scale profiles which we now consider.

Profiles
The SCET formalism we have used retains the leading order operators in the power counting parameter λ. This restricts the validity of the resummed cross-section to the regime where p ⊥ m h . At larger values of of p ⊥ /m h the resummation becomes superfluous. Moreover, the power corrections from subleading operators become increasingly relevant and to maintain accuracy it is vital to include these power corrections, as discussed above. So it is desirable to smoothly switch over from the resummed result to the full theory fixed order cross section when the singular and non-singular terms of the cross-section are comparable. Notice that formally, it is not necessary to turn off the resummation. As long as one makes the proper subtraction to avoid double counting, then the resummed terms that are kept should be sub-leading in the logarithmic power counting, since we are in the transition region. However, as discussed in the introduction, the fixed order result involves the cancellation between singular and non-singular terms. If the scale of the log

JHEP12(2015)097
is not chosen appropriately (of order p ⊥ ) then this cancellation is not effective and the cross section is over estimated beyond the canonical scale variation errors. 6 This will be demonstrated below.
The shutting down of the resummation is achieved by varying the low matching scale from its value 1/b 0 in the low p ⊥ region to m h in the high p ⊥ region. This has the effect of turning off the resummed exponent. To do this, we introduce profiles in µ, ν following the work done in [17,18]. We use three typical profiles for varying the renormalization scale from 1/b 0 to m h . Each profile is chosen as a linear combination of hyperbolic tangent functions which smoothly transitions from the resummation region to the fixed order one. The general equation for each profile R (where R can be µ, ν) can be written as where s determines the rate of transition, L = 2e −γ E /b = 1/b 0 is the initial value and H = m h is the final value. t is the value about which the transition is centered. Figure 3 shows three profiles with s = 2 and t = (35, 45, 55) GeV . One may worry that since the profile depends on the transverse momentum, the fourier transform is not truly a fourier transform. Since the profiles asymptote to constant functions, in both the peak region and the tail region, one effectively has a "pure" fourier transform in regions where there is a well-defined logarithmic power counting or the result is simply the fixed order cross-section. That is, at low values of p ⊥ all the logs log(m h b 0 ) are in the exponent, at large values they sit in the fixed order result. In the transition region, the profiles do exhibit dependence on p ⊥ , causing departures from a pure fourier transform. But this dependence formally cancels with the fixed-order terms/effective theory matrix elements to the logarithmic order we work, and the subleading residual effect is probed by our scale variation in the transition region. Thus departures from "purity" are included in the error budget. This cancellation is rather strong, as can be seen from the variations of the profiles, and this is also consequence of the fact that in the transition region, the large logs are ceasing to be large as one approaches the tail. Indeed the issue of a impure fourier transform also appear if we had scale set in p ⊥ space, due to the simple fact that the same expressions obtain whether one fourier transforms from b-space to p ⊥ space, then scale sets in p ⊥ , or scale sets with p ⊥ and then moves from b-space to p ⊥ space.

Error analysis
In our formalism, we have performed an expansion of the cross-section in m h /m t , p ⊥ /m h , Λ QCD /p ⊥ and the strong coupling α s . Therefore, the dominant error in our cross-section is due to the first sub-leading term that we drop in each of these parameters. When expanding in m h /m t , we retain only the leading order operator. A comparison of the cross section in this limit with the exact m t dependent result [20][21][22][23], reveals that this approximation works extremely well for p ⊥ < m t . For the range of p ⊥ discussed in this paper, the error due to these corrections is less than 1 %. For the parameter p ⊥ /m h , again we retain only JHEP12(2015)097 the leading order operator. This approximation works well for p ⊥ ≤ 30 GeV as the leading correction scales as p 2 ⊥ /m 2 h . As discussed earlier, for larger values of p ⊥ , we include the full NNLO expression which includes all orders in p ⊥ /m h .
For the third parameter, Λ QCD /p ⊥ one must consider the non-perturbative contributions. The reason is that we are working in SCETII where the collinear and soft function both have the same invariant mass scale. As such, there will be non-perturbative contributions coming from both of these sectors. As p ⊥ approaches Λ QCD , we need to systematically include higher dimensional operators in both the soft and the beam sectors. Thus there are multiple non-perturbative functions which contribute in this regime. In principle these matrix elements can be fit to the data. Work in this direction has been performed for vector boson production where ansätzes [59] for these matrix elements have been fit to the data [57]. However, since we are interested in a gluon initiated process we cannot hope to use any extraction from the quark initiated vector boson production process. Therefore we make no attempt to model these corrections and simply include a rough estimate of the errors due to these terms by assuming an error which scales as Λ 2 QCD /p 2 ⊥ . We use a conservative value of Λ QCD of 1 GeV for this term. The effective theory that we have formulated works at leading order in the λ ∼ Λ QCD /p ⊥ and hence is valid only in the regions where this parameter is small. For large values of p ⊥ , higher order corrections for this parameter are negligible but only become important as we approach p ⊥ ∼ Λ QCD . Many authors 7 do include a Gaussian non-perturbative model (with a single function of b to be fitted to data) to account for these corrections, which allow them to extrapolate the cross section to non-perturbative values of p ⊥ . This corresponds to augmenting with a shape function the soft function of our factorization theorem. But as stated earlier, in the region of interest, we have multiple functions beyond the pdf's which need to be taken in to account for a systematic analysis. As far as our results are concerned, we do not claim to have an accurate result for p ⊥ ≤ 3GeV . Furthermore, since the leading order correction JHEP12(2015)097 is of the form λ 2 , this amounts to a correction 2 − 10% in the range of 5-3 GeV, consistent with the perturbative accuracy of the cross section in the entire range of p ⊥ .
Finally we consider the errors due to higher order perturbative corrections. The accepted methodology for estimation of these errors is to vary the scale µ by factors of two. The idea is that such variations estimate the size of the constant terms which are not captured by the renormalization group, as well as the sensitivity to the perturbative truncation of the anomalous dimensions. Of course this is just a rough guess which falls into the rubric of "you do the best you can," though it is important to try to gauge all sources of large logarithms in these variations.
When working with the RRG we generalize this technique since we now have two distinct types of factorizations, virtuality and rapidity, and within the RRG there are therefore two distinct exponentiations. In both cases there is a choice as to what subleading terms should go into the exponent. Varying the scales µ, ν corresponds to varying the size of these sub-leading pieces. Given that virtuality and rapidity factorization are independent mechanisms we should vary both of these scales in order to estimate the errors in the choice of exponents. This variation is accomplished by adjusting the parameter L in each profile which sets the scale for the low scale matching of the jet and soft function that dominates the perturbative errors. In keeping with the canonical recipe adopted by the community we vary L by a factor of two about its central value (1/b 0 ). The corresponding modification in the profiles for either µ or ν is shown in the figure 4. When we vary these scales we must keep in mind the restriction that the argument of each large log is varied at most by a factor of two. This gives us a constraint 0.5 ≤ µ/ν ≤ 2 since there are logs of the form log(µ/ν) in the exponent. In all, this provides us with six different possibilities for each of the three profiles ( figure 3). We do a similar analysis for the end point region of each profile by varying the scale H ( figure 5). Since this region of transverse momentum is dominated by the fixed order cross section, the variation in H amounts to a fixed order JHEP12(2015)097  scale variation. We choose to keep the largest error band generated due to these variations. The plots which combine all the effects above, are shown in figure 6 and figure 7 for S = 8 and 14 TeV. We use the MSTW 2008 pdf at NNLO [63] and the value α s (m z ) = 0.1184.
We also consider how the cross section is typically affected by using other pdf sets. In particular we evaluate the cross section at 13 TeV with the CTEQ5 pdf.
The deviation form the MSTW2008 result is of the order of a few percent as shown in figure 9. We include these errors as well in the final plot for the cross section at 13 TeV figure 8.

Comparison to previous results
Let us now compare our results with the previous results in the literature which include the NNLL resummation. We will compare only with those papers which show plots which include all of the corrections considered here, i.e. [33,64]. We agree formally with the results of these papers. Of course formal agreement is not necessarily the relevant issue. Formal agreement means that the coefficients of the logs in the resummed result (as well as the non-singular pieces) agree. However, it does not mean that the argument of the logs are the same. Indeed for sufficiently different arguments the results can differ by an amount which is well outside any "reasonable" theoretical error estimation.  In the work [33] the authors calculate the spectrum up to 60 GeV and match onto the NLO result. Terms of order α 2 s are not included which is consistent as long as one presumes that αs π ln(60 2 /m 2 h ) ∼ 1. That is, as long as the logs still dominate all the way up to p ⊥ = 60 GeV, then keeping the terms of order α 2 s would not improve the accuracy of the calculation. The results in [33] do not track the rapidity scale. Meaning that an implicit fixed matching scale has been chosen. In contrast to [33], the NNLL error band in our analysis lies within the error band of NLL. The reason is mainly two fold. They have JHEP12(2015)097 no rapidity scale to vary and they have no errors due to non-perturbative corrections. This elucidates the importance of including the scale ν in error analysis for a better estimation of accuracy at each order in resummation. Further, the authors claim the existence of a non-perturbatively generated hard scale on the order of 8 GeV, which they state cuts off the non-perturbative corrections and therefore there results are valid down to vanishing p ⊥ . The underlying justification is based on the paper [62] where it was argued that since in QED the probability to emit a photon with p ⊥ ∼ 0 vanishes faster then any value of Q 2 , the only contribution to the rate must come from two photons whose total p ⊥ approximately vanishes but whose individual p ⊥ is large. The authors of [62] then conjecture that the same reasoning should hold in the non-Abelian case. Figure 10 shows a comparison with their central value at NNLL+NLO. We have not included their error bands in this case since they match only to NLO. For larger values of p ⊥ , the central value curve is below our estimates since the matching is done with fixed order NLO.

JHEP12(2015)097
Grazzini et al. [64] use the CSS formalism and plot the result out to 120 GeV, matching onto the fixed order result at NNLO. They do include an additional scale they call the resummation scale Q, in addition to the factorization scale of the PDFs and the renormalization scale of α. The variation of Q variation mixes what we would call µ and ν variation. This is discussed in detail in the appendix. To turn off the resummation in the exponent they make the replacement this has the effect of killing the resummation at asymptotically small b and insuring that the total inclusive cross section in reproduced. However, in the range m h /p ⊥ ∼ 1 the

JHEP12(2015)097
resummation is not shut off rapidly enough. 8 Indeed, our results disagree with [64] for the central values by an amount which is larger then the theory error bars in the tail region as can be seen in figure 11. The reasons for this overshoot, as discussed above, is that if one does not use a profile to shut off the resummation beyond the transition region, then the scale of the log in the large p ⊥ region is such that the singular and non-singular terms in the fixed order cross section no longer cancel as they do in the fixed order result. This effect is a recurring theme in the literature [56,60,61]. The authors of [64] also include the two loop matching to the hard function and show that its effects are at the one percent level. Moreover, without the other pieces of the calculation of this order, 4-loop cusp, 3 loop non-cusp, two loop soft and collinear matching, this inclusion does not increase the accuracy of the calculation. Regarding the scale variation, the authors vary three distinct scales as was performed in this paper. The variations mix rapidity and renormalization group scale dependence (see appendix) and one could argue that this is perhaps not as clean as varying µ and ν independently. But given that the scale variation technique is a blunt instrument to estimate perturbative uncertainty, it is not clear that distinguishing logs is quantitatively relevant, as long as all resummed logs are probed in the variation. The uncertainty bands in [64] are very close to those presented here. In the appendix we show how the work of Grazzini et al. can be parsed in terms of the RRG. Finally, a recent paper [66] working within the TMDPDF formalism of [36,37,39] published a cross-section for the Higgs transverse momentum distribution at NNLL, with no fixed order matching. The error bands are very small, with no overlap between the NLL and the NNLL resummation, which the authors noted as a sign that they had under-estimated the perturbative uncertainty. The bands were generated by varying the resummation scales generated by the µ evolution of the TMDPDFs, after the resummation of the rapidity logs, which followed the procedure of [68]. The rapidity resummation scales were taken as fixed quantities. Thus no attempt was made to directly gauge the perturbative uncertainty in the rapidity resummation, that is, the perturbative uncertainty in the Collins-Soper kernel, whose exponentiation also resums the rapidity logarithms. Interestingly, variation of the non-perturbative models for the Collins-Soper kernel produced a much larger impact on the cross-section, see figure 4 of [66]. Since these non-perturbative models are explicitly exponentiated with the perturbative Collins-Soper kernel, in our view, this variation of the non-perturbative model parameters then may not be an estimate of the non-perturbative physics at all, but a more indirect gauge of the perturbative truncation of the Collins-Soper kernel and the rapidity resummation. Separating out the effects of the perturbative uncertainty in the Collins-Soper kernel, which has not been studied in the literature, will prove important in isolating the non-perturbative physics of the TMDPDFs.

Conclusion
The next LHC runs promises to shed more light on the nature of the Higgs boson. Here we have attempted to give the most accurate result possible for the transverse Higgs spectrum 8 There exists similar prescriptions to (8.1) that allow for turning off the resummation more quickly.
These have been applied to the 0-jet bin distributions of [15,16], however, to our knowledge have not been used in a published Higgs transverse momentum spectrum.

JHEP12(2015)097
using the currently available theoretical calculations for matching coefficients, anomalous dimensions as well as the fixed order cross section. We resum the large logs in the low p ⊥ region to NNLL and match this result to the fixed order cross section for high values of transverse momentum. To maintain accuracy over the complete range of the transverse momentum, we have utilized profiles in both the virtuality (µ) and rapidity (ν) factorization scales which smoothly turns off the resummation and matches onto the NNLO fixed order cross section at large values of p ⊥ . This procedure prevents the problem of enhancing the fixed order result. We discuss the sources of error and ways to effectively estimate them using possible variations in the rapidity and renormalization scales for the profiles introduced. While we get a good agreement with previous results in the low p ⊥ region, we get a lower estimate for the cross section in the transition and tail region.
The code which generated the plots in this paper is available upon request.

A RRG resummation
In this appendix, we review the RRG formalism for the purpose of explicitly comparing to the method used in [43,44,65]. After integrating over the rapidity of the Higgs, the differential cross section is given by where we have written out the explicit momentum fraction dependence. Since we are focusing on the resummation, we will often suppress the other arguments of the renormalized functions. The soft and beam functions have the µ-anomalous dimension: Using the fact that the RG and RRG variations of the amplitude must commute, we can write the rapidity anomalous dimension as integral over the µ anomalous dimensions [26,27] and a constant piece Writing things in this way allows for there to be large µ dependent logs which may arise depending upon the choice of µ. As noted in section 4, Γ R is related to the cusp anomalous dimension. The non-cusp piece of the rapidity anomalous dimension γ R α s ( be γ E 2 −1 ) can be extracted from the calculation of the soft function at the renormalization point µ = be γ E 2 −1 . This is a choice of scheme which defines the scale in the Log which multiplies the cusp anomalous dimension. We run the beam/soft functions in ν from ν s and m h to ν respectively Note that we have explicitly chosen ν B = m h in the TMDPDF. The full RRG factor has the rapidity factorization scale ν cancel in the exponent, and gives: We also need the µ evolution of the TMDPDFs at ν = m h and the soft function at ν = ν s . With these scale choices, the TMDPDF has no double log-µ evolution, so one obtains: Combining these pieces together, we achieve the full resummation kernel for the crosssection: From the form of (A.11), we can directly see the variation of the independent resummation scales µ, µ i , ν s probe different exponentiated logs, thus giving estimates of the subleading terms in the perturbative expansions. The double logarithmic terms associated with µ i cancel manifestly in the exponent, 9 and thus the µ i variation is estimating the subleading JHEP12(2015)097 terms associated with γ s,f ⊥ . We can eliminate large logs of the impact parameter from the beam and soft function matrix elements by choosing the scales µ i , ν s appropriately. Four potential schemes for canonical scale choices are: As noted in [69], making scale choices in either momentum space or conjugate space can have a sizable numerical impact on the cross-section, as well as how accurately the resummation captures the higher order logs claimed by the resummation. Scale setting in momentum space, residual logs can exist in the plus distributions of the low scale theory that become apparent only after integrating the resummation kernel against the matrix elements, in this case, the TMDPDF and the soft function. Scale setting in conjugate space does not suffer from this ambiguity. Specifically, for the canonical choice ν s = µ i = be γ E 2 −1 , R becomes: (A.12) A.1 Comparision to Grazzini et al.
We now write out the resummation formalism used in [64], as derived in [43,44,65], suppressing the flavor sum in the PDFs: dσ dp 2 Then one switches to moment space by integrating 1 z = x 1 x 2 for W: The authors write this as Note that the factorization scale from the PDFs is in the function H. Then the form of the resummed exponent is:

JHEP12(2015)097
Note that B N has PDF running in it since it depends on the moment N . The function H N is claimed to have no large logarithms, and can be perturbatively calculated. To one loop: H N m H , µ R , µ F ; Q = 1 + α s π H (1) + 2C (1) C N is related to the order pieces of the DGLAP splitting kernels. H (1) is the hard matching, but also contain terms in the δ(1 − z) pieces of the TMDPDF. Note that µ R ∼ µ F ∼ Q ∼ m h , so that the PDFs are run from Λ QCD to the high scale, implicitly. By choosing the appropriate scales, we can connect the resummed formula of (A.16) to the resummation formulas of (A.11) and (A.12). We must also apply the resummation scale prescription of [65], where all large logs are split: and any log of ln m h Q is expanded out of the exponent and included in the hard function of (A.15). Setting µ = Q in (A.12), this gives: choosing Q = µ this result is identical to (A.12). The hard function of eq. (A.15) is explicitly given to two loops in [65] , and one can verify the factor ln m h Q γ R [α s be γ E 2 −1 ] in (A.23) is expanded out of the exponent and then included in the B (2) term of the hard matching. Indeed, from the form of the hard function given in [65], one can see the resummation scale Q probes both rapidity logs and RG logs since B (2) includes both types of anomalous dimensions. This explains the similar sized perturbative error estimation of this paper, and that of [64].

JHEP12(2015)097
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.