Resummation Improved Rapidity Spectrum for Gluon Fusion Higgs Production

Gluon-induced processes such as Higgs production typically exhibit large perturbative corrections. These partially arise from large virtual corrections to the gluon form factor, which at timelike momentum transfer contains Sudakov logarithms evaluated at negative arguments $\ln^2(-1) = -\pi^2$. It has been observed that resumming these terms in the timelike form factor leads to a much improved perturbative convergence for the total cross section. We discuss how to consistently incorporate the resummed form factor into the perturbative predictions for generic cross sections differential in the Born kinematics, including in particular the Higgs rapidity spectrum. We verify that this indeed improves the perturbative convergence, leading to smaller and more reliable perturbative uncertainties, and that this is not affected by cancellations between resummed and unresummed contributions. Combining both fixed-order and resummation uncertainties, the perturbative uncertainty for the total cross section at N$^3$LO$+$N$^3$LL$^\prime_\varphi$ is about a factor of two smaller than at N$^3$LO. The perturbative uncertainty of the rapidity spectrum at NNLO$+$NNLL$^\prime_\varphi$ is similarly reduced compared to NNLO. We also study the analogous resummation for quark-induced processes, namely Higgs production through bottom quark annihilation and the Drell-Yan rapidity spectrum. For the former the resummation leads to a small improvement, while for the latter it confirms the already small uncertainties of the fixed-order predictions.


Introduction
After the discovery of the Higgs boson [1,2], the LHC has entered an era of precision Higgs measurements. One important goal is the precise determination of the Higgs couplings in order to test the Standard Model and search for evidence of physics beyond it. Other important color-singlet processes like Drell-Yan production serve as standard candles that are used, for example, to constrain parton distribution functions (PDFs). In order to match the ever increasing level of experimental precision, precise theoretical predictions for the measured cross sections are needed. An important example is the dominant Higgs production via gluon fusion, which receives large perturbative corrections. This has led to the calculation of the total production cross section up to N 3 LO [3][4][5][6][7][8][9][10], and including the resummation of threshold logarithms up to N 3 LL [11][12][13][14][15][16]. However, due to the limited detector acceptance the experimental measurements cannot measure the cross section fully inclusively but only in a restricted kinematic range, in particular in a restricted range of Higgs rapidities. The interpretation of the experimental measurements thus fundamentally requires theoretical predictions differential in the Higgs kinematics. The essential nontrivial ingredient is the Higgs rapidity spectrum (or equivalently the cross section with a rapidity cut), which is so far known to NNLO [17][18][19][20][21].
A specific class of perturbative corrections to Drell-Yan-like color-singlet production arises from the associated quark and gluon form factors, which contain Sudakov logarithms ln 2 (−q 2 /µ 2 ), where q µ is the transferred hard momentum. For spacelike momentum transfer, q 2 = −Q 2 < 0 as in deep-inelastic scattering, these logarithms vanish with the standard choice µ 2 = Q 2 . For timelike production processes, the form factor enters the production cross section evaluated at timelike momentum transfer q 2 = Q 2 > 0. With the ordinary scale choice µ 2 = Q 2 , the form factor contains leftover Sudakov logarithms ln 2 (−1) = −π 2 , inducing large corrections at each order in the perturbative series. For simplicity, we will henceforth refer to these as "timelike" logarithms or contributions, as they arise in the ratio of the timelike and spacelike form factors. 1 This effect was first observed long ago in Drell-Yan production in ref. [22], where it was realized that the coefficients of these terms are directly related to infrared (IR) singularities. Due to the universal structure of IR singularities, these terms arise to all orders and their resummation is well known [23][24][25][26]. As discussed in ref. [27], the timelike logarithms are also present in the soft contributions to the pion electromagnetic form factor providing an enhancement compared to the spacelike case in agreement with the measured enhancement. The resummation of the timelike logarithms for gluon-fusion Higgs production was carried out in refs. [11,28] in the context of soft-gluon (threshold) resummation, where it was shown that it substantially reduces the large perturbative corrections to the total gg → H cross section.
The resummation of the timelike logarithms originating in the form factors has since been included in the resummation of various other exclusive color-singlet cross sections (see e.g. refs. [29][30][31][32][33][34][35][36][37]), leading to improvements in the perturbative uncertainties. In these contexts, the use of the resummed form factor is unambiguous, as it explicitly appears as an ingredient in the corresponding factorized cross section.
In this paper, we study in detail the utility of the resummed timelike form factors for predictions of inclusive color-singlet production cross sections. In the case of inclusive cross sections the benefit of the resummation is a priori not obvious, and its applicability has occasionally been called into question. For this reason, we discuss in some detail the arguments for it and its consistent application, as well as the potential pitfalls one might worry about. For our numerical analysis, we consider both gluon-induced and quark- 1 Since the resummed logarithms ln 2n (−1) happen to give factors of (−π 2 ) n , their resummation has been referred to as "π 2 -resummation". Since factors of π 2 from other (unrelated) sources are typical to appear in the perturbative coefficients as well, we will always refer to the resummed logarithms as "timelike logarithms", to avoid any possible confusion as to what is being resummed. induced processes. The cases we consider include a generic scalar resonance gg → X as a function of m X , gg → H as a function of the Higgs rapidity, bb → H, and Drell-Yan qq → Z as a function of the Z rapidity.
We find that in all cases the resummation of the timelike logarithms leads to stable perturbative predictions. For the gluon-induced cases it leads to a significantly improved convergence compared to the fixed-order predictions, as first pointed out in refs. [11,28]. This results in perturbative uncertainties that are both smaller and more reliable. In addition to the total cross section studied previously, we show how the resummation can be easily and consistently applied to generic inclusive cross sections differential in the Born kinematics. This allows us in particular to obtain the currently most precise predictions for the Higgs rapidity spectrum, or equivalently the inclusive cross section with a rapidity cut, with perturbative uncertainties that are reduced by almost a factor of two compared to NNLO. For the quark-induced processes, the improvement is not as dramatic. Here, the resummed and fixed-order results have a similar stability. With an optimal choice of µ F the resummation still provides some improvement in the perturbative convergence and uncertainties. This demonstrates that using the resummed form factor is also viable for quark-induced processes and provides additional confidence in the estimated perturbative uncertainties.
The remainder of the paper is structured as follows: The basic setup how to consistently incorporate the resummed form factors into the inclusive cross section is discussed in sec. 2. The application to gluon-fusion processes is then discussed in sec. 3, to Higgs production through bottom quark annihilation in sec. 4.1, and to Drell-Yan production in sec. 4.2. We conclude in sec. 5. For completeness all required perturbative ingredients for the resummed form factors are collected in appendix A.

Resummation framework
We consider the hadronic production gg → L or qq → L of a color-singlet final state L with total invariant mass Q 2 = q 2 > 0. The hard virtual corrections to these processes are described by the corresponding QCD form factors. The full form factors contain infrared divergences, which when combined into the full cross section cancel against the infrared divergences in the real corrections. Hence, what enters in the final cross section are the IR-finite parts of the form factor. In the context of soft-collinear effective theory (SCET) [38][39][40][41], these are equivalent to the Wilson coefficients from matching the QCD currents defining the form factors onto the corresponding SCET currents [42][43][44]. For the cases we consider, these are the gluon, quark vector, and quark scalar form factors. The corresponding matching conditions read schematically where the B n⊥ and χ n are collinear gluon and quark fields in SCET. (The exact matching conditions for the currents can be found e.g. in refs. [30,45].) The IR divergences in the full QCD form factors, given by the quark and gluon matrix elements of the left-hand side, are exactly reproduced by the corresponding matrix elements of the SCET operators on the right-hand side, such that the hard Wilson coefficients C ij are given in terms of the IR-finite parts of the form factors. The relevant object entering the cross section is the hard function given by the square of the Wilson coefficient, which we write as where by default we normalize H to unity at leading order, and H (n) denotes the O(α n s ) term. To all orders in perturbation theory, C and H depend on the hard momentum transfer q µ through logarithms L ≡ ln[(−q 2 − i0)/µ 2 ]. For spacelike processes, q 2 = −Q 2 < 0 such that L = ln(Q 2 /µ 2 ) = 2 ln(Q/µ), while for timelike processes q 2 = Q 2 > 0 such that L = 2 ln(−iQ/µ).
The Wilson coefficients in SCET obey the renormalization group equation (RGE) where Γ cusp (α s ) is the cusp anomalous dimension and γ H (α s ) the noncusp term. Integrating eq. (2.3) yields the solution (2.5) The explicit result for the evolution kernel U H is given in appendix A.4. By choosing the imaginary-valued scale µ H = −iQ, the hard function H(Q, µ H ) is free of logarithms and can be calculated in fixed-order perturbation theory, while the evolution kernel U H resums all logarithms ln(µ H /µ) = ln(−iQ/µ). The hard function explicitly appears in calculations of exclusive cross sections as Here T denotes a resolution variable, which resolves additional emissions, such that in the limit T Q the cross section is restricted to the soft-collinear regime. In this limit it is dominated by hard virtual corrections contained in H, and soft and collinear contributions (both real and virtual) at lower scales µ T ∼ T contained in SC, while hard real emissions are forbidden. At the partonic level, an example for T is the partonic threshold variable (1 − z)Q. More physical examples of T are beam thrust or the p T of the leading jet. The precise form of the soft-collinear contribution SC depends on the definition of T but is irrelevant for our discussion. For a given process always the same hard function appears independently of the precise choice of T . The factorization in eq. (2.6) implies that in the T Q limit H appears as a well-defined perturbative object (namely as a hard matching coefficient), which is fully factorized from the rest of the cross section. In particular, the only dependence on the hard timelike momentum transfer Q 2 resides in H, while SC only depends on parametrically smaller soft and collinear scales proportional to T . In practice, eq. (2.6) can be used to perform the resummation of logarithms of T in dσ/dT , which involves using eq. (2.4) to evolve H from its natural scale µ H = −iQ to the relevant lower scale µ T ∼ T .
We want to apply the resummed form factor to the inclusive cross section for colorsinglet production. Here, inclusive refers to the fact that the cross section is fully integrated over any additional QCD emissions, but it can still be differential in or contain cuts on any kinematic variables that are present at Born level and describe the produced color-singlet system, such as its total rapidity Y or total invariant mass Q. To do so, we can factor out the hard function from the inclusive cross section which defines the remainder R(X, µ FO ). Here, X denotes any dependence on Born variables or cuts. By definition, H only depends on the Born kinematics via Q, while the remainder R can depend on X.
We write the perturbative expansion of the remainder as where for convenience we pulled out the leading-order cross section σ (0) (X, µ FO ). The dependence on the factorization scale µ F related to the PDFs entirely cancels within R, and we will mostly suppress it. The µ FO scale in eqs. (2.7) and (2.8) is equivalent to the renormalization scale µ R in the fixed-order prediction, and its dependence explicitly cancels between H and R. The R (n) coefficients depend primarily only on the total colorsinglet invariant mass and rapidity, while any dependence on additional Born kinematics or cuts resides primarily in σ (0) . (This becomes exact for a scalar resonance in the narrowwidth approximation like the Higgs.) In the following we will for simplicity suppress the dependence on X and Q 2 . We also define the K factor which captures the total perturbative correction relative to the leading-order result. Expanding eq. (2.7) order by order in α s (µ), it is straightforward to obtain the fixed-order coefficients of R from those of K and H. Up to N 3 LO we have, (2.10) To resum the timelike logarithms from the form factor in the cross section we can simply take the resummed result for the hard function eq. (2.4) and use it in eq. (2.7), As indicated, the fixed-order expansions for H(µ H ) and R(µ FO ) are reexpanded against each other (but without reexpanding the α s (µ H ) inside the coefficients H (n) (µ H ) in terms of α s (µ FO )). This is analogous to the standard treatment in resummed predictions as would be used for example in eq. (2.6). This ensures that in the limit µ H = µ FO we exactly recover the usual fixed-order result without inducing any higher-order cross terms between H and R. Using the definition of R in eq. (2.7), the resummed cross section in eq. (2.11) can equivalently be written as where the brackets [. . . ] FO indicate the fixed-order reexpansion in powers of α s (µ FO ) and α s (µ H ), with σ FO the usual fixed-order cross section expanded in α s (µ FO ). Written in this way, the ratio of timelike to spacelike form factors is manifest. Equation (2.11) will be the basis of all our results. For consistency with the fixedorder limit, we always include H(µ H ) and R(µ FO ) to the same fixed order. Furthermore, we always combine the N n LO fixed-order contributions with the N n LL resummation for H, which corresponds to the primed resummation counting and ensures consistency with the exclusive resummations [30,33] based on eq. (2.6). We will denote the perturbative accuracy by N n LO+N n LL ϕ , where the subscript indicates that the resummed logarithms correspond to the complex phase ϕ of the hard scale in the form factor.
While the remainder R is uniquely defined by eq. (2.7), one should of course ask the question to what extent it is justified or meaningful to "brute-force" factorize the perturbative series for the inclusive cross section into those for H and R.
First, one might be worried by the fact that the remaining nonlogarithmic constant terms in the fixed-order expansion of H(µ H ) are scheme-dependent, i.e. they depend on the fact that H is renormalized in the MS scheme and using a different scheme would result in different constant terms. However, this fixed-order scheme dependence is canceled by R up to higher orders, and this cancellation is explicitly ensured in our implementation in eq. (2.11) by the fact that we always reproduce the exact fixed-order result, as discussed above. The cancellation can also be seen explicitly from eq. Therefore, the relevant question is whether the series of timelike Sudakov logarithms present in H can be considered to be independent from the perturbative series in R. This would not be the case if (and only if) R were to contain contributions at each order correlated with the timelike Sudakov series in H and of opposite sign, which would then lead to large cancellations between H and R at each order in perturbation theory. These cancellations would then be spoiled by resumming the timelike logarithms in H while keeping the corresponding pieces in R at fixed order. This would imply that the perturbative corrections for R would be noticeably larger than for the cross section itself, and since the resummation of H eliminates its large corrections, the larger perturbative corrections of R would result in the resummed cross section being worse behaved. In other words, the absence or presence of sizeable cancellations between the resummed terms and the unresummed fixed-order terms, is mathematically equivalent to whether the resummation improves the perturbative convergence of the cross section or not. This is of course easy to check up to the available order, and in all our applications we have verified that there are indeed no large cancellations that are being spoiled by the resummation.
The primary reason one could be worried about such cancellations is that this is actually what happens in the reverse timelike process, namely color-singlet decays such as H → gg, H → bb, e + e − → Z → qq, or hadronic τ decays. These processes involve the same timelike form factor, but their perturbative series is known to not contain timelike Sudakov logarithms. The relation to these processes was already discussed in some detail in ref. [11]. In these processes, timelike logarithms only appear as single logarithms (and thus only at higher orders) through the running of α s , for which analogous analytic continuation methods have been considered, e.g. for e + e − → hadrons in refs. [27,[46][47][48][49] and hadronic τ -decays e.g. in refs. [50][51][52] (see also refs. [53,54] and references therein).
However, the situation is fundamentally different when the hard partons appear in the initial vs. the final state. An explicit discussion how the timelike Sudakov logarithms cancel in the final-state case but not in the initial-state case can be found in ref. [27]. In the final-state case, the process can be written as the imaginary part of forward matrix elements summed over all possible cuts, in which case the whole calculation can be deformed into the Euclidean domain where the timelike logarithms never appear. That is, the timelike Sudakov logarithms fully cancel between all cuts, or equivalently between the virtual corrections to the form factor and the real corrections to the corresponding remainder. The same does not happen if the partons appear in the initial state, which simply cannot be obtained from cutting a diagram, i.e. the process with incoming partons is intrinsically more exclusive, which exposes the timelike Sudakov logarithms in the form factor. Note also that if the same cancellations as in the final-state case were present in the initial-state case, they would have to be present at each order starting at NLO. The fact that we do not observe this even in the first several orders of the perturbative series provides clear evidence that this is indeed not the case.
It is also easy to understand why one finds a substantial numerical improvement for inclusive Higgs production. Comparing eq. (2.7) with the exclusive cross section in eq. (2.6), in the soft-collinear limit the remainder R reduces to the soft-collinear contributions times power corrections, (2.14) Therefore, the factorization in eq. (2.7) also becomes formally justified when the inclusive cross section is numerically dominated by soft-collinear contributions. It is well known that a large portion of the Higgs cross section comes from the partonic threshold limit, in which the hard function factors out of the cross section as in eq. (2.6). One can also take the more physical limit and simply veto additional hard radiation (which is also a weaker limit as it allows both soft and collinear radiation). Going from this exclusive 0-jet region, to which eq. (2.6) strictly applies, to the inclusive cross section amounts to factoring out the form factor also from the nonsingular power corrections. As pointed out in refs. [30,33], using either beam thrust or the p T of the leading jet to veto hard radiation, one finds that utilizing the resummed form factor for both singular and nonsingular corrections, and hence for the full inclusive cross section, is actually important, since not doing so can easily lead to unphysical results with the inclusive cross section being smaller than the 0-jet cross section.
Finally, we note that it has been argued in ref. [10] on the basis of the coefficient of the δ(1 − z) term in the partonic cross section that the timelike logarithms are not a dominant source of higher-order corrections and in particular that their resummation fails to improve the results beyond NNLO. We need to disagree with this assessment, because this coefficient is strongly scheme dependent and not a very well-defined quantity. Rather the impact or improvement should be judged at the level of the physical cross section. A more detailed discussion on this is given in appendix B.

Perturbative uncertainties and numerical inputs
For our numerical predictions we consider the LHC at E cm = 13 TeV. We use the PDF4LHC nnlo 100 [55][56][57][58][59][60] NNLO PDFs with α s (m Z ) = 0.118. Since we are interested in the size of the coefficients in the perturbative series we always use this same PDF independent of the perturbative order in consideration. The numerical value of α s (µ R ) is obtained with the corresponding three-loop running, except for the total gluon-fusion cross section known at N 3 LO, where we use four-loop running (though the numerical differences are negligible). For bottom-quark annihilation we use the PDF sets from refs. [61,62], which are reevolved from PDF4LHC nnlo mc in order to allow varying the b-quark matching scale separately from the b-quark mass. The relevant masses entering our predictions are m H = 125 GeV, m t = 172.5 GeV, m b (m b ) = 4.18 GeV, and m Z = 91.1876 GeV.
Since we are primarily interested in investigating the perturbative structure, we do not consider parametric uncertainties due to PDFs and the value of α s (m Z ), which are straightforward to evaluate. They are essentially unaffected by the resummation of the form factor, since all PDF dependence, as well as the dominant overall dependence on α s (m Z ) in case of Higgs production, resides in the remainder R.
An important aspect of precision predictions is a reliable assessment of the theory uncertainties due to missing higher-order corrections. Our predictions in principle involve three scales that we can vary as a means to estimate the size of higher-order corrections: the factorization scale µ F probing collinear logarithms in the PDFs, the renormalization scale µ R probing higher orders in the fixed-order series, and the hard resummation scale µ H probing higher orders in the series of timelike Sudakov logarithms. We like to stress that these scales are unphysical parameters whose variations simply provide a convenient way to probe the "typical" size of the associated missing higher-order terms. The resulting variations in the cross section must be interpreted as such. In particular, we do not assign any meaning to accidentally small one-sided scale variations that yield asymmetric uncertainties, which are just artifacts of a nonlinear scale dependence, which is frequently encountered in predictions at higher orders or involving resummation. We therefore always consider the maximum absolute deviation from the central result at the chosen central scale as the (symmetric) uncertainty. To be explicit, an observed scale variation of +|x| and −|y| in the cross section is interpreted as a perturbative uncertainty of ± max{|x|, |y|}.
We parametrize the three scales as The choices for µ FO and κ F for the central value depend on the process we consider. For the resummed predictions we use the central choice ϕ = π/2 , while the fixed-order predictions correspond to taking ϕ = 0, which turns off the resummation. We explicitly distinguish two different sources of perturbative uncertainties, namely fixed-order and resummation uncertainties, that are associated to the two independent perturbative series involved. The fixed-order uncertainty, denoted as ∆ µ , is obtained via the conventional variations of µ R and µ F . This comprises a collective overall variation of µ FO by a factor of two around its central value, which is combined with an additional variation of κ F by a factor of two around its central value, without considering the extreme variations where both are varied up or down at the same time. That is, relative to the central values we consider the set of variations from which the fixed-order uncertainty ∆ µ is obtained as the maximum deviation from the central value In the limit where the resummation is turned off, this reproduces the perturbative uncertainty in the fixed-order predictions. For the resummed predictions, the magnitude of the hard scale by construction follows the µ FO variation, |µ H | = µ FO , as illustrated in fig. 1 on the left, such that the fixed-order variations do not change the resummed logarithms ln(µ H /µ FO ). For the resummation uncertainty, we vary the phase ϕ in the interval [π/4, 3π/4] around the central value of ϕ = π/2, while keeping µ FO at its central value, as illustrated in fig. 1 on the right. This probes the intrinsic size of the higher-order timelike logarithms. The phase variation by ±π/4 is chosen to be roughly equivalent to the usual factor of 2 for conventional logarithms since π/4 ln 2. The uncertainty ∆ ϕ is then obtained as Figure 1. Illustration of the scale variations used to estimate the perturbative uncertainties. Left: The overall variations of µ FO , which determines ∆ µ (in conjunction with the variation of κ F , which is not shown). Right: The phase variation for µ H for fixed µ FO , which determines the resummation uncertainty ∆ ϕ .
the maximum observed deviation from the central value (usually happening at one of the endpoints), such that This additional resummation uncertainty was not considered in earlier treatments, but has already been included in the resummed 0/1/2-jet-bin results reported in ref. [63]. The total perturbative uncertainty is obtained by adding the two independent sources, For bb → H we follow ref. [62] and consider the low-scale matching at µ b onto the b-quark PDFs as a third independent source of uncertainty ∆ b , which is estimated by varying µ b by a factor of two.

Gluon fusion
Gluon-fusion processes are well-known to contain large perturbative corrections, which are partially due to the timelike logarithms in the gluon form factor, as first demonstrated in ref. [28]. We first consider the total production cross section up to N 3 LO+N 3 LL ϕ for a generic scalar final state gg → X in sec. 3.1 and for the SM Higgs boson in the rEFT m t → ∞ limit in sec. 3.2. In sec. 3.3 we discuss how to incorporate quark-mass and electroweak effects into the resummed results. In sec. 3.4 we then present our results for the Higgs rapidity spectrum and the cross section with a rapidity cut to NNLO+NNLL ϕ .

Color-singlet production
We first consider the total production cross section from gluon fusion for a generic colorsinglet scalar X with mass m X . Its coupling to gluons at the scale µ ∼ m X can be expressed in terms of an effective Lagrangian as where Λ is a suitable high mass scale and C X is the Wilson coefficient from integrating out heavy particles that mediate the effective ggX interaction. This effective operator arises for SM Higgs production in the m t → ∞ limit, which we discuss in more detail in sec. 3.2.
Here, we use it as a simple case to study the effects of the resummation and its dependence on the mass over a wide range m X ∈ [100, 1000] GeV. For this purpose, the precise values of the effective coupling C X (µ = m X )/Λ need not be specified, as it drops out for the K-factor σ/σ (0) on which the resummation acts. We obtain the total gg → X cross section to N 3 LO from SusHi 1.6.0 [6, 10, 64-69]. Our central scale choices are µ FO = m X and κ F = 1, such that µ R = µ F = m X . Away from µ = m X , the perturbative running of C X (µ) induces logarithms of m X /µ at NNLO and N 3 LO. Their resummation is irrelevant and can be neglected, and they are instead included in the fixed-order cross section [66].
The gluon form factor is known up to three loops [70][71][72][73][74][75], and the Wilson coefficient C gg is explicitly extracted from it in ref. [75] (see also refs. [76,77]), where Γ g cusp (α s ) is the gluon cusp anomalous dimension and the last three terms are the total noncusp contribution. All the relevant ingredients are collected in appendix A.
The separation of the perturbative series for the K factor at fixed order into those of H and R is shown in fig. 2 as a function of m X . Half of the large NLO K factor comes from H and half from R, while beyond NLO the corrections in H are larger than for R. Hence, the large corrections to the K-factor present at each order are driven to a large extent (but also not entirely) by the corrections from H. In particular, the remainder R by itself has a much better behaved perturbative series than K, and there are clearly no cancellations between H and R. (Otherwise, as already explained in sec. 2.1, R would need to have negative corrections that are larger in size than those in K.) This pattern holds independently of m X . The visible increase in the corrections toward smaller m X is due to the running of α s (m X ).
The large perturbative corrections in H gg at the real scale µ H = m X are absent at the imaginary scale µ H = −im X , as shown by the long-dashed curve in the middle panel of fig. 2. To illustrate this more explicitly, the numerical values for an example mass of where each term is the contribution from a subsequent order in α s up to N 3 LO. Clearly, the large corrections to the gluon form factor at real scales are almost entirely due to the timelike Sudakov logarithms that are present for µ H = m X and are eliminated by taking µ H = −im X . Since the corrections in H gg at µ H = −im X are very small, the perturbative convergence of the resummed cross section will be essentially determined by that of the remainder R. In fig. 3, we compare the fixed-order and resummed cross sections as a function of m X , with the bands showing the total perturbative uncertainties evaluated as discussed in sec. 2.2. (Note that in case of gg → X and gg → H, the fixed-order uncertainties come from the variation of µ R for fixed µ F .) All results are normalized to the LO prediction σ (0) at fixed µ FO = m X . As expected, the absence of large corrections in the resummed hard function directly translates into a much faster convergence of the resummed cross section. Furthermore, the uncertainties in the resummed predictions at lower orders cover the higher-order bands much better than at fixed order, while at the same time being substantially reduced at higher orders. Hence, even at NNLO and N 3 LO, where the fixed-order results start to show convergence, the resummation noticeably improves the predictions. Due to their better convergence, the resummed predictions provide substantially improved The total cross section for gg → X at E cm = 13 TeV at fixed order (left) and including the resummation of timelike logarithms (right). All results are normalized to the central LO prediction at µ FO = m X . uncertainty estimates both in terms of their reliability and their size. In particular, we can be reasonably confident that the result at the next higher order will lie within the small N 3 LO+N 3 LL ϕ uncertainty band.

Inclusive Higgs production in the rEFT scheme
We now turn to the case of Higgs production through gluon fusion as an important application of the singlet production discussed above. For Higgs masses below the top threshold, m H < 2m t , the gluon-fusion cross section can be well approximated by an effective theory where the top quark is integrated out [78][79][80][81], giving rise to an effective Lagrangian analogous to eq. (3.1), In this case, the Wilson coefficient C t itself receives sizable QCD corrections, which have been calculated to N 4 LO in refs. [82][83][84]. The effective operator in eq. (3.5) is the same as in eq. (3.1), giving rise to the same gluon form factor and hard function H gg in eq. (3.2).
Rescaling the cross section σ EFT obtained from eq. (3.5) by the LO m t dependence [85] one obtains the inclusive cross section in the "rescaled EFT" scheme (rEFT), This rescaling is known to well reproduce the m t -exact result at NLO, and hence it is believed to be a useful approximation also at higher orders [5,[86][87][88][89][90][91][92]. The inclusion of further quark mass and electroweak effects will be discussed in sec. 3.3. We use SusHi 1.6.0 [6,[64][65][66][67] to compute the total cross section in the rEFT scheme to NNLO. For the N 3 LO contribution we use the results of ref. [10] as implemented in ggHiggs 3.5 [16]. 3  where again each term gives the contribution from a subsequent order in α s . The remainder R now includes the corrections to |C t | 2 . As before, its perturbative series is much better behaved than that of the cross section, whose large perturbative corrections are thus driven by the large corrections from timelike logarithms in H gg .
To illustrate the improved convergence of the resummed form factor, we consider the hard function H gg (m H , µ H ) at various scales µ H , Here, the uncertainty bands only show the fixed-order uncertainty ∆ µ . At ϕ = 0, σ res (ϕ) is just the fixed-order cross section. As ϕ → π/2, the timelike resummation is turned on, visibly improving the convergence of the cross section and providing better coverage of the uncertainty bands. The ϕ dependence becomes stationary at ϕ = π/2, where the timelike Sudakov logarithms exactly vanish. Beyond ϕ = π/2, powers of ϕ − π/2 start to enter again.
In the right panel of fig. 4, we compare the fixed-order results at the conventional scales of µ FO = m H and µ FO = m H /2 with the resummed results. The results are shown as relative corrections to our best prediction at N 3 LO+N 3 LL ϕ . For the resummed results, the inner uncertainty bars indicate ∆ ϕ alone, while the outer ones show ∆ µ ⊕ ∆ ϕ . While ∆ ϕ contributes to obtaining a more realistic uncertainty estimate at LO+LL ϕ (compared to LO), its impact is strongly reduced at higher orders. The overall picture and conclusions from the generic color-singlet case are unaffected by the presence of the Wilson coefficient |C t | 2 in the cross section. The resummation yields again a clear improvement in convergence and uncertainties, also compared to the fixed-order results at µ FO = m H /2, which are already better behaved than those at µ FO = m H . In particular, the NLO+NLL ϕ result already fully covers the highest-order result, which is not the case at fixed NLO, and the precision of the NNLO+NNLL ϕ result is roughly comparable to the fixed N 3 LO results. This gives us good confidence in the small remaining uncertainty at N 3 LO+N 3 LL ϕ , which is reduced by a factor of two compared to N 3 LO. The explicit numerical results at the highest order are variation in the resummed results gives ∆ µ = 0.67 pb, which combined with ∆ ϕ then yields a total perturbative uncertainty of 1.44%. Note also that using the threshold-expanded running in µ R and µ F as implemented in SusHi 1.6.0, the N 3 LO result at µ FO = m H /2 increases to (48.17±1.99 µ ) pb (4.14%), with a corresponding increase in ∆ µ since the result at µ FO = m H is unaffected.

Incorporating quark mass and electroweak effects beyond rEFT
While the previous section focused on the QCD corrections to Higgs production in the m t → ∞ limit, further corrections arise from finite quark-mass effects as well as electroweak contributions. Here we discuss how to consistently combine them with the resummation of timelike logarithms. The full dependence of the cross section on the heavy quark masses m t , m b , m c is fully known at NLO [5,64,86,87,[93][94][95]. We define δσ q (N)LO as the correction of the exact result relative to the rEFT result, On top of the exact NLO corrections, top-quark mass effects are also known in an asymptotic expansion in 1/m t at NNLO [88][89][90][91][92].
In the following we consider the top-mass effects in more detail. As discussed in refs. [66,[90][91][92], the asymptotic 1/m t corrections at NNLO cannot be expected to reliably improve over the m t → ∞ limit. Rather, they can serve to estimate the uncertainty due to the still unknown full NNLO m t corrections. For this reason we will only take into account the NLO corrections δσ t NLO . (The inclusion of the NNLO m t corrections would be completely analogous.) This is also consistent with our analysis of the rapidity spectrum in sec. 3.4, for which the m t -corrections are only known at NLO. For illustration, the numerical results for δσ t NLO are The finite m t contributions correspond to a correction to the m t → ∞ limit in eq. (3.5), from which the gluon form factor arises, and so a priori they do not involve the same local gluon form factor. Therefore, one option to include them in the resummed results is to simply add them to the rEFT results in eq. (3.10), which yields The complete results including those at lower orders are collected in table 1. Alternatively, following ref. [30] we can perform a one-step matching of the full Standard Model including the top quark onto SCET, simultaneously integrating out both the top quark and hard virtual corrections. The resulting hard function H t gg = |C t gg | 2 corresponds to the full SM gg → H form factor and includes all virtual finite-m t effects. It takes the form [30] H t gg (m t , m 2 where as before ρ ≡ m 2 H /(4m 2 t ). Compared to H gg = |C gg | 2 , the LO m t dependence F 0 (ρ) and the contributions from C t are now moved from the remainder into the hard function. The F 1 (ρ) contains the full virtual m t dependence at NLO and the O(ρ α 2 s ) terms denote the neglected NNLO virtual m t corrections. 5 Although H t gg is no longer normalized to unity at leading order, we can continue to use eq. (2.12) to obtain the resummed cross section. Compared to eq. (2.11), the result now contains an overall factor |α s (µ H )/α s (µ FO )| 2 from the ratio of hard functions, which replaces the α 2 s (µ FO ) inside the LO cross section by |α s (µ H )| 2 .
The RGE for C t gg is given by The noncusp terms in γ t gg differ from those in γ gg in eq. (3.3) due to the additional µ dependence of α s (µ)C t (µ), which is now included in the hard Wilson coefficient. The overall |α s (µ)| 2 |C t (m t , µ)| 2 in eq. (3.14) is now evaluated at µ H = −iµ FO and then evolved back to µ FO . For the overall α s (µ) this is largely irrelevant since it is ultimately evolved starting from α s (m Z ). For C t (µ), which is treated in fixed order, this induces different subleading timelike logarithms starting at NNLO compared to H gg . This is reflected in the noncusp terms differing by γ t , whose numerical effect however is not significant. Also, the perturbative convergence of |C t (µ)| 2 at µ = m H and µ = −im H (and at its natural scale µ = m t ) is practically the same.
The perturbative convergence of H t gg shows the same improvement as seen for H gg when evaluated at µ H = −im H rather than µ H = m H , The main difference compared to H gg are the additional constant terms from C t that are now included in H t gg . The finite-m t corrections have a very small effect on the NLO contribution, contributing a +0.005 to the above 0.82152 and 0.27631. 5 The −F1(0) here simply removes the leading mt → ∞ part of F1(ρ), which is already included via Ct.
We drop all cross terms of F1(ρ) − F1(0) with Cgg, which are of O(ρ α 2 s ) and higher, because these terms are also not included in the fixed-order cross section. This is equivalent to the results at µ H = −im H reported in ref. [63]. Including the full NLO m t dependence, we obtain The full set of results including the lower orders are shown in the last column of table 1. Comparing the last two columns of table 1, the resummed results using the two different ways to include the top-quark contributions are perfectly compatible with each other. The fixed-order uncertainty is essentially unaffected, because it is insensitive to the precise split of the constant terms into H and R due to the reexpansion of their fixed-order contributions [see eq. (2.11)]. The resummation uncertainty ∆ ϕ increases somewhat in the one-step matching, which reflects the fact that the C t contributions introduce an additional residual µ dependence and that they are evaluated at µ = −im H rather than their natural scale µ = m t . Overall, the numerical differences are however completely insignificant, which shows that the results are insensitive to the precise treatment of the top contributions. This also provides nontrivial verification that the scheme dependence in how the nonlogarithmic constant terms are split between H and R at each order is much smaller than the perturbative uncertainties and hence irrelevant.
A complete numerical inclusion of all known corrections beyond the rEFT limit is beyond the scope of this paper. The inclusion of b-quark and electroweak effects can proceed completely analogously to the treatment of the top contributions. Any multiplicative contributions can be trivially included, while additive corrections such as the NLO m bdependent terms can be treated analogously to the finite-m t corrections. For example, the dominant known electroweak corrections can be included by replacing [96] C t → C t + δ EW (1 + C 1w α s + · · · ) , where δ EW is the pure NLO electroweak correction to the LO cross section [97,98] and C 1w contains the mixed O(αα s ) correction calculated in ref. [96] by integrating out W -and Zbosons as an estimate of the full O(αα s ) corrections. These additional contributions will not affect the benefit of the resummation, in the same way the inclusion of the top corrections for gg → H did not affect the conclusions compared to the generic scalar gg → X case.

Higgs rapidity spectrum
As discussed in sec. 2.1, the resummed form factor can be incorporated in the same way as for the total production cross section into generic cross sections that are differential in or contain cuts on the Born kinematics. Here we consider the primary example of the rapidity spectrum as well as the cross section with a rapidity cut. For simplicity we do not consider additional fiducial cuts on the Higgs decay products here, but stress again that these are straightforward to include.
The rapidity spectrum for gluon-fusion Higgs production is known to NNLO [17][18][19][20][21], while the N 3 LO corrections are available in the threshold limit [99,100]. The resummation in the small-x limit is also known [101]. We obtain the fixed-order bin-integrated rapidity distribution for gg → H to NNLO with HNNLO 2.0 [20,21]. We use a binsize of ∆Y = 0.25 and for clarity in all plots interpolate the binned results.
We first consider the rEFT limit and exclude additional quark mass effects. In fig. 5, we display the perturbative remainder R(Y ) as a function of Y . Although it has some intrinsic nontrivial rapidity dependence, the overall behavior is as for the total cross section, namely it exhibits a noticeably better convergence than the full fixed-order spectrum. Hence, we expect a similar improvement from applying the resummation to the rapidity spectrum as for the total cross section.
The upper panels of fig. 6 show the fixed-order results at µ FO = m H and µ FO = m H /2, with the bands showing ∆ µ . The overall K factor at NLO and NNLO is roughly constant in the central rapidity range and similar to that of the total cross section. This is consistent with the fact that a large part of the K factor stems from the timelike logarithms in the gluon form factor, which is independent of the rapidity.
The resummed result including fixed-order and resummation uncertainties, ∆ µ ⊕∆ ϕ , is shown in the bottom panel of fig. 6. Clearly, resumming the timelike logarithms improves the perturbative convergence across the spectrum as it did for the total cross section. The NNLO+NNLL ϕ result has perturbative uncertainties that are almost a factor of two smaller than at NNLO. At the same time, the NNLO+NNLL ϕ result is well covered by the lower-order NLO+NLL ϕ uncertainty band, which is not the case at fixed order. Judging from the results for the total cross section, for which the full N 3 LO is known, we expect the precision of the NNLO+NNLL ϕ to be roughly comparable to what can be expected for the full N 3 LO result, and we can be confident that the corresponding N 3 LO+N 3 LL ϕ result will have further reduced uncertainties and will be well contained within the current NNLO+NNLL ϕ uncertainties.  Table 2. Cross section for gg → H with a rapidity cut of Y cut = 2.5. All results include the exact m t dependence δσ t (Y cut ) at NLO. The percent uncertainties for the resummed results correspond to the total uncertainty ∆ µ ⊕ ∆ ϕ . The approximate n = 3 results are constructed as explained in the text and are given for illustration only.
Next, we consider the cross section with a rapidity cut, We now also include the exact m t -dependence at NLO using HNNLO 2.0 [102], from which we extract the rapidity-dependent analogue of δσ t NLO in eq. (3.11). The latter is included in the resummed results as discussed in the previous subsection: It is either added to the resummed result based on H gg giving d(σ res + δσ t )/dY . Alternatively, we can apply the resummation to the full m t -exact spectrum d(σ rEFT + δσ t ) res /dY using H t gg . The obtained differential rapidity spectra look extremely similar to those in fig. 6, since the finite-m t correction yields a small negative shift −0.2 pb < d(δσ t )/dY < 0 throughout.
In table 2, we provide benchmark results for Y cut = 2.5 at fixed order and including resummation. As expected, the resummed results show the same improved convergence and perturbative uncertainties as the rapidity spectrum and the total cross section.
Since the full N 3 LO rapidity spectrum is not yet known, one might think about approximating it by rescaling the NNLO spectrum to the inclusive N 3 LO cross section, e.g., by taking σ N 3 LO (Y cut ) ≈ (σ N 3 LO /σ NNLO ) × σ NNLO (Y cut ). Although this likely improves the central value by moving it closer to the correct result, it is unclear what uncertainties one should assign. In fact, the resummed NNLO+NNLL ϕ result for the rapidity spectrum provides a clean way to essentially achieve this goal, because it includes the dominant part of the inclusive N 3 LO K factor, and it does so with the correct rapidity dependence and reliable uncertainties, which are already substantially reduced compared to NNLO.
To illustrate this, first note that the entire rapidity dependence is contained in the remainder R. Looking at fig. 5, one might assume that the N 3 LO contribution is roughly flat in rapidity such that taking provides a reasonable approximation, where R (3) /σ (0) is the N 3 LO correction for the total cross section. Using eq. (3.21) as input and combining it with the known H (3) and lower-order results we obtain the approximate numbers at N 3 LO and N 3 LO+N 3 LL ϕ shown in the last line of table 2. As expected, the NNLO+NNLL ϕ results indeed capture most of the expected shift from the NNLO result. Using the same procedure to evaluate the scale variations, we can project that the relative uncertainties will be similar to the inclusive cross section, namely around 4% at N 3 LO and around 2% at N 3 LO+N 3 LL ϕ . However, we caution again that it is debatable whether it is justified to assign these without knowing the exact result.

Quark annihilation
We now turn to qq annihilation processes, for which the perturbative corrections are typically much smaller than for gluon fusion. For these, timelike Sudakov logarithms still arise in the corresponding quark form factor at timelike momentum transfer and can be resummed to all orders. We apply the resummation to Higgs production through bottom quark annihilation in sec. 4.1 and to the Drell-Yan process in sec. 4.2.

Higgs production through bottom-quark annihilation
The cross section for Higgs production through bottom-quark annihilation, bb → H, is much smaller than for gg → H, but is important phenomenologically as it provides direct access to the bottom-quark Yukawa coupling and can be enhanced in theories beyond the Standard Model. The hard function for bb → H corresponds to the quark scalar form factor and is obtained from the SCET matching coefficient C S qq (Q 2 , µ) for the scalar current with two massless quarks, see eq. (2.1). The scalar form factor naturally arises in the five-flavor scheme calculation in which the b quark is treated as a massless quark at the hard matching scale, except for its Yukawa coupling y(µ) to the Higgs, which always has its physical value corresponding to m b (µ). We exclude the Yukawa coupling from the hard matching coefficient, just as we excluded the overall |α s C X | 2 for gluon fusion.
To the best of our knowledge, C S qq is not yet directly available in the literature. We extract it in appendix A.3.4 to N 3 LO from the three-loop massless quark scalar form factor calculated in ref. [103] (see also refs. [104][105][106] for the NNLO form factor, and ref. [107] for the NNLO form factor including the full mass dependence). The RGE of C S qq is given by where Γ q cusp is the quark cusp anomalous dimension, and the last two terms are the total noncusp contribution, where 2γ q C is the usual hard quark noncusp anomalous dimension (also appearing for the vector current) and γ m is the mass anomalous dimension of the Yukawa coupling.
We obtain the bb → H total cross section to NNLO in the five-flavor scheme [108,109] from SusHi 1.6.0 [6,[65][66][67]104]. The Yukawa coupling is evaluated at µ FO and obtained In contrast to the gluon-initiated case, here the µ F dependence of the cross section plays the dominant role, while the µ R dependence is much less important. For our central scales we use µ FO = m H and κ F = 1/4 corresponding to (µ R , µ F ) = (m H , m H /4), as typically adopted in the five-flavor scheme calculation. This low central value for µ F is motivated by the observation that it roughly minimizes the collinear logarithms in the partonic cross section and leads to more stable predictions, see e.g. refs. [61,104,[110][111][112][113][114]. This is also seen in the left panel of fig. 7, which shows the fixed-order results as a function of the central value used for κ F in the range κ F = µ F /m H ∈ [1/8, 1] (always using µ FO = µ R = m H for the central value). The bands show the fixed-order uncertainty ∆ µ , which is obtained from the µ FO and κ F variations as discussed in sec. 2.2.
The perturbative series of the fixed-order cross section and its separation into H S qq and R are given by   showing that the corrections at real scales µ H = m H are mostly due to the timelike logarithms, which are eliminated at µ H = −im H . The resummed results are shown in the right panel of fig. 7, again as a function of the central choice for κ F . From the above observations, we expect the resummed results to have an improved convergence at µ F = m H /4. This is clearly seen, as the resummed predictions at the different orders roughly intersect at µ F = m H /4, with the corrections beyond LO+LL ϕ amounting to less than 10%. The explicit numerical results at each order including a breakdown of the perturbative uncertainties are collected in table 3.
Note that the remainder R carries the full µ F dependence of the cross section. Evaluating it at κ F = 1, corresponding to µ F = m H , its corrections are substantially larger than at µ F = m H /4, We also checked that in the considered range of µ F the perturbative coefficients of the remainder are indeed minimized around µ F = m H /4. This explains why the resummed results in fig. 7 do not improve the fixed-order results above κ F = 1/4, since there the cross section becomes dominated by large negative corrections to the remainder and is also affected by accidental numerical cancellations between H and R. Hence, the central choice κ F = 1/4 is also quite optimal from this point of view, as it leads to a clear "division of labor" between the remainder and the hard function in capturing the perturbative corrections to bb → H: By evaluating the former at an appropriate µ F the collinear logarithms arising from initial-state g → bb splittings are resummed into the b-quark PDFs, while the latter captures the hard virtual contributions to the bbH vertex, the bulk of which are enhanced at timelike kinematics and are resummed by setting µ H = −im H . The resummation also reduces the fixed-order uncertainty ∆ µ by reducing the µ FO dependence and by eliminating large cross terms, which can also reduce the impact of the large µ F dependence from lower-order contributions. This is seen in the slightly reduced µ F dependence in fig. 7 at NLO+NLL ϕ and NNLO+NNLL ϕ compared to their fixedorder counterparts. However, the reduction of ∆ µ is not nearly as dramatic as for gg → H, because the µ F dependence plays a much bigger role here. (In principle, the µ F dependence and resulting uncertainties can be reduced further by reorganizing the perturbative series as discussed in ref. [61], which is however beyond our scope here.) In table 3 we also include the ∆ b uncertainty from the low-scale b-quark PDF. (The parametric uncertainty in m b is much smaller than ∆ b and not considered.) The resummation uncertainty ∆ ϕ is completely negligible compared to ∆ µ , thanks to the very stable resummed hard function.
Overall, we find that the NNLO and NNLO+NNLL ϕ results are very similar, which is reassuring, and that the resummation of timelike logarithms provides a useful tool to accelerate the convergence of the bb → H cross section and to reduce its perturbative uncertainties.
(For the inclusion of threshold resummation effects see e.g. refs. [24,[121][122][123][124][125][126].) The necessary quark vector form factor is known to three loops [71][72][73][74][75][127][128][129][130]. The corresponding hard matching coefficient C V qq for the quark vector current to N 3 LO was obtained in refs. [75,131]. We also need the matching coefficient C A qq for the axial-vector current, which is equal to the vector coefficient up to singlet corrections that start entering at O(α 2 s ) [45]. At NNLO the axial-vector coefficient receives a nonvanishing singlet contribution from the axial-vector anomaly due to the large bottom-top mass splitting, but these contributions have been found to be small at cross section level [132,133]. Since they are also not implemented in the program Vrap 0.9 [116], which we use for our fixed-order predictions, we set C A qq = C V qq , dropping any singlet terms, and take H qq ≡ |C V qq | 2 for our analysis. The RGE in either case is identical, since the axially anomalous terms are nonlogarithmic, and reads We consider the double differential cross section d 2 σ/dY dQ, where Y and Q are the rapidity and invariant mass of the produced lepton pair, and for simplicity focus on the Z pole, Q = m Z . Note that perturbative corrections to the quark current can only have a weak logarithmic dependence on Q, so our observations also hold away from the Z pole. We obtain the rapidity distribution from Vrap 0.9 [116], taking into account V = Z, γ * and their interference terms.
As for bottom-quark fusion, we can expect that for Drell-Yan the choice of µ F will be important to improve the predictions by minimizing the perturbative corrections to the remainder. We first consider the rapidity spectrum at fixed Y = 0. (Alternatively, one could integrate over rapidity, which yields similar conclusions.) Figure 8 shows the dependence on the factorization scale µ F at fixed order (top left panel) and including the  Table 4.
Cross section for pp → Z/γ * → + − at Q = m Z and Y = 0 at the LHC with E cm = 13 TeV. The central scale is always µ FO = µ R = m Z . The percent uncertainties for the resummed results correspond to the total uncertainty ∆ µ ⊕ ∆ ϕ . resummation of timelike logarithms in the quark form factor (top right panel), always using µ FO = µ R = m Z . In the bottom panel of fig. 8 we directly compare the NNLO (gray) and NNLO+NNLL ϕ (orange) results. The fixed-order results show the best convergence around µ F = m Z , which is the typical choice adopted for Drell-Yan. In contrast, the resummed results clearly favor lower values of µ F . We therefore take µ F = m Z /2 as our central choice, as in this region the lower-order uncertainties provide the best coverage of the higher-order results, while at the same time the µ F dependence at the highest order is the most stable. (In principle, one might even consider going as low as µ F ≈ m Z /4 close to where the NLO+NLL ϕ and NNLO+NNLL ϕ central values coincide. However, the µ F dependence is much larger there, and we also prefer to stay within a factor of two of the canonical value µ F = m Z .) The dot-dashed line in fig. 8 shows the results at the highest order only including the qq-initiated contributions. At µ F = m Z /2, the cross section is dominated by the qq contribution, while at µ F = m Z the total NNLO corrections are small only due to substantial cancellations between the qq and remaining channels. This supports our central choice of µ F = m Z /2 in the resummed results as the quark form factor is associated primarily with the qq channel. The explicit results at Y = 0 are given in table 4, including a breakdown of the uncertainties.
The perturbative series of the fixed-order cross section at Q = m Z and Y = 0 and its decomposition into hard function H qq and remainder R at µ F = m Z and µ F = m Z /2 is given by  In the bottom panel, the highest-order results at NNLO and NNLO+NNLL ϕ are directly compared with a further zoomed-in y-axis. The uncertainty bands show ∆ µ (fixed order) and ∆ µ ⊕ ∆ ϕ (resummed). discussed in sec. 2, there is a priori no reason to expect that this continues to happen at higher orders. Also, the NLO and NNLO contributions for R are of the same size, indicating that the NLO contribution is artificially small or the NNLO contribution unusually large or a mixture of both. For µ F = m Z /2, the NNLO contribution to R is very small and the NNLO contribution to the cross section primarily comes from H qq . It will be interesting to see how this pattern continues at N 3 LO.
The hard function itself shows again a notably improved convergence at µ H = −im Z compared to µ H = m Z , As for bbH, the improvement is not as dramatic as for gluon fusion due to the reduced color factor for quarks vs. gluons. Nevertheless, by eliminating the timelike logarithms in H qq , the higher-order contributions are substantially reduced, except that the NLO contribution actually gets larger (in contrast to bbH). The reason for this is an accidental numerical cancellation in the one-loop matching coefficient 8) where the rather large nonlogarithmic constant term of −8 + π 2 /6 partially cancels the − ln 2 (−1) = π 2 when H qq is evaluated at µ H = m Z . As discussed in sec. 2, the separation of the nonlogarithmic constant terms between H and R amounts to a scheme choice and only their sum is ultimately relevant. Hence, this large NLO constant term is a scheme-dependent artifact and in fact cancels most of the equally large NLO contribution in R(µ F = m Z /2). This also means that the +0.263 NLO contribution in the cross section at µ F = m Z /2 in eq. (4.6) does in fact primarily come from the NLO timelike logarithm in H qq , which gives a contribution of +0.247 to it, even though this is not immediately obvious from eq. (4.6). This explains the much improved convergence of the resummed results at µ F = m Z /2 compared to the fixed-order results at the same scale.
It was already noted in ref. [25] that the constant terms in C V qq are scheme dependent and hence not physical, unlike the ratio of form factors. Since the constant terms in H and R are evaluated at different scales, there is a residual scheme dependence, which is analogous to a scale choice in that it affects the numerical results but is formally of higher order [see eq. (2.13)]. To check this, we can consider an alternative renormalization scheme for the Wilson coefficientC V qq , for which all constant terms exactly vanish. That is, the corresponding hard function solely consists of timelike logarithms,H V qq (m 2 Z , µ H = m Z ) = 1 + 0.247 + 0.073 + 0.013, whileH V qq (m 2 Z , µ H = −im Z ) = 1 + 0 + 0 + 0. Hence, the constant terms are moved entirely into the remainder. In this scheme, the resummed result at NNLO+NNLL ϕ is d 2σ res (Y = 0, Q = m Z , µ R = m Z , µ F = m Z /2) = (71.4±0.7 µ ± 0.2 ϕ ) pb/GeV (1.0%). Compared to the last line of table 4, the difference of ±0.8 pb is of the same size as the uncertainties and thus of the typical size we expect for an O(α 3 s ) effect. We now discuss the effect of the resummation on the rapidity spectrum. In fig. 9, we show the remainder R(Y ) normalized to the Born cross section as a function of Y for µ F = m Z (left) and µ F = m Z /2 (right). The behavior discussed for Y = 0 above is similar throughout most of the spectrum. For µ F = m Z , the NNLO corrections are comparably large and of the same size and opposite sign as the NLO contributions, while the NNLO corrections are almost negligible at µ F = m Z /2 over most of the rapidity range. We see again that at µ F = m Z the NNLO corrections involve substantial cancellations between the qq and non-qq channels, which individually are very large. In contrast, at µ F = m Z /2 also the individual corrections to the remainder are very small, again supporting this central choice when including the resummation. (A large part of the rapidity-independent constant shift at NLO will again be canceled by the constant term in H qq .) In fig. 10, we compare the rapidity spectrum at Q = m Z , where the fixed-order calculation at its optimal µ F = m Z is shown for Y < 0 and the resummed result at µ F = m Z /2 for Y > 0.   Figure 10. Rapidity spectrum for pp → Z/γ * → + − at Q = m Z . The fixed-order results are shown for Y < 0 and the resummed results for Y > 0. For the central scale we use µ FO = m Z in both cases, while κ F = 1 (fixed order) and κ F = 1/2 (resummed). The uncertainty bands indicate ∆ µ (fixed order) and ∆ µ ⊕ ∆ ϕ (resummed).
Overall, we find that the NNLO and NNLO+NNLL ϕ predictions provide very similar results. On the one hand, this is reassuring, as it shows that the good convergence of the fixed-order series is not spoiled by the resummation. On the other hand, given the extreme reduction of the perturbative uncertainties in the fixed-order results at the conventional choice of µ F = m Z by a factor of five when going from NLO to NNLO and the substantially larger uncertainties at µ F = m Z /2, one might perhaps be worried that the fixed-order uncertainties are somewhat underestimated, in part due to the accidentally small NNLO contribution. In this respect, the resummed results provide a useful confirmation and increased confidence in the very small perturbative uncertainties in the Drell-Yan predictions.

Conclusion
We have investigated in detail the resummation of timelike logarithms ln 2 (−1) = −π 2 that arise to all orders in perturbation theory and are an important source of perturbative corrections in s-channel color-singlet production processes, which involve a timelike hard momentum transfer. These logarithms can be resummed to all orders using the RG evolution of the corresponding quark or gluon form factors from spacelike to timelike scales.
We have shown how to incorporate the resummed form factor in a completely straightforward manner into predictions for generic inclusive cross sections with arbitrary dependence or cuts on the Born kinematics. We have verified that this does not spoil the perturbative series in all considered cases. We have also discussed the assessment of the uncertainties intrinsic to the resummation.
We first revisited the resummation for the total gluon-fusion cross section, for which it has been discussed before, considering both the production of a generic scalar as well as the SM Higgs boson in the m t → ∞ limit up to N 3 LO+N 3 LL ϕ . For the latter we have also shown how to incorporate quark-mass and electroweak effects. We confirm that the resummation significantly improves the perturbative series, and find that it reduces the perturbative uncertainties at the highest orders by about a factor of two.
For the Higgs rapidity spectrum as well as the cross section with a cut on the Higgs rapidity we obtain results at NNLO+NNLL ϕ , which provide the currently most precise predictions with central values close to what might be expected at N 3 LO, and perturbative uncertainties of ∼ 6%, which are almost a factor of two smaller than at NNLO. Once N 3 LO results for the rapidity dependence become available, we project that the corresponding resummation at N 3 LO+N 3 LL ϕ will provide a similar improvement.
We also studied the resummation of timelike logarithms for quark-induced processes, namely Higgs production through bottom-quark annihilation and the Drell-Yan rapidity spectrum. For the former, the resummation provides a small improvement in the perturbative convergence and resulting uncertainties. For Drell-Yan production, the resummation provides no clear improvement but also no worsening of the predictions, due to the already fast convergence of the fixed-order perturbative series. In this case it provides a useful confirmation of the very small residual perturbative uncertainties.
We conclude that utilizing the resummed timelike quark and gluon form factors is viable and beneficial for obtaining precise and reliable predictions for s-channel color-singlet production processes. the PIER Helmholtz Graduate school. J. M. thanks DESY for hospitality and gratefully acknowledges support by Münster University funds designated for student research.

A Perturbative ingredients
A.1 Master formula for hard Wilson coefficients to three loops The hard matching coefficients C satisfy an RGE of the form which allows us to completely predict the logarithmic structure in terms of the cusp and noncusp anomalous dimension coefficients. We write the perturbative expansion of the hard coefficient as Normalizing C such that the tree-level result is C 0 = 1, the perturbative solution of eq. (A.1) to N 3 LO is given by Here, β n are the beta-function coefficients, Γ n ≡ Γ i n the appropriate quark or gluon cusp anomalous dimensions coefficients, and γ n are the coefficients of the total noncusp anomalous dimension γ in eq. (A.1) as appropriate for the hard coefficient of interest. All required anomalous dimension coefficients are given below in appendix A.2. The results for the nonlogarithmic constant terms C n for the different Wilson coefficients are given below in appendix A. 3 The full expression for the hard function is obtained by squaring C, accounting for cross terms. In the case of H t gg defined in eq. (3.14) the product of C t C gg is reexpanded.

A.2 Anomalous dimensions
We expand the β function of QCD as The coefficients up to four loops in the MS scheme are [134][135][136][137] The coefficients of the MS cusp anomalous dimension to three loops are [138][139][140] Γ q n = C F Γ n , Γ g n = C A Γ n , (for n = 0, 1, 2) , T F n f = 4 3 (4 − π 2 )C A + 5β 0 , The resummation at N 3 LL formally also requires the yet unknown four-loop coefficient Γ i 3 , which we estimate as usual by the Padé approximation and explicitly verify that a variation ±200% only affects the hard evolution kernel U H (and thus the resummed cross section) at the sub-permille level. We therefore neglect this source of theory uncertainty. The gluon noncusp anomalous dimension γ g C enters the RGE for the gluon-to-scalar matching coefficients C gg and C gg in eqs. (3.3) and (3.15). The coefficients in MS up to three loops are [72,76,77] γ g C 0 = −β 0 , The evolution of C t gg in the one-step matching also requires the anomalous dimension γ t of the Wilson coefficient C t arising from integrating out the top quark. It is given by The quark noncusp anomalous dimension γ q C enters the RG eqs. (4.1) and (4.5) for both quark-induced processes we consider. The coefficients in MS up to three loops are [72,77,130,141] The evolution of C S qq also requires the anomalous dimension of the quark Yukawa coupling, which is equivalent to the quark mass anomalous dimension γ m , It is known to five loops [142][143][144][145][146][147][148]. For our main analysis at NNLL we only require the two-loop result, while the three-loop coefficient γ m 2 serves to verify our N 3 LO result for C S qq . The results are

A.3 Constant terms to three loops
In the following, we provide the process-specific nonlogarithmic constant terms C n for the various hard matching coefficients. For C gg , C t , and C V qq , we can collect the results from the literature. The result for C S qq we have extracted from the three-loop scalar quark form factor. By convention, we normalize all coefficients to unity at LO, Note that for all coefficients quoted here, we closely follow the notation from the original publications. For this reason, we set T F = 1/2 in the following results.

A.3.1 Gluon matching coefficient
The finite terms of C gg can be read off from the full result given in ref. [75], available to NNLO, the three-loop coefficient never enters our resummed predictions. For the illustrative values of H V qq given in eq. (4.7) we set N F,V = 0 for the sake of comparison. The explicit three-loop results for C V qq were also extracted in ref. [131] from the threeloop form factor in ref. [74]. We verified that the above results agree with the numerical results for the vector hard function H V qq given in ref. [131] after setting N F,V = 0.
A.3.4 Quark scalar-current matching coefficient As far as we are aware, a result for C S qq has not been given explicitly in the literature so far. The quark scalar form factor F in QCD has been computed to O(α 3 s ) in ref. [103], from which we can extract C S qq . A slight difficulty arises as F is only given at timelike kinematics and fixed µ = Q in ref. [103]. To obtain the full dependence on L = 2 ln(−iQ/µ), we start from the bare form factor F given in ref. [103] and perform its UV-renormalization at an arbitrary MS renormalization point µ. We explicitly checked that the ratio of the timelike to spacelike form factor is IR-finite as required. We then proceed by subtracting the IR poles in F in MS by a multiplicative renormalization factor, In SCET with pure dimensional regularization, the 1/ IR poles in F ( , µ) are the UV poles of the bare Wilson coefficient, so eq. (A.20) is equivalent to the MS renormalization of C S qq . Here we have made explicit that the renormalized quark Yukawa coupling y(µ) is excluded from C S qq . We have also verified that the obtained renormalization factor Z reproduces the correct anomalous dimension for C S qq , i.e. that it satisfies order-by-order in α s , which provides a strong check on the pole structure of F . Equivalently, we also checked that the full result for C S qq (µ) obtained from eq. (A.20) agrees with eq. (A.3) (with Γ cusp ≡ Γ q cusp and γ = 2γ q C − γ m ). For the nonlogarithmic constant terms of C S qq we obtain C S qq 1 = C F (−2 + ζ 2 ) , C S qq 2 = C F C F 6 + 14ζ 2 −

B Fixed-order estimates from resummed timelike logarithms
It is instructive to compare explicitly the fixed-order contributions induced purely by the timelike logarithms in the form factor with the full fixed-order result to assess whether they are indeed a dominant part of the perturbative corrections. However, we also stress that this is not a good way for judging the usefulness of the resummation as a whole, since it does not capture the full resummed result and in particular does not take into account the improvements in perturbative convergence and uncertainties. In ref. [10], such an analysis was carried out for gg → H for the coefficient C δ of the δ(1 − z) term in the partonic cross section. This coefficient is fully determined by the partonic threshold limit z = m 2 H /ŝ → 1, whereŝ is the partonic center-of-mass energy, and factorizes as in eq. (2.6) into the product of the gluon form factor and purely soft contributions. Ref. [10] found that C δ is poorly predicted from the timelike logarithms alone. However, this can be very misleading since the δ(1 − z) coefficient is strongly scheme dependent and not a physical quantity. Rather, this type of analysis should be carried out at the level of the cross section, which is a scheme-independent physical observable. To illustrate this, we repeat the analysis of ref. [10] and compare it to a different convention for the soft function, as well as considering the hadronic K factor.
Following ref. [10], we choose µ FO = m H and work in the pure EFT limit of sec. 3.1 with C X (µ FO ) = 1. (The result for Higgs production is easily restored by reexpanding against |C t | 2 .) The relevant hard function is hence H gg in eq. The result for S δ to O(α 3 s ) can be obtained from ref. [13], which writes the soft function in terms of the standard (plus) distributions Evidently, the resummed results approximate the higher fixed-order terms in the K factor very well, except at NLO, where H (1) and R (1) each contribute about half of the full K factor. This is precisely equivalent to our discussion in sec. 3.1 that the large corrections to the K factor are primarily driven by the timelike logarithms in H, while the nonlogarithmic constant terms, R (n) and H (n) (µ H = −im X ), are much smaller.