Power Corrections for N-Jettiness Subtractions at ${\cal O}(α_s)$
Markus A. Ebert Ian Moult Iain W. Stewart Frank J. Tackmann Gherardo Vita Hua Xing Zhu







MIT-CTP 5007DESY 18-120



1]Markus A. Ebert,
ebert@mit.edu

2,3]Ian Moult,
ianmoult@lbl.gov

1]Iain W. Stewart,
iains@mit.edu

4]Frank J. Tackmann,
frank.tackmann@desy.de

1]Gherardo Vita,
vita@mit.edu

5]and Hua Xing Zhu
zhuhx@zju.edu.cn



[1]Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

[2]Berkeley Center for Theoretical Physics, University of California, Berkeley, CA 94720, USA

[3]Theoretical Physics Group, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA

[4]Theory Group, Deutsches Elektronen-Synchrotron (DESY), D-22607 Hamburg, Germany

[5]Department of Physics, Zhejiang University, Hangzhou, Zhejiang 310027, China






We continue the study of power corrections for N-jettiness
subtractions by analytically computing the complete next-to-leading power corrections at (α_s) for color-singlet production. This includes
all nonlogarithmic terms and all partonic channels for Drell-Yan and gluon-fusion Higgs production.
These terms are important to further improve the numerical performance of the
subtractions, and to better understand the structure of power corrections
beyond their leading logarithms, in particular their universality.  
We emphasize the importance of computing the power corrections
differential in both the invariant mass, Q, and rapidity, Y, of the color-singlet system,
which is necessary to account for the rapidity dependence in the subtractions.
This also clarifies apparent disagreements in the literature. 
Performing a detailed numerical study, we find excellent agreement of our analytic results with a previous numerical extraction.




Power Corrections for N-Jettiness Subtractions at O(α_s)
    [
    
========================================================






§ INTRODUCTION






The precision study of the Standard Model at the LHC, as well as increasingly sophisticated searches for physics beyond the Standard Model, require precision predictions for processes in a complicated hadron collider environment.
When calculating higher-order QCD corrections, the presence of infrared divergences require techniques to isolate and cancel all divergences. Completely analytic calculations are only possible for some of the simplest cases, e.g. <cit.>, while for more complicated processes, in particular those involving jets in the final state, numerical techniques are typically required.

At next-to-leading order (NLO) the FKS <cit.> and CS <cit.> subtraction schemes provide generic subtractions for arbitrary processes, and have been used to great success. At next-to-next-to-leading order (NNLO), due to the more complicated structure of infrared singularities, the development of general subtraction schemes has proven more difficult. While subtraction schemes have been demonstrated both for color-singlet production <cit.>, as well as for several processes involving jets in the final state <cit.>, significant work is still required before efficient NNLO subtractions can be achieved for arbitrary colored final states.

N-jettiness subtractions <cit.> are based on the N-jettiness resolution variable Τ_N <cit.>, and are applicable to generic N-jet final states. They have successfully been applied to NNLO calculations for a variety of color-singlet final states <cit.>, as well as final states involving a single jet <cit.>. They are also a key ingredient in one of the first methods for combining NNLO calculations with parton showers <cit.>. 
The leading-power subtraction terms are given by an all-orders factorization formula derived in Stewart:2009yx, Stewart:2010tn using soft-collinear effective theory (SCET) <cit.>. Required ingredients are explicitly known to NNLO with up to a single jet in the final state  <cit.>.

An important feature of N-jettiness subtractions is that power corrections in the resolution variable can be calculated in an expansion about the soft and collinear limits, allowing the numerical performance of the subtractions to be systematically improved. Recently there has been significant interest in understanding subleading power corrections to collider cross sections <cit.>.
Advances in the understanding of subleading power limits using SCET <cit.> have allowed the leading logarithmic (LL) next-to-leading power (NLP) corrections to be computed at NLO and NNLO <cit.>, with independent calculations of the same terms done by a second group in Boughezal:2016zws,Boughezal:2018mvf. The leading logarithms have also been resummed to all orders for pure glue QCD for 2-jettiness in H→ gg <cit.>.

[2]

The inclusion of the leading logarithms was found to improve the numerical accuracy, and thereby the computational efficiency, of N-jettiness subtractions for color-singlet production by up to an order of magnitude <cit.>. The analytic calculation of the next-to-leading logarithmic (NLL) power corrections is important for several reasons. Theoretically, it greatly furthers our understanding of the power corrections, since the perturbative structure becomes significantly more nontrivial at NLL as compared to at LL. From a practical perspective, they provide further substantial improvements in the numerical performance of the subtractions. In particular, they make the subtractions much more robust in cases where there are accidental cancellations between different channels or where the NLL terms are numerically enhanced relative to the LL terms.

In this paper, we compute the full NLP corrections at (α_s) for both Drell-Yan and gluon-fusion Higgs production in all partonic channels, including the nonlogarithmic terms (which are the NLL terms at (α_s)).
One focus of this paper is to derive a master formula for the NLP corrections to 0-jettiness at (α_s) and to discuss in detail the subleading-power calculation, in particular the treatment of measurements at subleading power, which in our case are the invariant mass Q and rapidity Y of the color-singlet system.
More generally, one can consider measuring any observable that does not vanish for the partonic process at lowest order in perturbation theory,
which we refer to as Born measurements.
Our analysis lays the ground
for extending the calculation of the NLP corrections to higher powers, higher orders in α_s,
and to more complicated processes.

We also perform a detailed numerical study comparing our analytic results with the full nonsingular result extracted numerically from MCFM8 <cit.>, which allows us both to verify our analytic calculation and to probe the typical size of the higher-power corrections in various different partonic channels.
We find that the NLL power corrections can exhibit a much more pronounced Y dependence from PDF effects than is the case at LL, which demonstrates the importance of calculating the power corrections fully differential in the Born phase space.

Our discussion of the Born measurements also allows us to clarify an apparent disagreement in the recent literature regarding the LL power corrections. As we discuss in more detail in compare_lit,
since the calculations in Boughezal:2016zws, Boughezal:2018mvf are not differential in the color-singlet rapidity Y, their results can only be used fully integrated over Y.
[We thank the authors of Boughezal:2016zws, Boughezal:2018mvf for discussions
and confirmation of this point.]
In contrast, the results computed here, which agree with those previously derived by a subset of the present authors in Moult:2016fqy, Moult:2017jsg, are differential in Y.
When integrating over all Y, integration by parts is used to explicitly show that these LL results are equivalent to those of Boughezal:2016zws, Boughezal:2018mvf. In compare_lit we also compare our new differential NLL results with the integrated results of Boughezal:2018mvf.

The outline of this paper is as follows: In Njettiness, we briefly review N-jettiness subtractions for color-singlet production and define our notation. In mom_route, we discuss in some detail the treatment of Born measurements at subleading power. In master_formula, we derive a formula to NLP for the soft and collinear power corrections for 0-jettiness. Although our primary focus is on NLO, the general strategy is valid to higher orders as well. In results_singlet, we use our master formula to derive explicit results for the NLP power corrections at NLO for both Drell-Yan and gluon-fusion Higgs production. In compare_lit, we provide a detailed comparison with the literature for those partonic channels where results are available. In results, we present a detailed numerical study, and compare our analytic results with a previous numerical extraction. We conclude in conclusions.






§ N-JETTINESS SUBTRACTIONS, DEFINITIONS AND NOTATION





In this section we briefly review N-jettiness subtractions <cit.> in the context of color-singlet production, and discuss the structure of the power corrections to the subtraction scheme. This also allows us to define the notation that will be used in the rest of this paper. For a detailed discussion, we refer the reader to Gaunt:2015pea.

To compute a cross section for color-singlet production σ(X), where X denotes some set of cuts on the Born phase space, we write the cross section as an integral over the differential cross section in the resolution variable Τ_0


    σ(X)
       = σ(X, Τ_) + ∫_Τ_Τ_0 σ(X)/Τ_0 ,


where


    σ(X, Τ_)    ≡∫^Τ_Τ_0 σ(X)/Τ_0 .



For a general measure, the 0-jettiness Τ_0 can be defined as <cit.>


    Τ_0 = ∑_i min{2 q_a · k_i/Q_a , 2 q_b · k_i/Q_b} ,


where the sum runs over all hadronic momenta k_i in the final state.
Here, q_a,b are projected Born momenta (referred to as label momenta in SCET), which are given in terms of the
total leptonic invariant mass Q and rapidity Y as


    q_a^μ = x_a n^μ/2 = Q e^Y n^μ/2 ,    
    q_b^μ = x_b ^μ/2 = Q e^-Y^μ/2 ,


where

    n^μ = (1,0,0,1)  ,  ^μ = (1,0,0,-1)

are lightlike vectors along the beam directions.
The choice in Born corresponds to parameterizing the Born phase space in terms of Q and Y,
and this choice already enters in the leading-power factorization theorem, where the beam functions
are evaluated at x_a,b = Q e^± Y.

The Q_a,b measures in Tau0_0 determine the different definitions of 0-jettiness.
Two different definitions, originally introduced in Stewart:2009yx, Berger:2010xi
as beam thrust, are the leptonic and hadronic definitions given by

 
    leptonic:      Q_a    = Q_b = Q  ,        Τ_0^   = ∑_i min{ x_a /Q n ·k_i  ,  x_b /Q ·k_i }

             
     = ∑_i min{ e^Y n ·k_i  ,  e^-Y ·k_i }


    hadronic:        Q_a,b    = x_a,b  ,
     Τ_0^   = ∑_i min{ n ·k_i  ,  ·k_i }
 .

It has been shown <cit.> that the power corrections for the hadronic definition are poorly behaved, and grow exponentially with rapidity, while the e^± Y factor in the measure for the leptonic definition exactly avoids this effect.

For later convenience, we write the dimensionful and dimensionless 0-jettiness resolution variables
in terms of a Y-dependent parameter ρ(Y) as


    Τ_0^x
       = ∑_i min{ρ_x  k_i^+, ρ_x^-1 k_i^- } ,     τ^x ≡Τ_0^x/Q ,


with


    leptonic:         ρ_  = e^Y ,     τ^=Τ_0^/Q , 
    hadronic:         ρ_ = 1 ,     τ^=Τ_0^/Q .


In the following, we will mostly drop the subscript 0 on Τ_0, since there should be no cause for confusion that our results are for 0-jettiness.
For generic results that apply to both the leptonic and hadronic definitions we also drop the superscript and simply use Τ and τ, keeping a generic parameter ρ when necessary.

To implement the N-jettiness subtractions, we now add and subtract a subtraction term to the cross section (suppressing the dependence on the Born measurements X for simplicity)


    σ   =
    σ^(Τ_)
    + ∫_Τ_Τ_0 σ/Τ_0
    + [σ(Τ_) - σ^(Τ_) ]
    
       ≡σ^(Τ_) + ∫_Τ_Τ_0 σ/Τ_0
    + Δσ(Τ_)
     .


Since Τ_0 is a zero jet resolution variable, for τ = Τ_0/Q → 0 we can expand the differential cross section σ/τ and its cumulative σ(τ_) about the soft and collinear limits from τ→ 0 and τ_→ 0 as


    σ/τ   = σ^(0)/τ + σ^(2)/τ+ σ^(4)/τ + … , 
    σ(τ_)
       = σ^(0)(τ_) + σ^(2)(τ_) + σ^(4)(τ_) + … .


Here σ^(0)/τ and σ^(0)(τ_) contain all leading-power terms,


    σ^(0)/τ∼δ(τ)+ [ln^jτ/τ]_+
     ,     σ^(0)(τ_) ∼ln^jτ_ .


These terms must be included in the subtraction term to obtain a finite result, namely


    σ^(Τ_) = σ^(0)(τ_ = Τ_/Q)  [1 + τ_]
     .


The further terms in the series expansion in xsec_expand are suppressed by powers of τ


    τσ^(2k)/τ   ∼τ^k ln^jτ ,     σ^(2k)(τ_) ∼τ_^k ln^jτ_ .


While these terms with k≥ 1 do not need to be included in the subtraction term, the size of the neglected term, Δσ(τ_) is determined by the leading-power corrections that are left out of σ^ sub. Therefore, including additional power corrections in σ^ sub can significantly improve the performance of the subtraction. Indeed, general scaling arguments imply that up to an order of magnitude in performance can be gained for each subleading power logarithm that is included in the subtractions <cit.>. For the leading logarithms, this was explicitly confirmed for most partonic channels in the numerical studies in Moult:2016fqy, Moult:2017jsg. Here, we extend the calculation to the NLL terms at (α_s), which yields the nonlogarithmic terms, hence giving the complete NLP result. The remaining NLO power corrections then scale at worst as α_s τ_^2 logτ_, and will be very small, as we will see in our numerical studies.





§ BORN MEASUREMENTS AT SUBLEADING POWER




We begin by discussing in some detail the treatment of the Born measurements, Q^2 and Y, which plays an important role at subleading powers. We will use the soft and collinear expansions from SCET, which  provide a convenient language when discussing the power expansion of QCD amplitudes at fixed order. We will not need to employ any of the field theory technology from SCET for our analysis here.





 §.§ General Setup and Notation



Consider the production of a color-singlet final state L of fixed invariant mass Q and rapidity Y,
together with an arbitrary measurement Τ
that only acts on hadronic radiation and gives Τ=0 at Born level.
Since the observable Τ resolves soft and collinear emissions
it will induce large logarithms ln(Τ/Q).
Our goal is to expand the cross section in Τ (or τ = Τ/Q)
in order to systematically understand its logarithmic structure.

Consider proton-proton scattering with the underlying partonic process

    a(p_a) + b(p_b) → L(p_1, ⋯) + X(k_1, ⋯)
     ,

where L is the leptonic (color-singlet) final state and X denotes additional QCD radiation.
Its cross section reads

    σ/ Q^2  Y Τ   
     = ∫_0^1ζ_a ζ_b f_a(ζ_a)  f_b(ζ_b)/2 ζ_a ζ_b ^2∫(∏_i ^d k_i/(2π)^d (2π) δ_+(k_i^2) )
       ∫^d q/(2π)^d  |(p_a, p_b; {k_i}, q)|^2
       
         ×
       (2π)^d δ^(d)(p_a + p_b - k - q)  δ(Q^2 - q^2)
        δ(Y - 1/2lnq^-/q^+)
        δ[Τ - Τ̂({k_i})]
     .

Here, the incoming momenta are given by

    p_a^μ = ζ_a n^μ/2 ,    
     p_b^μ = ζ_b ^μ/2 ,

k = ∑_i k_i is the total outgoing hadronic momentum, and q is the total leptonic momentum.
Since our measurements are not sensitive to the details of the leptonic final state,
we have absorbed the leptonic phase space integral into the matrix element,

    |(p_a, p_b; {k_i}, q)|^2    = ∫Φ_L(q)   |(p_a, p_b; {k_i}, {p_j})|^2
     ,
    Φ_L(q)    = ∏_j ^d p_j/(2π)^d (2π) δ_+(p_j^2 - m_j^2)
      (2π)^d δ^(d)(q - ∑_j p_j)
     .

The matrix element  contains the renormalization scale μ^2,
which as always is associated with the renormalized coupling α_s(μ), and may also contain virtual corrections.
For now the measurement function Τ̂({k_i}) is kept arbitrary.

We can now solve the Q^2 and Y measurements to fix the incoming momenta as

    ζ_a(k)    = 1/(k^- +  e^+Y√(Q^2 + k_T^2))
     ,
    ζ_b(k)    = 1/(k^+ +  e^-Y√(Q^2 + k_T^2))
     .

Taking the Jacobian factors from solving the δ functions into account, sigma1 becomes

    σ/ Q^2  Y Τ   
     = ∫(∏_i ^d k_i/(2π)^d (2π) δ_+(k_i^2) )
       f_a(ζ_a)  f_b(ζ_b)/2 ζ_a ζ_b ^4(Q, Y; {k_i}) δ[Τ - Τ̂({k_i})]
     ,

where we defined

    (Q, Y; {k_i}) ≡ |(p_a, p_b, {k_i}, q=p_a+p_b - k)|^2

to stress that the squared matrix element only depends on the Born measurements Q and Y,
which fix the incoming momenta through p_abzeta_ab,
and the emission momenta k_i.
Note that we have left implicit in our notation in sigma2 the dependence of ζ_a,b on k through zeta_ab.
They are restricted to ζ_a,b∈ [0,1], which is implicit in the support of the proton PDFs.





 §.§ Power Expansion in Soft and Collinear Limits



Instead of solving the Τ measurement function to express (some of the) k_i in terms of Τ, we find that a convenient strategy to organize the expansion in Τ is to multipole expand the final state momenta.
At this stage we need only assume that Τ is a observable,
which is true for many definitions of N-jettiness. For such observables, it is known from that
we can organize the cross section in terms of a power counting parameter λ∼√(τ).
All momenta k_i can then be categorized as either collinear or soft modes
(since we work in these are often called ultrasoft, although we will not make this distinction),
whose momenta scale as

    n-collinear:    k_n ∼ Q   (λ^2, 1, λ)
     ,
       -collinear:    k_∼ Q   (1, λ^2, λ)
     ,
       soft:  k_s ∼ Q   (λ^2, λ^2, λ^2)
     ,

where we decomposed each momentum into lightcone coordinates

    k^μ = k^- n^μ/2 + k^+ ^μ/2 + k_⊥^μ≡ (k^+, k^-, k_⊥)
     .

Here n and n̅ are lightlike vectors satisfying n·=2. The components of the momenta that scale like λ^2 are referred to as residual momenta.
The soft momenta are homogeneous, and have purely residual scaling. Overlap between the soft and collinear modes occurring in integrals over final state momenta is removed by the zero-bin subtraction procedure <cit.>.

The benefit of this decomposition is that it allows one to expand sigma2 in λ,
agnostic of the actual measurement Τ.
The LP result is then simply obtained by expanding the cross section through λ^0,
the NLP result by expanding through λ^2, etc.
Note that when performing this expansion, all other factors, such as Q, E_ cm∼λ^0.

While the expansion of the matrix element is of course process dependent,
we can give general expressions for the incoming momentum fractions zeta_ab,
independent of the process and observable Τ.

If k is a soft momentum, then the expansion required at NLP is given by

    ζ_a(k)    = x_a (1 +  k^- e^-Y/Q) + (λ^4)
     ,
    ζ_b(k)    = x_b (1 +  k^+ e^+Y/Q) + (λ^4)
     ,

where we factored out the Born momentum fractions

    x_a = Q e^+Y/ ,    
     x_b = Q e^-Y/ .

In the n-collinear limit, we obtain

    ζ_a(k)    = x_a [ (1 + k^- e^-Y/Q) + k_T^2/2 Q^2] + (λ^4)
     ,
    ζ_b(k)    = x_b [ 1 + (k^+ e^+Y/Q + k_T^2/2 Q^2) ] + (λ^4)
     .

For clarity, we have grouped terms of the same power counting in round brackets.
Similarly, one can obtain the -collinear limit,
or any combination as might appear when combining multiple emissions.





§ MASTER FORMULA FOR POWER CORRECTIONS TO NEXT-TO-LEADING POWER




In this section we derive a master formula for the NLP corrections.
This formula applies to any observable in color-singlet production.
In results_singlet, we will apply it to derive explicit results for Drell-Yan and gluon-fusion Higgs production.





 §.§ General Setup for Color-Singlet SCET I Observables



For reference, we start with the LO cross section for the production
of a color-singlet final state L of invariant mass Q^2 and rapidity Y,
together with an (up to now arbitrary) measurement Τ acting only on hadronic radiation,

    σ^/ Q^2  Y Τ   
     = f_a(x_a)  f_b(x_b)/2 x_a x_b ^4^(Q, Y) δ(Τ)
     ,

where x_a,b = Q/ e^± Y and ^ is the squared matrix element
in the Born kinematics, see Msquared.
For future reference, we also define the LO partonic cross section, σ̂^ LO(Q,Y), by

    σ^/ Q^2  Y Τ   
     = σ̂^(Q,Y)   f_a(x_a)  f_b(x_b)  δ(Τ)
     ,    σ̂^(Q,Y) = ^(Q, Y)/2 x_a x_b ^4 .


Next, consider an additional real emission to the Born process.
Eq. (<ref>) yields

    σ/ Q^2  Y Τ   
     = ∫^d k/(2π)^d  (2π) δ_+(k^2) f_a(ζ_a)  f_b(ζ_b)/2 ζ_a ζ_b ^4 (Q, Y; {k}) δ[Τ - Τ̂({k})]
     ,

where we remind the reader that the incoming momenta p_a,b are given by zeta_ab,

    p_a^μ   = ζ_a(k) n^μ/2
             = (k^- +  e^+Y√(Q^2 + k_T^2)) n^μ/2 ,
    
     p_b^μ   = ζ_b(k) ^μ/2
             = (k^+ +  e^-Y√(Q^2 + k_T^2)) ^μ/2 .

From these solutions, we see the interesting feature that at subleading power, regardless of the type of final-state emission, the momenta entering both PDFs are modified.

Since we do not measure the azimuthal angle of k, it can be integrated over,

    ∫^dk/(2π)^d  (2π) δ_+(k^2)
    = Ω_2-2/4 (2π)^d-1∫_0^∞ k^+  k^-/(k^+ k^-)^
    = (4π)^-2+/Γ(1-)∫_0^∞ k^+  k^-/(k^+ k^-)^ .

Eq. (<ref>) simplifies to

    σ/ Q^2  Y Τ   
     = ∫_0^∞ k^+  k^-/(k^+ k^-)^f_a(ζ_a)  f_b(ζ_b)/(4π)^2 2 ζ_a ζ_b ^4 (4π)^/Γ(1-) (Q, Y; {k}) δ[Τ - Τ̂({k})]
     .

So far, this expression is still exact.
In the next step, we wish to expand the NLO cross section in λ∼Τ/Q.
When Τ is a observable, we can use the EFT knowledge from to expand the momentum k in both collinear and soft limits, as discussed in mom_route.

At NLP, we need to expand sigmaNLO_2 consistently through (λ^2).
The (λ^2) power corrections arise from the following sources:

 
  * The incoming momenta ζ_a,b: While collinear and soft limits yield
 quite different power expansions, both give a well-defined expansion in λ.
 We thus simply define the expansion
 
    ζ_a,b = x_a,b[ 1/z_a,b + Δ_a,b^(2) + (λ^4) ]
      ,

 where z_a,b∼λ^0 and Δ_a,b^(2)∼λ^2.
 We have pulled out the Born momentum fractions x_a,b and written 1/z_a,b as a fraction
 for later convenience.
 Explicit expressions can be obtained from zeta_ab_softzeta_ab_coll,
 and will be given below.

 
  * PDFs: Since the momenta ζ_a,b enter the PDFs, these also have to be power expanded,
 
    f_a,b(ζ_a,b)
         = f_a,b(x_a,b/z_a,b) + x_a,bΔ_a,b^(2)  f'_a,b(x_a,b/z_a,b) + (λ^4)
      
       ≡ f_a,b + x_a,bΔ_a,b^(2)  f'_a,b + (λ^4)
      .


 
  * Flux factor: Similar to the PDFs, we have to expand the flux factor
 
    1/ζ_a ζ_b = z_a z_b/x_a x_b[ 1 - z_a Δ^(2) - z_b Δ^(2)_b ]
      .


 
  * Matrix element: The expansion of the matrix element depends both on the process and the considered limit.
 Here, we define the LP and NLP expansions by
 
    (Q, Y; {k}) = ^(0)(Q,Y; {k})) + ^(2)(Q,Y; {k}) + ⋯ .

 In the soft limit ^(0)∼λ^-4 and ^(2)∼λ^-2,
 while in the collinear limit ^(0)∼λ^-2 and ^(2)∼λ^0. In both cases we have the scaling ∫ k^+  k^- A^(2j)∼λ^2j, which is why the soft and collinear corrections enter at the same order. 


 For example, for a 2 → 2 process, the matrix element can be written in terms
 of the Mandelstam variables
 
    s_ab   = 2 p_a · p_b = Q^2/z_az_b[ 1 + z_a Δ^(2)_a + z_b Δ^(2)_b  + (λ^4) ]
       ,
    
      s_ak   = -2 p_a · k = -k^+ Q e^+Y(1/z_a +Δ^(2)_a  + (λ^4) )
       ,
    
      s_bk   = -2 p_b · k = -k^- Q e^-Y(1/z_b +Δ^(2)_b  + (λ^4) )
      .

 Since these terms now have a definite power counting, one can simply insert mandelstamgeneral
 into (Q, Y; {k}) and expand to the required order in λ.

 
  * Measurement: Depending on the observable, the measurement function Τ̂ may also
 receive power corrections. Since 0-jettiness is defined in terms of n, n̅, Q and Y, none of which receive power corrections in our approach, we do not have such corrections, and will therefore not write them explicitly in the following formulae.  More generally, these could be obtained from expanding δ[Τ - Τ̂({k})].

Inserting all these expansions into sigmaNLO_2 and expanding consistently to (λ^2), we obtain the LP result

    σ^(0)/ Q^2  Y Τ   
     = ∫_0^∞ k^+  k^-/(k^+ k^-)^z_a z_b   f_a f_b/(4π)^2 2 x_a x_b ^4(4π)^/Γ(1-)^(0)(Q, Y; {k})  δ[Τ - Τ̂({k})]

and the NLP master formula

    σ^(2)/ Q^2  Y Τ   
     = ∫_0^∞ k^+  k^-/(k^+ k^-)^z_a z_b/(4π)^2 2 x_a x_b ^4(4π)^/Γ(1-)δ[Τ - Τ̂({k})]
       
         ×{^(0)(Q, Y; {k}) [
        f_a f_b (-z_a Δ_a^(2) - z_b Δ_b^(2))
       + x_a Δ_a^(2) f'_a f_b + x_b Δ_b^(2) f_a f'_b ]
       
           
       + f_a f_b  ^(2)(Q, Y; {k})
       } ,

where z_a,b≡ z_a,b(k) and Δ_a,b^(2)≡Δ_a,b^(2)(k) are defined by zeta_ab_expansion.
Note that the LP limits of the matrix elements are universal,
and hence sigmaLP holds independently of the process,
i.e. it only depends on the observable Τ.
Although the focus of this paper is on the power corrections,
in app_LP we provide a brief derivation of the leading-power terms.
Likewise, the A^(0) term together with the square bracketed factor on the second line of sigmaNLP is universal, such that all the process dependence arises from the last A^(2) term.
We will discuss this in more detail in universality.

In the following, we will evaluate sigmaNLP in both the soft and collinear limit
for 0-jettiness, Tau0, whose measurement function for one emission is given by

    δ[Τ - Τ̂(k)]
     = Θ(ρ k^+ - ρ^-1 k^-) δ(Τ - ρ^-1 k^-)
     + Θ(ρ^-1 k^- - ρ k^+) δ(Τ - ρ k^+)
     .

The value of ρ depends on the specific definition of Τ, as given in Njet_def.





 §.§ Collinear Master Formula for 0-Jettiness




The expansion of the incoming momenta ζ_a,b for an n-collinear emission
k ∼ (λ^2, 1, λ) has been given in zeta_ab_coll,

    ζ_a(k)    = x_a [ (1 + k^- e^-Y/Q) + k_T^2/2 Q^2] + (λ^4)
     ,
    ζ_b(k)    = x_b [ 1 + (k^+ e^+Y/Q + k_T^2/2 Q^2) ] + (λ^4)
     ,

so the explicit expressions for the expansion zeta_ab_expansion are
 
     z_a       = (1 + k^- e^-Y/Q)^-1
 ,         Δ^(2)_a = k_T^2/2 Q^2
 ,

     z_b       = 1
 ,         Δ^(2)_b = (k^+ e^+Y/Q + k_T^2/2 Q^2 )
 .

Since an n-collinear emission satisfies k^- ≫ k^+,
the 0-jettiness measurement delta_Tau simplifies to

    δ[Τ - Τ̂(k)] = δ(Τ - ρ k^+)
     .

Note that the integration in sigmaNLP also includes the region k^- → 0,
where the assumption k^- ≫ k^+ is invalid.
Indeed, this region corresponds to the soft expansion.
It is guaranteed by the zero-bin subtraction procedure that
this overlap regime between the soft and collinear limits
is not double counted <cit.>.
An important benefit of 0-jettiness and our setup here is that the
zero-bin contribution that removes the overlap is scaleless and vanishes in
pure dimensional regularization, such that we do not need to consider it further.

Eq. <ref> fixes the k^+ integral in sigmaNLP.
It is also useful to write the remaining k^- integration in terms of z_a using auxiliary_coll, giving

    k^- = Q e^Y 1 - z_a/z_a ,    ∫_0^∞ k^- = ∫_x_a^1  z_a/z_a^2  Q e^Y
     .

Here the lower bound on the integration follows from the physical support of the PDF
f_a(x_a/z_a).
Plugging back into sigmaNLP, we obtain the n-collinear master formula

    σ_n^(2)/ Q^2  Y Τ   
     = ∫_x_a^1  z_a/z_a 1/2 x_a x_b ^4Q e^Y/ρ(Q Τ e^Y/ρ)^-z_a^/(1-z_a)^(4π)^/(4π)^2 Γ(1-)
         ×{Τ e^Y/Q ρ^(0)(Q, Y, {k}) [
        f_a  f_b (1-z_a)^2 - 2/2 z_a
       + 1-z_a/2z_a  x_a  f'_a  f_b + 1+z_a/2z_a  f_a  x_b  f'_b ]
       
           
       + f_a  f_b  ^(2)(Q, Y, {k})
       } ,

where k is given by

    k^μ = Q e^Y 1-z_a/z_an^μ/2 + Τ/ρ^μ/2
           + √(Q Τe^Y/ρ1-z_a/z_a) n_⊥^μ .

It only remains to plug in the expansions of the matrix element A^(0) and A^(2) and to expand in .
Note that even in the n-collinear case considered here, both the n and -collinear incoming momenta
receive power corrections, see p_ab_2, leading to derivatives of both PDFs in sigma_NLO_NLP_coll_FULL.

The analogous results for the -collinear limit are obtained in the same manner, giving

    σ_^(2)/ Q^2  Y Τ   
    = ∫_x_b^1  z_b/z_b 1/2 x_a x_b ^4Q ρ/e^Y(Q Τρ/e^Y)^-z_b^/(1-z_b)^(4π)^/(4π)^2 Γ(1-)
         ×{Τρ/Q e^Y^(0)(Q, Y, {k}) [
       f_a f_b (1-z_b)^2 - 2/2 z_b
       + 1+z_b/2z_b  x_a f'_a  f_b
       + 1-z_b/2z_b  f_a  x_b f'_b ]
       
           
       + f_a  f_b  ^(2)(Q, Y, {k})
       } ,

where k is given by

    k^μ = Τρ n^μ/2 + Q/e^Y1-z_b/z_b^μ/2
           + √(Q Τρ/e^Y1-z_b/z_b) n_⊥^μ .







 §.§ Soft Master Formula for 0-Jettiness




The expansion of the incoming momenta ζ_a,b for a soft emission
k ∼ (λ^2, λ^2, λ^2) has been given in zeta_ab_soft,

    ζ_a(k)    = x_a (1 +  k^- e^-Y/Q) + (λ^4)
     ,
    ζ_b(k)    = x_b (1 +  k^+ e^+Y/Q) + (λ^4)
     ,

so the explicit expressions for the expansion zeta_ab_expansion are
 
     z_a       = 1
 ,           Δ^(2)_a = k^- e^-Y/Q
 ,

     z_b       = 1
 ,           Δ^(2)_b = k^+ e^+Y/Q
 .
Plugging back into sigmaNLP, we get

    σ_s^(2)/ Q^2  Y Τ   = ∫_0^∞ k^+  k^-/(k^+ k^-)^1/(4π)^2 2 x_a x_b ^4(4π)^/Γ(1-)δ[Τ - Τ̂({k})]
    {1/Q^(0)(Q, Y, {k})
    
           ×[
       f_a  f_b  (-k^- e^-Y - k^+ e^+Y)
       + k^- e^-Y  x_a f'_a  f_b + k^+ e^+Y  f_a  x_b f'_b ]
    
         
       + f_a  f_b  ^(2)(Q, Y, {k})
    } .

Here, the measurement is given by delta_Tau,

    δ[Τ - Τ̂(k)]
     = Θ(ρ k^+ - ρ^-1 k^-) δ(Τ - ρ^-1 k^-)
     + Θ(ρ^-1 k^- - ρ k^+) δ(Τ - ρ k^+)
     .


We can further simplify softmaster by utilizing the fact that
^(0) and ^(2) have a well defined dependence on k^+ or k^-
because of power counting and mass dimension,

    ^(0)(Q, Y, {k}) = ^(0)(Q,Y)/k^+ k^- ,  ^(2)(Q, Y, {k}) = _+^(2)(Q,Y)/k^+ + _-^(2)(Q,Y)/k^- .

Here the 's are process-dependent expressions that depend on the Born measurements Q and Y,
but are independent of both k^+ and k^-.
This implies that the k^± integrals in softmaster have the generic structure

    ∫_0^∞ k^+  k^-/(k^+ k^-)^δ[Τ - Τ̂(k)]/(k^+)^α (k^-)^β
     = ρ^α-βΤ^1-α-β-2(1/+α-1 + 1/+β-1)
     .

We then find the soft NLP master formula

    σ_s^(2)/ Q^2  Y Τ   = 1/(4π)^2 2 x_a x_b ^4(4π)^/Γ(1-)1/Τ^-2/Q1-2/1- {^(0)(Q,Y)
    
           ×[ f_a(x_a)  f_b(x_b) (- ρ/e^Y -e^Y/ρ)
       + ρ/e^Y  x_a f'_a(x_a)  f_b(x_b) + e^Y/ρ  f_a(x_a)  x_b f'_b(x_b) ]
    
         
       +  f_a(x_a)  f_b(x_b)  [ ρ Q  _+^(2)(Q,Y) + Q/ρ _-^(2)(Q,Y) ]
    } .






 §.§ Universality of Power Corrections for 0-Jettiness




Having derived our master formulas, in this section we comment on the universality of the power corrections.
In both the collinear and soft limits the power corrections arising from the derivatives of the PDFs
and from the expansion of the flux factor are proportional to the LP matrix element ^(0)(Q,Y),
see sigma_NLO_NLP_coll_FULLsigma_NLO_NLP_soft.
Since the factorization properties of ^(0)(Q,Y)
are universal, most of the NLP corrections are universal as well,
in the sense that they essentially reduce to the LO cross section times a universal factor,
as we will make explicit below.
The only process-dependent piece arises from the NLP expansion ^(2)(Q,Y) of the matrix element.
We stress that this limit is defined in our particular choice of Born measurements Q^2 and Y.
Using different observables, e.g. q^±, the NLP corrections to ζ_a,b in zeta_ab would change,
inducing also a change of the NLP matrix element.





  §.§.§ Universality of Collinear Limit




We begin by considering the n-collinear limit of a real emission amplitude in detail.
We consider the Born process

    κ_a(q_a) + κ_b(q_b) → L(q_a + q_b)
     ,

where κ_i denotes all quantum numbers, including flavor,
of the incoming partons, and L is the leptonic final state of momentum q = q_a + q_b.
The incoming momenta for the hard collision are given by

    q_a^μ = x_a n^μ/2 = Q e^+Yn^μ/2 ,    
     q_b^μ = x_b ^μ/2 = Q e^-Y^μ/2 .

Now consider that parton a arises from an n-collinear splitting of a parton with flavor a',

    κ'_a(q'_a) + κ_b(q_b) → L(q'_a + q_b - k) + κ_1(k)
     .

To describe this at leading power, we only need the (λ^0) relations
for the momenta of the incoming partons, which can be read off from zeta_ab_coll2auxiliary_coll,

    q'_a^μ   = q_a^μ/z_a + (λ^2) = Q e^Y/z_an^μ/2 + (λ^2)
     .

The n-collinear emission is given by k_collinear,

    k^μ = Q e^Y 1-z_a/z_an^μ/2 + Τ/ρ^μ/2
           + √(Q Τe^Y/ρ1-z_a/z_a) n_⊥^μ .

It follows that the leptonic momentum q' = q'_a + q'_b - k = q + (λ^2)
is equal to the Born momentum q = q_a + q_b, and hence the collinear splitting does not
affect the leptonic phase space at LP.

The LP limit only exists if the splitting κ'_a →κ_a + κ'_1 is allowed,
in which case it is given by the (λ^-2) piece of the squared amplitude,

    _a' b → L k(Q, Y, {k})
     = 8π^2/Q e^Y k^+ P_a a'(z_a, ) _a b → L^(Q, Y) + (λ^0)
     ,

where the 1/k^+ gives rise to the λ^-2 behavior of the amplitude.
Here the P_aa' are the -dimensional splitting functions which are summarized in app_LP.

Recall that in our case, the measurement fixes k^+ = Τ/ρ.
In the notation of expansion_M2, we hence have for the LP matrix element

    _a' b → L k^(0)(Q, Y, {k})
     = 8π^2ρ/Q Τ e^Y P_a a'(z_a, ) _a b → L^(Q, Y)
     .

These results enable us to explicitly give the universal part of the NLP result in the collinear limit.
Inserting into the collinear master formula sigma_NLO_NLP_coll_FULL
and converting to the MS scheme, we find

    σ_n^(2)/ Q^2  Y Τ
       =
     σ̂^(Q,Y) /4 πe^Y/Q ρ×(Q Τ/μ^2e^Y/ρ)^-e^γ_E/Γ(1-)∫_x_a^1  z_a/z_a z_a^/(1-z_a)^ P_aa'(z_a,)
       
         ×[ f_a'(x_a/z_a) f_b(x_b) (1-z_a)^2 - 2/z_a
       + 1-z_a/z_a x_a f'_a'(x_a/z_a) f_b(x_b)
       + 1+z_a/z_a f_a'(x_a/z_a) x_b f'_b(x_b) ]
    
         
     + ∫_x_a^1  z_a/z_a f_a(x_a/z_a)  f_b(x_b)/2 x_a x_b ^4Q e^Y/ρ(Q Τ e^Y/ρ)^-z_a^/(1-z_a)^(4π)^/(4π)^2 Γ(1-) ^(2)(Q, Y, {k})
     .

Here, we factored out the LO partonic cross section sigmaLO_2,
which is only possible because the collinear splitting leaves the leptonic
momentum invariant at LP.
We have made explicit the universal piece and nonuniversal components.
As already discussed, the full nonuniversal structure arises from the NLP matrix element
^(2)(Q, Y, {k}) in the last line.

It would also be interesting to understand if there is a universal structure to ^(2)(Q,Y).
This has recently been studied in  Nandan:2016ohb for pure n gluon scattering amplitudes at the level of the Cachazo-He-Yuan scattering equations <cit.>, where it was proven that in the subleading power collinear limits the tree-level amplitude factorizes into a convolution of the n-1 gluon integrand and a universal collinear kernel. It would be interesting to understand this at the level of the amplitude itself, as well as for fermions. Unlike at leading power, we do not expect that there are universal subleading power splitting functions that are simply functions of z, but there may exist splitting functions that involve differential or integral operators, as occurs in the soft limit at subleading power <cit.>. Understanding this will be particularly important for generalizing the calculation of the power corrections to more complicated processes.





  §.§.§ Universality of Soft Limit




As for the collinear case, the LP soft limit of the matrix element is universal.
Following similar steps as in univ_coll, one can express the LP soft limit by

    _ab → L k^(0)(Q, Y; {k})
     = 16 π^2𝐂/k^+ k^-×_a b → L^(Q, Y)
     ,

which only exists for ab=gg, qq̅ and where 𝐂 = C_A, C_F is the appropriate Casimir constant.
We thus obtain

    σ_s^(2)/ Q^2  Y Τ   
     = σ̂^/Q𝐂/π[ 1/ - lnΤ^2/μ^2 - 1 ]
       [ f_a(x_a) f_b(x_b) (- ρ/e^Y -e^Y/ρ)
       
       
       + ρ/e^Y  x_a  f'_a(x_a)  f_b(x_b) + e^Y/ρ  f_a(x_a)  x_b f'_b(x_b) ]
    
         
       + f_a(x_a)  f_b(x_b)/(4π)^2 2 x_a x_b ^4(4π)^/Γ(1-)1/Τ^-21-2/1-[ ρ |_+^(2)(Q,Y)|^2 + 1/ρ |_-^(2)(Q,Y)|^2 ]
     .

As for the collinear case, this emphasizes that the terms arising
from the expansion of the PDFs and flux factor are universal,
in the sense that they only depend on the universal LP soft limit of the amplitude.
The only nonuniversal contributions are |_±^(2)|^2.
However, these terms can in fact be derived from universal formulae <cit.> involving differential operators.
This has been recently studied in the threshold limit,
where one only requires soft contributions <cit.>.
However, when one is away from the threshold limit as considered here, one in general requires collinear contributions,
which as discussed above, are not (yet) known to be universal.







§ POWER CORRECTIONS AT NLO FOR COLOR SINGLET PRODUCTION




In this section we give explicit results for the full NLP correction for 0-jettiness at NLO
for Higgs and Drell-Yan production in all partonic channels.

Since we only consider cases that are s-channel processes at Born level, the LO matrix element only depends on Q and one can factor out the LO partonic cross section σ̂^(Q).
We write the NLP cross section as

    σ^(2,n)/ Q^2  Y Τ   = σ̂^(Q)  (/4π)^n
        ∫_x_a^1  z_a/z_a∫_x_b^1  z_b/z_b[
    f_i(x_a/z_a) f_j(x_b/z_b) C_f_i f_j^(2,n)(z_a, z_b, Τ)
    
         
    + x_a/z_a f'_i(x_a/z_a) f_j(x_b/z_b) C_f_i' f_j^(2,n)(z_a, z_b, Τ)
    + f_i(x_a/z_a) x_b/z_b f'_j(x_b/z_b) C_f_i f_j'^(2,n)(z_a, z_b, Τ)
    ]
     ,

where as always

    x_a = Q e^Y/ ,   x_b = Q e^-Y/ .

We will always express the real emission amplitudes in terms of the Mandelstam variables

    s_ab = 2 p_a · p_b
     ,  
     s_ak = -2 p_a · k
     ,  
     s_bk = -2 p_b · k
     .

This allows us to straightforwardly obtain the LP and NLP expansion using mandelstamgeneral.
We will give an explicit example of the derivation of the soft and collinear master formulas
for the gg → H g channel, and only summarize the results in the other channels.





 §.§ Gluon-Fusion Higgs Production




We begin by considering Higgs production in gluon fusion in the m_t →∞ limit.
At NLP, there are three different partonic channels, gg→ Hg, qq̅→ H g and q g → H q,
which we consider separately. The calculation for gg → Hg is shown in full detail
as an illustration of our master formulae.
The LL power corrections were computed in <cit.>. Ref. <cit.> also computed the qq̅→ Hg NLL power corrections. 
The NLL power corrections for all partonic channels for gluon fusion Higgs were computed in <cit.>.
We will compare with these results in compare_lit.

Throughout this section we consider on-shell Higgs production,
for which the partonic cross section is given by

    σ̂^(Q,Y)
     = ^(Q,Y)/2 x_a x_b ^4
     = 2πδ(Q^2 - m_H^2) |_gg→ H^(Q)|^2/2 Q^2 ^2 .

The LO matrix element in d=4-2ϵ dimensions is given by
<cit.>

    |^_gg→ H(Q)|^2    = ^2 Q^4/576 π^2 v^2(4π^2/m_t^2)^2Γ^2(1+)/1- .





  §.§.§ gg → Hg




The spin- and color-averaged squared amplitude for g(p_a) + g(p_b) → H(q) + g(k) is given by <cit.>

    _gg → H g(Q, Y, {k})    
     = ^_gg→ H(Q) ×8 π C_A ^2/ Q^4 (1-)
         ×[ (1-2) Q^8 + s_ab^4 + s_ak^4 + s_bk^4/s_ab s_ak s_bk
       + /2(Q^4 + s_ab^2 + s_ak^2 + s_bk^2)^2/s_ab s_ak s_bk]
     .




  
n-Collinear Limit
Expanding M2_ggHg using mandelstamgeneralauxiliary_coll,
the LP and NLP limits of the matrix element are obtained as

    ^(0)_gg → Hg(Q, Y; {k})    
     = 16 πα_s C_A ^2 ϵ^_gg→ H(Q)  (1-z_a+z_a^2)^2 /(1-z_a) z_aρ e^-Y/Q Τ ,
    ^(2)_gg → Hg(Q, Y; {k})    
     = 16 πα_s  C_A ^2 ϵ^_gg→ H(Q)  1/Q^2 z_a^2[ 1 + 5 z_a^2 - z_a^3 + 2 z_a^4 - z_a^5 - 2 z_a^2 1/1-]
     .

Since our scaling variable is λ∼√(Τ/Q),
we clearly see that ^(0)∼λ^-2 and ^(2)∼λ^0,
as required at LP and NLP.

Inserting these expansions into sigma_NLO_NLP_coll_FULL and converting to the MS scheme yields 

    σ_n^(2)/ Q^2  Y Τ   
     = σ̂^_gg→ H(Q) × C_A/π∫_x_a^1  z_a/z_a1/Qe^Y/ρ( Τ Q/μ^2e^Y/ρ)^-z_a^/(1-z_a)^e^γ_E/Γ(1-)
         ×{ (1-z_a+z_a^2)^2 /(1-z_a) z_a[ f_a   f_b (1 - z_a)^2 - 2/2 z_a
       + x_a f'_a   f_b  1-z_a/2 z_a
       + f_a   x_b f'_b  1+z_a/2 z_a]
       
       
       + f_a   f_b 1/z_a^2[ 1 + 5 z_a^2 - z_a^3 + 2 z_a^4 - z_a^5
      - 2 z_a^2 1/1-]
      } .

To expand this in , we collect all powers of (1-z_a) and then use the distributional identity

    (1-z_a)^-1-ϵ = -δ(1-z_a)/ϵ + _0(1-z_a) + (ϵ)
     ,

where _0(1-z)=1/(1-z)_+ is the usual plus distribution.
We also combine the two separate f_a f_b pieces, as at this level there is no use to further
distinguish the universal and process dependent pieces. This yields

    σ_n^(2)/ Q^2  Y Τ   
     = σ̂^_gg→ H(Q) ×/4 π 
       4 C_A e^Y/Qρ∫_x_a^1  z_a/z_a
         ×{
         f_g(x_a/z_a)   f_g(x_b)
         [ (1/ - lnQ Τ e^Y/μ^2 ρ) δ(1-z_a)
         
       
          + 1 - 2 z_a + 8 z_a^2 - 14 z_a^3 + 12 z_a^4 - 10 z_a^5 + 3 z_a^6/2z_a^2_0(1-z_a)
          ]
       
           
       + x_a/z_a f'_g(x_a/z_a)   f_g(x_b)  (1-z_a+z_a^2)^2/2 z_a
           
       + f_g(x_a/z_a)   x_b f'_g(x_b)  [
       (-1/+ lnQ Τ e^Y/μ^2 ρ) δ(1-z_a)
       
       
       + (1 + z_a)(1 - z_a + z_a^2)^2/2 z_a^2 _0(1-z_a) ]
      } .

Comparing to sigma_NLP, we can read off the n-collinear kernels,

    C_f_g f_g,n^(2,1)(z_a, z_b, Τ)
       = 4C_A e^Y/Qρ[
    (1/ - lnQ Τ e^Y/μ^2 ρ) δ(1-z_a)
    
         
    + 1 - 2 z_a + 8 z_a^2 - 14 z_a^3 + 12 z_a^4 - 10 z_a^5 + 3 z_a^6/2z_a^2 _0(1-z_a) ]  δ(1-z_b)
     ,
    
    C_f'_g f_g,n^(2,1)(z_a, z_b, Τ)
       = 4C_A e^Y/Qρ (1-z_a + z_a^2)^2/2 z_a δ(1-z_b)
     ,
    
    C_f_g f'_g,n^(2,1)(z_a, z_b, Τ)
       = 4C_A e^Y/Qρ [
    ( -1/ + lnQ Τ e^Y/μ^2 ρ) δ(1-z_a)
    
         
    + (1+z_a)(1-z_a+z_a^2)^2/2z_a^2 _0(1-z_a) ] δ(1-z_b)
     .




  
Soft Limit
To expand the matrix element in the soft limit, we use mandelstamgeneralauxiliary_soft to obtain

    _gg→ H g(Q, Y, {k}) = ^_gg→ H(Q) ×
      16 πα_s C_A ^2 ϵ/k^+ k^- + (λ^0)
     .

Note that the first term scales as (k^+ k^-)^-1∼λ^-4,
while there is no (λ^-2) component.
The NLP term in the expansion of the amplitude thus vanishes, and in the notation of altMdef we have

    _gg → H g^(0)(Q) = ^_gg→ H(Q) × 16 πα_s C_A ^2 ϵ ,    _gg → H g^(2)(Q) = 0
     .

Inserting into sigma_NLO_NLP_soft and converting to the MS scheme yields

    σ_s^(2)/ Q^2  Y Τ   = σ̂_gg→ H^(Q) α_s/4π× 4C_A 1/Q( 1/ - lnΤ^2/μ^2 - 1 )
    
         ×[ f_g(x_a) f_g(x_b) (- ρ/e^Y -e^Y/ρ)
       + ρ/e^Y  x_a f'_g(x_a)  f_g(x_b)
       + e^Y/ρ  f_g(x_a)  x_b f'_g(x_b) ]
     .

Since there is no NLP matrix element, one can also obtain this from the
universal expression for the soft limit in sigma_NLO_NLP_soft_UNIV.
Comparing to sigma_NLP, we can read off the soft kernel,

    C_f_g f_g,s^(2,1)(z_a, z_b, Τ)
       = 4C_A 1/Q(-e^Y/ρ - ρ/e^Y) ( 1/ - lnΤ^2/μ^2 - 1 ) δ(1-z_a) δ(1-z_b)  ,
    
    C_f'_g f_g,s^(2,1)(z_a, z_b, Τ)
       = 4C_A 1/Qρ/e^Y ( 1/ - lnΤ^2/μ^2 - 1 ) δ(1-z_a) δ(1-z_b)  ,
    
    C_f_g f'_g,s^(2,1)(z_a, z_b, Τ)
       = 4C_A 1/Qe^Y/ρ ( 1/ - lnΤ^2/μ^2 - 1 ) δ(1-z_a) δ(1-z_b)
     .




  
Final Result
Adding the n-collinear kernel kernel_Cgg_n, the -collinear kernel which follows from symmetry,
and the soft kernel kernel_Cgg_s, all poles in  cancel as expected, and we obtain

    C_f_g f_g^(2,1)(z_a, z_b, Τ)
       = 4 C_A e^Y/Qρ[ (lnΤρ/Q e^Y + 1 ) δ(1-z_a)
    
         
       + 1 - 2 z_a + 8 z_a^2 - 14 z_a^3 + 12 z_a^4 - 10 z_a^5 + 3 z_a^6/2z_a^2 _0(1-z_a)
       ] δ(1-z_b)
    
         
       + ( a ↔ b  , ρ/e^Y→e^Y/ρ)
     ,
    
    C_f'_g f_g^(2,1)(z_a, z_b, Τ)
       = 4 C_A ρ/Qe^Y δ(1-z_a) [
       ( - lnΤ e^Y/Q ρ - 1  ) δ(1-z_b)
       
       
       + (1+z_b)(1-z_b+z_b^2)^2/2z_b^2 _0(1-z_b) ]
    
         
    + 4C_A e^Y/Qρ(1-z_a + z_a^2)^2/2 z_aδ(1-z_b)
     ,
    
    C_f_g f_g'^(2,1)(z_a, z_b, Τ)
       = 4 C_A e^Y/Q ρ[
       ( - lnΤρ/Q e^Y - 1  ) δ(1-z_a)
       
       
       + (1+z_a)(1-z_a+z_a^2)^2/2z_a^2 _0(1-z_a) ] δ(1-z_b)
    
         
    + 4C_A ρ/Q e^Yδ(1-z_a) (1-z_b + z_b^2)^2/2 z_b .

Substituting these results into sigma_NLP yields the NLP cross section for gg→ Hg at NLO.





  §.§.§ gq → Hq



The gq→ Hq channel has power corrections at both LL and NLL.
The spin- and color-averaged squared amplitude for g(p_a) + q(p_b) → H(q) + q(k) is given by <cit.>

    _gq → Hq(Q, Y, {k})
     = - ^_gg→ H(Q) × 8 π C_F ^21/Q^4 s_bk[ s_ab^2 + s_ak^2 -  (s_ab+s_ak)^2 ]
     .




  
Soft Limit
The LP soft limit vanishes, since a leading-power soft interaction (which is eikonal) cannot change a n-collinear quark into a n-collinear gluon and soft quark. However this does occur at NLP in the soft expansion and yields

    ^(2)_gq → Hq(Q, Y, {k})    
     = ^_gg→ H(Q) × 8 π C_F ^21-/Q k^- e^-Y ,

and the soft kernel is given by

    C_f_g f_q , s^(1,2)(z_a, z_b, Τ)
    = 2 C_F e^Y/Qρ ( 1/ - lnΤ^2/μ^2 - 2 )
     δ(1-z_a) δ(1-z_b)
     .




  
-Collinear Limit
The -collinear limit has both a LP and NLP contribution, given by

    ^(0)_gq → Hq(Q, Y, {k})    
     = ^_gg→ H(Q)  × 8 π C_F e^Y/ρ1 + (1-z_b)^2 -  z_b^2/Q Τ z_b , 
    ^(2)_gq → Hq(Q, Y, {k})    
     = ^_gg→ H(Q) × 4 π C_F ^24 - z_b^3 + z_b^4 -  z_b^2 (4 - z_b + z_b^2)/Q^2 z_b^2 .

The -collinear kernel is obtained as

    C_f_g f_q , ^(1,2)(z_a, z_b, Τ)
       = C_F ρ/Qe^Y δ(1-z_a) 2 - 2 z_b + 5 z_b^2 - 5 z_b^3 + 2 z_b^4/z_b^2 ,
    
    C_f'_g f_q , ^(1,2)(z_a, z_b, Τ)
       = C_F ρ/Qe^Y δ(1-z_a) (1+z_b)[1 + (1-z_b)^2]/z_b^2 ,
    
    C_f_g f_q' , ^(1,2)(z_a, z_b, Τ)
       = C_F ρ/Qe^Y δ(1-z_a) (1-z_b)[1 + (1-z_b)^2]/z_b .





  
n-Collinear Limit
The n-collinear limit vanishes at LP. The NLP expansion of the matrix element gives

    ^(2)_gq → Hq(Q, Y, {k})    
     = ^_gg→ H(Q)  ×8 π C_F ^2 (1-)/Q^2 (1-z_a) z_a .

The only nonvanishing kernel is

    C_f_g f_q , n^(1,2)(z_a, z_b, Τ)
       = 2C_F e^Y/Qρ[
    ( -1/ + lnQ Τ e^Y/μ^2 ρ + 1 ) δ(1-z_a)
    + _0(1-z_a)/z_a] δ(1-z_b)
     .




  
Final Result
Combining the n-collinear, -collinear, and soft kernels, the 1/ pole vanishes,
and we obtain the final results,

    C_f_g f_q^(1,2)(z_a, z_b, Τ)
       = 2C_F e^Y/Qρ[
       ( -lnΤρ/Q e^Y - 1) δ(1-z_a)
       + _0(1-z_a)/z_a] δ(1-z_b)
    
         
    + C_F ρ/Qe^Y δ(1-z_a) 2 - 2 z_b + 5 z_b^2 - 5 z_b^3 + 2 z_b^4/z_b^2 ,
    
    C_f'_g f_q^(1,2)(z_a, z_b, Τ)
       = C_F ρ/Qe^Y δ(1-z_a) (1+z_b)[1 + (1-z_b)^2]/z_b^2 ,
    
    C_f_g f_q'^(1,2)(z_a, z_b, Τ)
       = C_F ρ/Qe^Y δ(1-z_a) (1-z_b)[1 + (1-z_b)^2]/z_b .

Substituting these results into sigma_NLP yields the NLP cross section for gq→ Hq at NLO.





  §.§.§ qg → Hq



The final results needed in sigma_NLP for qg → H q follow from kernel_Cgq
by flipping a ↔ b, e^Y/ρ↔ρ/e^Y and f_g ↔ f_q,

    C_f_q f_g^(1,2)(z_a, z_b, Τ)
       = 2C_F ρ/Qe^Y δ(1-z_a) [
    ( - lnΤ e^Y/Q ρ - 1 ) δ(1-z_b) + _0(1-z_b)/z_b]
    
         
    + C_F e^Y/Qρ 2 - 2 z_a + 5 z_a^2 - 5 z_a^3 + 2 z_a^4/z_a^2 δ(1-z_b)
     ,
    
    C_f'_q f_g^(1,2)(z_a, z_b, Τ)
       = C_F e^Y/Qρ (1-z_a)[1 + (1-z_a)^2]/z_a δ(1-z_b)
     ,
    
    C_f_q f_g'^(1,2)(z_a, z_b, Τ)
       = C_F e^Y/Qρ (1+z_a)[1 + (1-z_a)^2]/z_a^2 δ(1-z_b)
     .






  §.§.§ qq̅→ Hg



The qq̅→ Hg channel first contributes at NLL.
It was first given in <cit.> and then in <cit.>, which agreed,
but we reproduce it here for completeness.
The squared matrix element, including the average on the initial state spin and colors, is  given by <cit.>

    _qq̅→ Hg(Q, Y, {k})    
     = ^_gg→ H(Q) ×64π/3 C_F ^21-/Q^4 s_ab[ s_ak^2 + s_bk^2 -  (s_ak+s_bk)^2 ]
     ,

With our choice of Born measurements, the soft limit vanishes both at LP and NLP, leaving only the collinear NLP correction.
The LP collinear limit also vanishes, leaving only the NLP n-collinear limit

    ^(2)_qq̅→ Hg(Q, Y, {k}) = ^_gg→ H(Q) ×64π/3 C_F ^2 (1-)^2 (1-z_a)^2/Q^2 z_a ,

and the -collinear result is obtained by replacing z_a ↔ z_b.
Combining both, we obtain the kernel for sigma_NLP

    C_f_q f_q̅^(2,1)(z_a, z_b, Τ)
    = 16 C_F/3 1/Q[ e^Y/ρ (1-z_a)^2/z_a δ(1-z_b)
       + ρ/e^Y δ(1-z_a) (1-z_b)^2/z_b]
     .






 §.§ Drell-Yan Production



We now consider the Drell-Yan process p p → Z/γ^* → l^+ l^- , and for brevity denote it as p p → V.
At NLO we have the partonic channels qq̅→ Vg and qg → Vq.
The LL power corrections for these channels were calculated to NNLO in <cit.>.

For Drell-Yan, it is important to be able to include off-shell effects.
The LO partonic cross section as a function of the leptonic invariant mass Q is given by

    σ̂^(Q) = 4πα_em^2/3 N_c Q^2 ^2[ Q_q^2 + (v_q^2 + a_q^2)(v_l^2 + a_l^2) - 2 Q_q v_q v_l (1 - m_Z^2/Q^2)/(1-m_Z^2/Q^2)^2
             + m_Z^2 Γ_Z^2 / Q^4]
     .

Here, v_l,q and a_l,q are the standard vector and axial couplings of the leptons and quarks to the Z boson, and we have integrated over the l^+ l^- phase space.




  §.§.§ qq̅→ Vg



We first consider the partonic channel qq̅→ Vg. The squared amplitude is given by <cit.>

    |_qq̅→ Vg|^2    
     = |_qq̅→ V|^2 ×8 π C_F ^2/Q^2[ (1-ϵ) (s_ak/s_bk+s_bk/s_ak)+2 s_ab Q^2/s_ak s_bk -2ϵ]
     .




  
Soft Limit
With our setup, the soft limit of the matrix element has no NLP correction,

    _qq̅→ Vg(Q,Y,{k})
     = ^_q q̅→ V(Q) ×16 π C_F ^2/k^+ k^- + (λ^0)
     ,

and the soft kernels are given by

    C_f_q f_q̅, s^(2,1)(z_a, z_b, Τ)    
     = 4 C_F (-e^Y/Q ρ - ρ/Q e^Y)
       ( 1/ - lnΤ^2/μ^2 - 1 )
       δ(1-z_a) δ(1-z_b)
     ,
    
      C_f'_q f_q̅, s^(2,1)(z_a, z_b, Τ)    
     = 4 C_F ρ/Q e^Y ( 1/ - lnΤ^2/μ^2 - 1 )
       δ(1-z_a) δ(1-z_b)
     ,
    
      C_f_q f'_q̅, s^(2,1)(z_a, z_b, Τ)    
     = 4 C_F e^Y/Q ρ ( 1/ - lnΤ^2/μ^2 - 1 )
       δ(1-z_a) δ(1-z_b)
     .




  
Collinear Limit
The n-collinear expansion of the matrix element yields
(at NLP, we only need →0)

    ^(0)_qq̅→ Vg(Q,Y,{k})
        = ^_q q̅→ V(Q) × 8 π C_F ^2ρ/e^Y1 + z_a^2 - (1-z_a)^2/Q Τ (1-z_a) ,
    ^(2)_qq̅→ Vg(Q,Y,{k})
        = ^_q q̅→ V(Q) × 4 π C_F 1 - z_a + z_a^2 - z_a^3/Q^2 z_a .

The n-collinear kernel is

    C_f_q f_q̅, n^(2,1)(z_a, z_b, Τ)    
      = 4 C_F e^Y/Q ρ [
        (1/ - lnQ Τ e^Y/μ^2 ρ)
        δ(1-z_a)
      
       
        + 1/2(z_a - 2)(1+z_a^2) _0(1-z_a) ] δ(1-z_b)
     ,
    
     C_f'_q f_q̅, n^(2,1)(z_a, z_b, Τ)    
      = 4 C_F e^Y/Q ρ 1+z_a^2/4δ(1-z_b)
     ,
    
     C_f_q f'_q̅, n^(2,1)(z_a, z_b, Τ)    
      =  4 C_F e^Y/Q ρ [
        ( - 1/ + lnQ Τ e^Y/μ^2 ρ)
        δ(1-z_a)
      
       
       + (1+z_a)(1+z_a^2)/4z_a _0(1-z_a) ] δ(1-z_b)
     .




  
Final Result
Adding the n,  and s kernel, we get

    C_f_q f_q̅^(2,1)(z_a, z_b, Τ)    
      = 4 C_F e^Y/Q ρ [ (lnΤρ/Q e^Y + 1 ) δ(1-z_a)
        + 1/2 (z_a - 2)(1+z_a^2) _0(1-z_a) ] δ(1-z_b)
       
           + ( e^Y/ρ→ρ/e^Y , a ↔ b )
     ,
    
     C_f'_q f_q̅^(2,1)(z_a, z_b, Τ)    
      = 4 C_F ρ/Q e^Y δ(1-z_a)[ ( -lnΤ e^Y/Q ρ - 1) δ(1-z_b)
      + (1+z_b)(1+z_b^2)/4z_b _0(1-z_b) ]
      
         
      + 4 C_F e^Y/Q ρ 1+z_a^2/4 δ(1-z_b)
     ,
    
     C_f_q f'_q̅^(2,1)(z_a, z_b, Τ)    
     = 4 C_F e^Y/Q ρ[ ( -lnΤρ/Q e^Y - 1) δ(1-z_a)
      + (1+z_a)(1+z_a^2)/4z_a _0(1-z_a) ] δ(1-z_b)
      
         
      + 4 C_F ρ/Q e^Y δ(1-z_a) 1+z_b^2/4 .

Substituting these results into sigma_NLP yields the NLP cross section for qq̅→ Vg at NLO.





  §.§.§ qg→ Vq



Next we consider the partonic channel qg→ Vq. The squared amplitude is given by <cit.>

    _qg → Vq(Q,Y,{k})
     = - ^_q q̅→ V(Q) ×8 π T_F ^2/Q^2 (1-)[ (1-ϵ) (s_ab/s_bk+s_bk/s_ab)+2 s_ak Q^2/s_ab s_bk -2ϵ]
     .




  
Soft Limit
The LP soft limit vanishes, and the NLP soft expansion is given by

    ^(2)_qg → Vq(Q,Y,{k})
     =  ^_q q̅→ V(Q) × 8 π T_F ^2e^Y/Q k^- .

The soft kernel is given by

    C_f_q f_g,s^(2,1) =  2 T_F e^Y/Q ρ( 1/ - lnΤ^2/μ^2 - 1 ) δ(1-z_a) δ(1-z_b)
     .




  
n-Collinear Limit
The n-collinear limit does not contribute at LP, since the LP interaction
can not change the -collinear gluon into a -collinear antiquark.
The NLP matrix element is given by

    ^(2)_qg → Vq(Q,Y,{k})    
     =  ^_q q̅→ V(Q)  × 8 π T_F
     1 + (1-z_a)^2 - z_a^2/(1-) (1-z_a)Q^2 ,

and the collinear kernel is

    C_f_q f_g,n^(2,1)(z_a, z_b, Τ)    
     = 2 T_F e^Y/Q ρ [ (-1/ + lnQ Τ e^Y/μ^2 ρ) δ(1-z_a)
       
       
       + [1 + (1-z_a)^2] _0(1-z_a) ] δ(1-z_b)
     .





  
-Collinear Limit
The -collinear limit is IR finite, so we work in d=4,

    ^(0)_qg → Vq(Q,Y,{k})    
     =  ^_q q̅→ V(Q)  × 8π T_F e^Y/ρ1 - 2 z_b + 2 z_b^2/Q Τ ,
    ^(2)_qg → Vq(Q,Y,{k})    
     =  ^_q q̅→ V(Q)  × 4π T_F 1+z_b+4z_b^2-8z_b^3+4z_b^4/Q^2 z_b .

The -collinear kernel is given by

    C_f_q f_g,^(2,1)(z_a, z_b, Τ)    = T_F ρ/Q e^Y δ(1-z_a)   (1-z_b)(1 + 8 z_b - 6 z_b^2)
     ,
    
     C_f'_q f_g,^(2,1)(z_a, z_b, Τ)    = T_F ρ/Q e^Y δ(1-z_a)  (1+z_b)(1 - 2 z_b + 2 z_b^2)/z_b ,
    
     C_f_q f_g',^(2,1)(z_a, z_b, Τ)    = T_F ρ/Q e^Y δ(1-z_a)   (1 - z_b)(1-2 z_b + 2 z_b^2)
     .




  
Final Result
Adding the s,n, kernels, the pole in  cancels and we get

    C_f_q f_g^(2,1)(z_a, z_b, Τ)    =
      2 T_F e^Y/Q ρ [
      ( -lnΤρ/Q e^Y - 1 ) δ(1-z_a) + [1 + (1-z_a)^2] _0(1-z_a) ] δ(1-z_b)
      
         
      + T_F ρ/Q e^Y δ(1-z_a)   (1-z_b)(1 + 8 z_b - 6 z_b^2)
     ,
    
     C_f'_q f_g^(2,1)(z_a, z_b, Τ)
        = T_F ρ/Q e^Y δ(1-z_a)  (1+z_b)(1 - 2 z_b + 2 z_b^2)/z_b ,
    
     C_f_q f_g'^(2,1)(z_a, z_b, Τ)
        = T_F ρ/Q e^Y δ(1-z_a)   (1 - z_b)(1-2 z_b + 2 z_b^2)
     .

Substituting these results into sigma_NLP yields the NLP cross section for qg→ Vq at NLO.





  §.§.§ g q→ Vq



For completeness, we also give the explicit results for the gq→ Vq channel,
which can easily be obtained from kernel_Cqg_DY by flipping a ↔ b,
e^Y/ρ↔ρ/e^Y and f_q ↔ f_g,

    C_f_g f_q^(2,1)(z_a, z_b, Τ)    =
      2 T_F ρ/Q e^Y δ(1-z_a) [
      ( -lnΤ e^Y/Q ρ - 1 ) δ(1-z_b) + [1 + (1-z_b)^2] _0(1-z_b) ]
      
         
      + T_F e^Y/Q ρ  (1-z_a)(1 + 8 z_a - 6 z_a^2) δ(1-z_b)
     ,
    
     C_f'_g f_q^(2,1)(z_a, z_b, Τ)
        = T_F e^Y/Q ρ  (1 - z_a)(1-2 z_a + 2 z_a^2) δ(1-z_b)
     ,
    
     C_f_g f'_q^(2,1)(z_a, z_b, Τ)
        = T_F e^Y/Q ρ (1+z_a)(1 - 2 z_a + 2 z_a^2)/z_a δ(1-z_b)
     .






§ COMPARISON WITH INTEGRATED RESULTS IN THE LITERATURE




In this section, we compare our NLO results to previous results in the literature.

The LL results presented by a subset of the present authors in Moult:2016fqy,Moult:2017jsg fully agree with the results obtained in this paper.

The results in Boughezal:2016zws,Boughezal:2018mvf are given only integrated over the color-singlet rapidity Y, and hence take quite a different form at the integrand level.
To compare to them, we integrate our results over Y,
which allows us to use integration by parts to bring our results into the same integrated
form as those in Boughezal:2016zws,Boughezal:2018mvf.
For leptonic Τ, whose definition involves Y, we find that Boughezal:2018mvf uses a different definition, and hence we cannot make a meaningful comparison.
For hadronic Τ, whose definition is independent of Y, we find explicit agreement
for the LL results after integrating over Y.

At NLL, the results obtained here for the power corrections differential in Y, for both the leptonic and hadronic definitions and all partonic channels, are new. After integrating over Y we find almost complete
agreement with the hadronic results of Boughezal:2018mvf, up to a relatively simple
term.
[This missing term has been confirmed by the authors of Boughezal:2018mvf
and was corrected in their version 2.]

Since there are a number of differences in our treatment compared to Boughezal:2016zws,Boughezal:2018mvf, we provide a detailed comparison in this section.
In phasespace_NLO we discuss our different treatments of the NLO phase space and of the Born measurements, and show that the rapidity dependence cannot be easily reconstructed from the results in Boughezal:2016zws,Boughezal:2018mvf.
In compare_results_lit we provide an explicit comparison of the results for the gg → Hg channel integrated over rapidity at LL and NLL, both analytically and numerically.





 §.§ Treatment of the NLO Phase Space




The derivation in Boughezal:2018mvf differs from ours here (and that in Moult:2016fqy,Moult:2017jsg) in that it is not differential in the rapidity Y.
To explore the differences arising from this, we give a brief derivation of the NLO phase space following the same steps as Boughezal:2018mvf.
Note that in the following we always work with an on-shell process,
in contrast to our more general setup in master_formula.
We also only consider the case k^+ < k^-, since the case k^+ > k^- follows by symmetry.

We start with the expression for the NLO phase space as given in Boughezal:2018mvf,

    PS_NLO/Τ   = Τ^-ϵ (4π^2)^-ϵ/8πΓ(1-ϵ)∫ξ_a ξ_b   f_g(ξ_a) f_g(ξ_b)/2ξ_a ξ_b^2(  Q_a ξ_a/x_a)^1-ϵ
         ×∫ z_a (1-z_a)^-ϵδ( ξ_a ξ_b z_a ^2  -m_H^2 -Q_a ξ_a/x_aΤ)
     ,

where s = ^2, Q_a is defined in Tau0_2 and x_a arises from the Τ measurement.

We can derive a similar expression in our notation, including in addition the rapidity measurement
as done in our main derivation. Denoting the incoming momenta at NLO by q'_a,b, we have from sigma1

    PS^/ Y Τ   
     = ∫_0^1 ξ_a ξ_b f_a(ξ_a) f_b(ξ_b)/2 ξ_a ξ_b ^2^2∫^d k/(2π)^d (2π) δ_+(k^2)
       ∫^d q/(2π)^d (2π) δ_+(q^2 - Q^2)
       
         ×  (2π)^d δ(q'_a + q'_b - q - k)
       δ(Y - 1/2lnq^-/q^+)
       δ[Τ - Τ̂(k)]
    
       
     = 1/8π(4π^2)^/Γ(1-)∫_0^1 ξ_a ξ_b f_a(ξ_a) f_b(ξ_b)/2 ξ_a ξ_b ^2∫_0^∞ k^+  k^-/(k^+ k^-)^δ[Τ - Τ̂(k)]
       
         ×δ(ξ_a ξ_b ^2 - ξ_a  k^+ - ξ_b  k^- - Q^2)  δ(Y - 1/2lnξ_a  - k^-/ξ_b  - k^+)
     .

As in PS_NLO_FP, we assume that k^+ < k^- to set Τ̂(k) = ρ k^+, which gives

    PS^/ Y Τ   
     = 1/8 π(4π^2)^/Γ(1-)∫_0^1 ξ_a ξ_b f_a(ξ_a) f_b(ξ_b)/2 ξ_a ξ_b ^2∫_0^∞ k^-/ρ( ρ/Τ k^-)^
         ×δ(ξ_a ξ_b ^2 - Q^2 - ξ_a Τ/ρ - ξ_b  k^-)  δ(Y - 1/2lnξ_a  - k^-/ξ_b  - Τ/ρ)
     .

Following Boughezal:2018mvf, we now change variables via k^- = ξ_a  (1-z_a),

    PS^/ Y Τ   
     = Τ^-/8π(4π^2)^/Γ(1-)∫_0^1 ξ_a ξ_b f_a(ξ_a) f_b(ξ_b)/2 ξ_a ξ_b ^2( ξ_a /ρ)^1-∫ z_a   (1-z_a)^-
         ×δ(z_a ξ_a ξ_b ^2 - ξ_a Τ/ρ - Q^2)  δ(Y - 1/2lnz_a ξ_a/ξ_b - Τ/ρ)
     .

Up to the rapidity measurement from the final δ function, we find complete agreement with PS_NLO_FP
if we identify

    ρ≡ρ(Y) = x_a /Q_a(x_a) .


At this step, our treatment differs from the one in Boughezal:2018mvf.
Since we explicitly implement measurement δ functions for both Q and Y, we can uniquely solve
for ξ_a and ξ_b in terms of Q and Y or equivalently x_a and x_b,

    ξ_a    =  e^+Y/2 z_a^2 [  Τ e^Y/ρ (1-z_a) + √((Τ e^Y/ρ)^2 (1-z_a)^2 + 4 Q^2 z_a^2)]
     ,
    ξ_b    =  e^-Y/2 z_a [  Τ e^Y/ρ (1+z_a) + √((Τ e^Y/ρ)^2  (1-z_a)^2 + 4 Q^2 z_a^2)]
     .


This holds for both ρ = 1 and ρ = e^Y.
This is equivalent to zeta_ab (where we used the notation ζ_a,b instead of ξ_a,b here).
The reason this expression looks different is just because in zeta_ab we performed this step before fixing k^+ in terms of
Τ and before changing variables from k^- to z_a via k^- = ξ_a  (1-z_a).

Following a similar strategy as in master_formula, one can now replace ξ_a,b in PS_NLO_1 by the solution xi_ab,
take the Jacobian from solving the δ functions into account, and then simply expand in Τ.
The main difference to the derivation in master_formula is that here, one directly expands the phase space in Τ, while in master_formula we expanded in terms of the generic power-counting parameter λ.

In Boughezal:2018mvf, there is only the Q^2 measurement but no rapidity measurement, i.e. Y is implicitly integrated over.
Hence, there is only one constraint for the two variables ξ_a,ξ_b, whose solution is not unique.
They choose to perform the variable transformation from ξ_a,b to new variables x̃_a,b defined by

    ξ_a = x̃_a^2 x̃_b ^2/z_a x̃_a x̃_b ^2 - Q_a(x̃_a) Τ ,     ξ_b    = x̃_b
     .


We write x̃_a,b here to distinguish these from the Born variables x_a,b = Q e^± Y/ that appear in the Born-projected momenta in Born.
While they satisfy x̃_a x̃_b = Q^2/^2 due to the Q^2 measurement constraint, (1/2)ln(x̃_a / x̃_b) is not equal to the rapidity Y, which would require the solution
in xi_ab.

In Boughezal:2018mvf, the x̃_a,b defined by xi_a
also enter in the definition of the 0-jettiness measure in Tau0_0 in place of x_a,b.
As a result, the nonhadronic Τ definition in Boughezal:2018mvf is not the same as the usual leptonic Τ with ρ = e^Y that we use.
Their hadronic definition is the same as ours, as it has no x_a,b dependence.
Therefore in the following we restrict our comparison to the hadronic definition.

We also note that one cannot easily recover
the rapidity dependence from the integrands of the final results in Boughezal:2018mvf.
To see this explicitly, consider inserting the rapidity measurement by comparing PS_NLO_FPPS_NLO_1, which gives

    δ(Y - 1/2lnz_a ξ_a/ξ_b - Τ/ρ)
        = δ(Y - 1/2lnz_a ξ_a^(0)/ξ_b^(0))
     
         
       + Τ/2δ'(Y - 1/2lnz_a ξ_a^(0)/ξ_b^(0))
         ( ξ'_b^(0)/ξ_b^(0) - ξ'_a^(0)/ξ_a^(0)
                - 1/ξ_b^(0)ρ)
     + (Τ^2)
     .

On the right-hand side we have carried out the power expansion about Τ→ 0 and the superscript ^(0) denotes the results for these variables at LP, while ξ_a^'(0) = dξ_a/dΤ|_Τ→ 0, etc. This accounts for the fact that in general the ξ_a,b can depend on Τ themselves.

Equation (<ref>) shows that one cannot use the LP expression
δ[Y - (1/2)ln(z_a ξ_a^(0)/ξ_b^(0))] = δ[Y - (1/2)ln(x̃_a/x̃_b)]
to recover the Y dependence from the x̃_a,b dependence of the results in Boughezal:2018mvf, as this does not account for the additional power corrections
induced by the Y measurement in the second line of Y_insert.
This implies that the results in Boughezal:2018mvf and also those in Boughezal:2016zws
cannot be used when being differential in rapidity or integrated over bins of rapidity, but only integrated over all Y. This was also confirmed to us by the authors.





 §.§ Explicit Comparison to Results in the Literature for gg→ Hg




Our final results take a quite different form than those in Boughezal:2016zws,Boughezal:2018mvf.
For us, both ξ_a and ξ_b receive power corrections resulting in derivatives for both PDFs.
In contrast, the variable transformation in xi_a for the case of k^+ < k^-
does not yield power corrections for ξ_b and hence no derivatives of f_b,
while the expansion of ξ_a yields derivatives of f_a (and vice versa for k^+ > k^-).
Due to this different form, one cannot directly compare the integrands of the two results,
but one needs to use integration by parts to bring the results into the same form,
as we will now show explicitly. In particular, we will show that the results of Moult:2016fqy,Moult:2017jsg, obtained also here, do agree with the results of Boughezal:2016zws, Boughezal:2018mvf at LL when integrating over all Y.

Integrating our result over Y, and transforming the integration variables to x_a,b = Q e^± Y/, we obtain from sigma_NLP

    σ^(2,1)/Τ   
     = /4π∫_0^1  x_a  x_b  
       2πδ(x_a x_b ^2 - m_H^2) |_gg→ H^(m_H)|^2/2 x_a x_b ^2∫_x_a^1  z_a/z_a∫_x_b^1  z_b/z_b 
         ×[
        f_i(x_a/z_a) f_j(x_b/z_b) C_f_i f_j^(2,1)(z_a, z_b, Τ)
       + x_a/z_a f'_i(x_a/z_a) f_j(x_b/z_b) C_f_i' f_j^(2,1)(z_a, z_b, Τ)
       
       
       + x_b/z_b f_i(x_a/z_a) f'_j(x_b/z_b) C_f_i f_j'^(2,1)(z_a, z_b, Τ)
       ]
     .


We will show the integration by parts explicitly for the f_i f'_j piece.
Let us denote the piece we wish to integrate by parts by D^(2,1),
which can be chosen freely.
To integrate over Y_1 < Y < Y_2, we switch the integration variables x_a, x_b back to Q^2 and Y, use that

    x_b/z_b f'_j(x_b/z_b)
     = Q e^-Y/ z_b f'_j(Q e^-Y/ z_b)
     = - / Y f_j(Q e^-Y/ z_b)
     ,

and integrate by parts with respect to Y.
Combining the resulting pieces with those in sigma_NLP_NLO_1, we find

    σ^(2,1)/Τ   
     = /4π∫_0^1  x_a  x_b  
       2πδ(x_a x_b ^2 - m_H^2) |_gg→ H^(m_H)|^2/2 x_a x_b ^2∫_x_a^1  z_a/z_a∫_x_b^1  z_b/z_b 
         ×{
        f_i(x_a/z_a) f_j(x_b/z_b)
        [ C_f_i f_j^(2,1)(z_a, z_b, Τ) + / Y D^(2,1)(z_a, z_b, Τ) ]
       
       
       + x_a/z_a f'_i(x_a/z_a) f_j(x_b/z_b)
        [ C_f_i' f_j^(2,1)(z_a, z_b, Τ) + D^(2,1)(z_a, z_b, Τ) ]
       
       
       + x_b/z_b f_i(x_a/z_a) f'_j(x_b/z_b)
         [ C_f_i f_j'^(2,1)(z_a, z_b, Τ) - D^(2,1)(z_a, z_b, Τ) ]
       }
       
    - /4π∫ Q^2  
       2πδ(Q^2 - m_H^2) |_gg→ H^(m_H)|^2/2 Q^2^2∫_Q e^Y/^1  z_a/z_a∫_Q e^-Y/^1  z_b/z_b
           ×
       f_i(Q e^Y/ z_a)  f_j(Q e^-Y/ z_b)
       D^(2,1)(z_a, z_b, Τ) |_Y=Y_1^Y=Y_2 .

The dependence on D^(2,1) exactly cancels in this expression.
We can choose D^(2,1) freely to obtain different forms of the Y-integrated  result.
The last term in sigma_NLP_NLO_2 is the boundary contribution, which vanishes as Y_1,2→±∞, i.e. only if one is fully inclusive in Y.
They do in general contribute when placing acceptance cuts on Y.

We now work out explicitly the required integration by parts both at LL and NLL
to bring our results into the integrated form as given in Boughezal:2016zws,Boughezal:2018mvf.
For concreteness, we focus on the gg→ Hg channel. For the reasons mentioned
earlier, we can only compare the results for the hadronic Τ definition.





  §.§.§ Comparison at LL



At LL, our results in kernel_Cgg simplify to

    C_f_g f_g^(2,1),LL(z_a, z_b, Τ)    
     = 4 C_A [
       e^Y/Q ρlnΤρ/Q e^Y +
       ρ/Q e^YlnΤ e^Y/Q ρ]
       δ(1-z_a) δ(1-z_b)
     ,
    
     C_f'_g f_g^(2,1),LL(z_a, z_b, Τ)    
     =  - 4 C_A ρ/Q e^YlnΤ e^Y/Q ρδ(1-z_a) δ(1-z_b)
     ,
    
     C_f_g f_g'^(2,1),LL(z_a, z_b, Τ)    
     = - 4 C_A e^Y/Q ρlnΤρ/Q e^Yδ(1-z_a) δ(1-z_b)
     .

These agree with the earlier results obtained by a subset of the current authors in Moult:2016fqy,Moult:2017jsg.
Note that for strict LL accuracy, one can also write the logarithms as
ln(Τ/Q) ±ln (ρ/e^Y) and only keep the ln(Τ/Q) at LL, while
including the ±ln(ρ/e^Y) pieces in the NLL contributions.
(This is the convention used in Moult:2016fqy, Moult:2017jsg and
in results.)
Here, we keep them as part of the LL result, as they are relevant for the comparison with Boughezal:2018mvf.

Up to a trivial change in notation, the LL result given in Boughezal:2018mvf
for hadronic Τ is

    σ^NLP<cit.>_LL/Τ^   
     =  C_A/π∫_0^1 x̃_a x̃_b   2πδ(x̃_a x̃_b ^2 - m_H^2) |^_gg→ H(m_H)|^2/2 x̃_a x̃_b ^2
         ×[
        - x̃_a f'_g(x̃_a) f_g(x̃_b)  x̃_a /m_H^2lnΤ^/x̃_a 
        - x̃_b f_g(x̃_a) f'_g(x̃_b)  x̃_b /m_H^2lnΤ^/x̃_b ]
     .

As discussed before, the x̃_a,b here are not equal to the Born variables x_a,b.

Inserting our LL result in kernel_Cgg_LL with ρ = 1 into sigma_NLP_NLO_1, we have

    σ^(2,1)_LL/Τ^   
     =  C_A/π∫_0^1  x_a  x_b   2πδ(x_a x_b ^2 - m_H^2) |^_gg→ H(m_H)|^2/2 x_a x_b ^2
         ×[
        f_g(x_a) f_g(x_b) ( e^Y/m_HlnΤ^ e^-Y/m_H + e^-Y/m_HlnΤ^ e^Y/m_H)
       
           
       - x_a f'_g(x_a) f_g(x_b) e^-Y/m_HlnΤ^ e^Y/m_H
       - x_b f_g(x_a) f'_g(x_b) e^Y/m_HlnΤ^ e^-Y/m_H]
     ,

where e^Y = √(x_a/x_b).
At the integrand level, the two results clearly have a different form, as was also remarked in Moult:2017jsg,Boughezal:2018mvf.

To show explicitly that sigma_ggHg_LL_had_litsigma_ggHg_had_LL_1 do agree,
we integrate by parts to move the f_g f_g contribution in sigma_ggHg_had_LL_1
into the f_g f'_g and f'_g f_g terms. Using sigma_NLP_NLO_2, we can achieve this by choosing

    D^(2,1)(z_a, z_b, Τ^)
       = 4 C_A ( - e^Y/m_HlnΤ^ e^-Y/m_H
                     + e^-Y/m_HlnΤ^ e^Y/m_H)
       δ(1-z_a) δ(1-z_b)
     .

Integrating over Y_1 < Y < Y_2 and using that e^Y = √(x_a / x_b) and m_H = √(x_a x_b), sigma_ggHg_had_LL_1 becomes

    σ^(2,1)_LL(Y_1, Y_2)/Τ^   
     =  C_A/π∫_0^1  x_a  x_b
         2πδ(x_a x_b ^2 - m_H^2) |_gg→ H^(m_H)|^2/2 x_a x_b ^2
           ×[
       - x_a f'_g(x_a) f_g(x_b) x_a /m_H^2lnΤ^/x_a 
       - x_b f_g(x_a) f'_g(x_b) x_b /m_H^2lnΤ^/x_b ]
    
         
    +  C_A/π2π|_gg→ H^(m_H)|^2/2 m_H^2 ^2
       f_g(m_H e^Y/) f_g(m_H e^-Y/)
    
           ×[ e^Y/m_HlnΤ^ e^-Y/m_H -
         e^-Y/m_HlnΤ^ e^Y/m_H]
       |_Y_1^Y_2
         
     +  C_A/π∫_0^1 x_a  x_b
         2πδ(x_a x_b ^2 - m_H^2) |_gg→ H^(m_H)|^2/2 x_a x_b ^2
       f_g(x_a) f_g(x_b) e^Y + e^-Y/m_H .

The first two lines exactly reproduce sigma_ggHg_LL_had_lit.
The following two lines are the boundary term from integration by parts,
which vanishes as Y_1,2→±∞.
The last line is a NLL effect and can be neglected for the LL comparison.
(It is induced by the integration by parts acting on the
Y dependence kept inside the argument of the logarithms.)
Therefore, the two expressions in sigma_ggHg_LL_had_litsigma_ggHg_had_LL_1
agree at LL and at integrated level if and only if one integrates over all rapidity.



To illustrate this numerically, the Y-integrated results are compared in comparison_cumulant_Taucm.
First note that the hadronic LL results in sigma_ggHg_had_LL_1 do not exactly correspond to those previously given in Moult:2016fqy,Moult:2017jsg.
This is due to the formally NLL terms proportional to ln(ρ/e^Y), discussed below kernel_Cgg_LL, which are dropped in the strict LL results in Moult:2016fqy,Moult:2017jsg,
but are kept in sigma_ggHg_had_LL_1. The analogous NLL terms proportional to ln(x̃_a,b)
are also kept in Boughezal:2016zws,Boughezal:2018mvf and sigma_ggHg_LL_had_lit.
Dropping these NLL terms in sigma_ggHg_LL_had_litsigma_ggHg_had_LL_1, our and their LL results defined in terms of the same ln(Τ/m_H) agree exactly,
as shown by the solid red and blue dashed curves in comparison_cumulant_Taucm.
[
In the first version of Boughezal:2018mvf an analogous numerical comparison showed a disagreement
between their integrated LL results and our corresponding result from Moult:2017jsg. This was only due to an incorrect comparison. We thank the authors of Boughezal:2018mvf for confirming this.]
The long-dashed orange and dotted blue curves in comparison_cumulant_Taucm show the
results when using instead ln(Τ e^± Y/m_H) or ln[Τ/(x̃_a,b)] to multiply
the LL coefficients. The observed difference to the solid red/dashed blue strict LL result has the size of a typical NLL contribution. There is also a very small difference between the long-dashed orange and dotted blue results due to the fact that e^± Y/m_H is not exactly the same as x̃_a,b. This difference is exactly accounted for by the last line in sigma_ggHg_had_LL.






  §.§.§ Comparison at NLL




We now extend our comparison of the Y-integrated results to NLL,
focusing again only on the g g → H g channel,
which contains all possible complications.
The full NLL result of Boughezal:2018mvf can be written as

    σ^NLP<cit.>/Τ^   
     =  C_A/π∫_0^1 x̃_a x̃_b   2πδ(x̃_a x̃_b ^2 - m_H^2)
       |^_gg→ H(m_H)|^2/2 x̃_a x̃_b ^2∫_x̃_a^1  z_a/z_ax̃_a /m_H^2
         ×{
       f_g(x̃_a/z_a) f_g(x̃_b) [
        ((1-z_a+z_a^2)^2/z_a^2 - 1 ) _0(1-z_a)
        + 3 z_a^2 + 1 - z_a + z_a^3/z_a^2]
       
       
       + x̃_a/z_a f'_g(x̃_a/z_a) f_g(x̃_b) [
         - δ(1-z_a) lnΤ^/x̃_a 
         + (1-z_a + z_a^2)^2/z_a^2_0(1-z_a)
       ]
       
       
     - f_g(x̃_a/z_a)  x̃_b f'_g(x̃_b) δ(1-z_a) }
         
       + (a ↔ b)
     .

To bring our result into this same form,
we need to integrate by parts twice, first with respect to Y as shown in sigma_NLP_NLO_2,
and then with respect to z_a.
The details of this calculation are given in ggHg_NLPs.
The final result is shown in sigma_NLL_comparison_4 and is given by the result of
Boughezal:2018mvf in sigma_ggHg_NLL_lit plus an extra contribution,

    σ^(2,1)/Τ^   = σ^NLP<cit.>/Τ^
         
     + 2 C_A/π∫_0^1  x_a  x_b  
       2πδ(x_a x_b ^2 - m_H^2) |_gg→ H^(m_H)|^2/2 x_a x_b ^2
       f_g(x_a) f_g(x_b) e^Y + e^-Y/m_H .

The two results should agree exactly upon integration, and we have not been able to find a source for this discrepancy. 
As discussed in more detail in
results, the numerical comparison with MCFM provides a strong confirmation of our result.
The numerical extraction of the integrated NLL coefficient yields -0.460 ± 0.026, which agrees
well with our analytic predicted value of -0.466 (see table <ref> below).
Dropping the term in the final line of sigma_NLL_comparison would instead predict the value -1.669.
[
We recently received confirmation from the authors of Boughezal:2018mvf that after rechecking their calculation they identified a missing term, and now agree with our result for σ^(2,1)/Τ^.]





§ NUMERICAL RESULTS




In this section we study our results numerically, including the size of the power corrections and the rapidity dependence. We also compare our analytic results for the (α_s) NLP power corrections with the full nonsingular spectrum obtained numerically from the LO V+jet and H+jet calculations in MCFM8 <cit.>.
In Moult:2016fqy,Moult:2017jsg, the NLP corrections were extracted numerically
by using a fit of the known form of their logarithmic structure to the nonsingular spectrum from MCFM8.
In Moult:2016fqy,Moult:2017jsg, these fits were carried out for the leptonic definition.
Here, we have in addition performed the fits also for the hadronic definition.
We find excellent agreement between the analytically predicted values
and the numerically extracted values for all coefficients, i.e., for the LL and NLL coefficients
in all partonic channels for both the leptonic and hadronic definition.
This provides a strong and independent cross check for the correctness of the analytic NLL results obtained here. By comparing the complete nonsingular spectrum with our NLP result, we can also assess the importance of power corrections beyond NLP.

The NLO power corrections for each partonic channel are extracted from the nonsingular spectrum by using the fit function

    F_NLO(τ)
       = /lnτ{τ[(a_1 + b_1 τ + c_1 τ^2) lnτ
    + a_0 + b_0τ + c_0τ^2 ] } ,


with τ≡Τ_0/m_Z for Z production and τ≡Τ_0/m_H for Higgs production. Details of the fitting procedure have been
described already in Moult:2016fqy,Moult:2017jsg, so we do not repeat them here.
A key point is that in order to obtain a precise and unbiased fit result for the to-be extracted a_i coefficients, it is crucial to include the higher-power b_i and c_i terms in fitfun, and to carefully choose the fit range and verify the stability of the fit, as was
done in Moult:2016fqy,Moult:2017jsg. At the level of precision the a_i are extracted,
this is essential since the full nonsingular cross section includes the complete set
of power corrections and if the b_i and c_i terms were
neglected, these higher-power corrections would be absorbed by the a_i terms
in the fit, rendering their numerically extracted values meaningless.
To obtain a precise extraction of the NLL coefficient a_0, we fix the LL coefficient
a_1 in the fit to its analytic result.

The relevant coefficients for our NLP comparison at NLO are the LL coefficient a_1 and the NLL coefficient a_0. For leptonic Τ they
were extracted for Drell-Yan in Moult:2016fqy and for gluon-fusion Higgs in Moult:2017jsg and for the hadronic Τ we have obtained them here.
Depending on the partonic channel, the uncertainties on the fitted
coefficients range from 0.08% to 2.3% for leptonic Τ and from 0.6% to 5.7%
for hadronic Τ. The latter has larger uncertainties because its power corrections are
larger, requiring the fit to be restricted to smaller Τ values where the
uncertainties in the nonsingular data are larger.




 §.§ Drell-Yan Production



We first consider Drell-Yan production, taking pp→ Z/γ^*
at E_ cm = 13. We use the MMHT2014 NNLO PDFs <cit.> with
fixed scales μ_r = μ_f = m_Z, and α_s(m_Z) = 0.118. We fix Q=m_Z,
integrate over the vector-boson rapidity,
and work in the narrow-width approximation for the Z-boson. The NLP corrections
for the leptonic Τ definition were numerically extracted in Moult:2016fqy.
The results for both the leptonic and hadronic definitions for all partonic channels are collected and
compared to our analytic predictions in table <ref>.
We find excellent agreement within the fit uncertainties in all cases.














In fitNLO we show the complete NLO nonsingular contributions as black dots, as well as a fit to their form with the solid red curve. Given the agreement in table <ref> between our analytic a_0 and the earlier fit result for a_0, we have fixed a_0 to the analytic result, and redone the fit using fitfun to obtain this red curve. The red curve from this fit is fully consistent with the earlier fit result from Moult:2016fqy. The dashed orange curve in fitNLO is the extension of the fit function beyond its fit range. In dotted green and dashed blue we show our analytic predictions. We see that with the inclusion of the NLL power
corrections, we obtain an excellent description of the full nonsingular cross
section up to nearly Τ_0 ∼ 1 GeV. This is quite remarkable, and shows
that additional higher-order power correction terms are truly suppressed.

In cumulantNLO we show a plot of the corresponding residual power corrections for
the cumulant, Δσ(τ_), on both a linear scale (left)
and logarithm scale (right). The solid red curve shows the full power corrections,
the solid green curves show the remaining power
corrections after including a_1 in the subtractions,
and the solid blue curve those after including a_1 and a_0 in the subtractions.
We see that with the inclusion of the full NLL power corrections, we achieve
more than a factor of 100 reduction in the residual power corrections as
compared with the leading-power result at NLO. Both partonic
channels have similarly sized power corrections and show a fast convergence of
the power expansion. The fact that the blue curve in the logarithmic plot exhibits
a steeper slope than the red and green curves is due to its τ_^2 scaling
corresponding to a next-to-next-to-leading power correction. This provides a nice
visualization that our results correctly capture the complete NLP contribution.

The analogous results for the fitted nonsingular spectrum and the residual
power corrections Δσ(τ^_) for the hadronic Τ definition
are shown in fitNLOcmcumulantNLOcm.
As expected, the power corrections
are substantially larger for Τ^ than for the leptonic definition. To obtain
similarly sized power corrections, one has to go to about an order of magnitude
smaller values of Τ^. Apart from the overall enhancement, the qualitative
behavior of the LL and NLL contributions and the different partonic channels is
the same. This is expected from our analytic results, which show that the coefficients
for both definitions have essentially the same structure and primarily differ in the overall
factors of e^± Y leading to the rapidity enhancement for the hadronic definition
already observed in Moult:2016fqy, Moult:2017jsg.

In dsigdY_DY we show the rapidity dependence of the NLP corrections at fixed
τ_ = 10^-3 for both leptonic and hadronic Τ normalized to the LO
rapidity spectrum. We can clearly see the exponential enhancement for the hadronic definition at large Y. For the qg channel,
the asymmetric behavior in rapidity is expected from its analytic result.
The result for the gq channel corresponds to taking Y→ -Y, such that their sum
is symmetric in rapidity. While the leptonic definition does not suffer from the
exponential enhancement of the hadronic definition, it still exhibits a substantial increase at
large positive Y in the qg channel, as well as a suppression at
large negative Y.  This is due to the substantially different x-dependence
of the quark-gluon luminosity (and its derivative) compared to the
qq̅ luminosity in the LO result to which we normalize. Knowing the
NLL contribution to the power corrections differential in rapidity enables one
to explicitly account for this effect in the subtractions.







 §.§ Gluon-Fusion Higgs Production



Next, we consider gluon-fusion Higgs production. We take pp→ H
at E_ cm = 13 with an on-shell, stable Higgs boson with m_H = 125,
integrated over all Y. We use the MMHT2014 NNLO PDFs <cit.>, with
fixed scales μ_r = μ_f = m_H, and α_s(m_H) = 0.1126428. The NLP power corrections for this configuration for the leptonic Τ definition were extracted numerically in Moult:2017jsg.
The results for both leptonic and hadronic definitions for all partonic channels
are collected and compared to our analytic predictions in table <ref>.
In all cases, excellent agreement is observed within the fit uncertainties.


















In fitNLO_gg we show as the solid red curve a fit to the full nonsingular
result  at NLO (black points), which is compared with the LL and NLL
predictions in dashed green and dashed blue, respectively. Once again this
solid red fit curve is obtained using the form in fitfun with a_1 and a_0 fixed
by the analytic result in table <ref>, and agrees very well
with the corresponding result obtained in Moult:2017jsg where a_0 was a parameter in the fit.
In all cases, we find that the NLL result provides a good description of the full nonsingular cross section. This is expected since the NLL results includes all NLP terms in the NLO cross section.
We see, however, that particularly for the gq+qg channel, the NLL result for a_0 is required to get a good description, and the LL power correction a_1 alone is not sufficient. Thus the gq+qg channel provides an example where simply looking at the size of the residual nonsingular result after subtracting the a_1 term does not suffice to validate the value of this coefficient.

In cumulantNLO_gg, we show a plot of the corresponding power
corrections for the cumulant, Δσ(τ_), on both a linear scale (left)
and logarithm scale (right). Here we
more easily see that the inclusion of the NLL power corrections significantly
reduces the residual power corrections for the subtractions. For the dominant
gg→ Hg channel at a typical value of τ_∼ 10^-3 approximately one order
of magnitude is gained at each logarithmic order that the power corrections are
computed.
From table <ref>, we see that for the gq+qg→ Hq channel, the LL
coefficient is numerically suppressed, while in contrast its NLL coefficient is
quite larger. Due to this unusual behavior, the NLL result is
required to consistently reduce the power corrections as compared with the
leading-power result. In the qq̅ channel there is no a_1 term, and
significant improvement is apparent from including a_0.

The analogous results for the fitted nonsingular spectrum and the residual
power corrections Δσ(τ^_) for the hadronic Τ definition
are shown in fitNLOcm_ggcumulantNLOcm_gg.
The power corrections are noticeably larger, though the effect of the rapidity
enhancement is not as pronounced as for Drell-Yan, since here the PDFs suppress the
cross section contributions at larger rapidities. For the dominant gg→ Hg channel
there are also numerical cancellations in the NLL coefficient. More precisely the value for a_0
in table <ref> arises as a_0 = 2.356 + (-2.822) = -0.466,
where the first term corresponds to the rapidity-enhanced version of the leptonic
a_0 while the second term is the NLL contribution arising from the additional
rapidity dependence in the argument of the leading logarithm discussed below kernel_Cgg_LL.
As a result of this cancellation, including only a_1 in the subtractions
leads to slightly smaller power correction above τ_ > 10^-3
than subtracting both a_1 and a_0
(compare the green and blue solid lines in the top row of cumulantNLOcm_gg).
If the second NLL contribution were included as part of the LL result, the latter would
provide a much poorer approximation and including the remaining NLL contribution
would provide a substantial improvement.
Either way, the remaining power corrections after subtracting the full NLL result
shows a much steeper slope, which is as expected from its τ_^2 scaling.
This provides another example where considering only the overall size of the improvement
can be potentially misleading. The gq+qg→ Hq channel shows a similarly unusual
behavior as for the leptonic definition.

In dsigdY_H we show the rapidity dependence of the NLP corrections at fixed
τ_ = 10^-3 for both leptonic and hadronic Τ normalized to the LO
rapidity spectrum. The exponential enhancement for the hadronic definition at large Y is again apparent in the LL results. The NLL coefficients again exhibit an enhancement already for
the leptonic definition at large Y. This is again due to the different x dependence of the
quark PDF and the PDF derivatives compared to the LO gg luminosity to which we normalize.
The quark PDF contributions are also the main
reason why the NLL term for the gq channel (a_0 in table <ref>) is much
larger than the LL contribution.
For the hadronic definition, the e^± Y factors from the observable definition turn
out to partially compensate these PDF effects. This is best visible in the gq channel,
where the PDF enhanced terms at negative (LL) or positive (NLL) Y get reduced by a e^± Y factor from the observable definition. The same effect is also present in the gg channel at NLL.
This is the reason why the a_0 term for the hadronic definition in the gq channel
turns out to be even slightly smaller than for the leptonic definition.







§ CONCLUSIONS




In this paper, we have computed the next-to-leading power corrections in the N-jettiness resolution variable for Drell-Yan and gluon-fusion Higgs production at NLO. This builds on our previous work by computing the non-logarithmically enhanced terms at this order. These results enable the performance of the N-jettiness subtraction method to be improved, and provide important information on the structure of subleading power corrections beyond the leading logarithms.
Our calculation is based on a master formula applicable to observables, and highlights a large degree of universality of these power corrections.

We explained in detail the issue of the treatment of Born measurements at subleading power. We have shown that an apparent disagreement in the literature arises due to the fact that the representation used to obtain the power corrections in Boughezal:2016zws,Boughezal:2018mvf is only valid when integrated over all rapidities, and therefore cannot be directly compared with the results of Moult:2016fqy,Moult:2017jsg and those in the present paper, which are differential in rapidity.
We show that after integration over rapidity the LL results agree.
Further details can be found in compare_lit.

We find that the rapidity dependence of the NLL terms is quite sizeable and is therefore important
to know to be able to improve the subtractions.
One reason for this effect is the different x-dependence of parton luminosities
or derivatives of PDFs appearing in the power corrections as compared to the Born-level
parton luminosity. Hence, one can expect this to be a generic feature of subleading
power corrections.

We also compared our analytic NLL results for gluon fusion Higgs production and Drell-Yan to numerical predictions for these NLO power corrections  obtained from a fit to data from MCFM.  In all cases, excellent agreement was found. In addition we studied the extent to which the inclusion of the NLL power corrections improves the subtraction. At NLO, the inclusion of the NLL power corrections completely captures the (τ) terms. Numerically, summing over production channels, the inclusion of these results reduces the size of the power corrections by two orders of magnitude in Higgs production and three orders of magnitude in Drell-Yan production.

There are a number of directions for future work. It will be interesting to extend the calculation of the NLL power corrections to NNLO. Generically we expect up to an order of magnitude improvement could also be obtained by extending the known LL power corrections at this order to NLL.  Beyond fixed order, the derivation of subleading power renormalization group evolution equations at NLL would allow for the all-orders prediction of the NLL terms. Finally, while we have focused here on color-singlet production, our results provide an important step toward the calculation of the NLP corrections at higher orders and for more complicated processes.






We thank Radja Boughezal, Andrea Isgro, Xiaohui Liu, and Frank Petriello for discussions
regarding the results of Boughezal:2016zws, Boughezal:2018mvf.
This work was supported in part by the Office of Nuclear Physics of the U.S.
Department of Energy under Contract No. DE-SC0011090, by the Office of High Energy Physics of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231, and
the Simons Foundation Investigator Grant No. 327942. We thank the Aspen Center for Theoretical Physics, which is supported by National Science Foundation grant PHY-1607611, for hospitality while portions of this work were completed. We also thank DESY and MIT for hospitality.










§ DERIVATION OF NLO LEADING POWER RESULTS




At leading power the singular terms for N-jettiness are most easily obtained from known factorization formulas <cit.>, which describe the singular behavior of the observable to all orders. The fixed order approach of this paper is therefore most useful when such factorization formula are not available, or well understood, such as at subleading power. However, it can also be applied to reproduce the LP results. In this appendix we illustrate this at NLO, by reproducing the one-loop beam and soft functions for beam thrust.

To relate the beam and soft functions as defined in SCET to our calculation in this work,
recall the LP factorization formula for beam thrust <cit.>,

    σ^(0)/ Q^2  Y Τ
     = σ̂^(Q,Y)   H_ab(Q, μ) ∫ t_a  t_b  
       B_a(t_a, x_a, μ) B_b(t_b, x_b, μ) S(Τ - t_a/Q - t_b/Q, μ)
     ,

where the superscript ^(0) refers to LP, and as before x_a,b = Q/ e^± Y.
The hard function H_ab describes virtual corrections to the hard process
a b → L, B_a,b are the two beam functions and S is the soft function.
The beam functions can be further matched onto normal PDFs,

    B_i(t, x, μ) = ∑_j ∫_x^1  z/z_ij(t, z, μ) f_j(x/z, μ)
     .


All of these functions have definitions as field theory matrix elements in the EFT. 
Their fixed order definitions give rise to UV divergences, which are as usual removed
by a renormalization procedure, which in turn gives rise to RGEs that can be used
to resum large logarithms of Τ.
In the approach presented in this paper, the same divergences appear as 1/ϵ IR divergences in the soft and collinear limits of QCD amplitudes.

At LO, we have

    H^_ij(Q^2,μ) = 1
     ,  _ij^(t, x, μ) = δ_ijδ(t)
     ,  
     S^(Τ, μ) = δ(Τ)
     .

At one loop, the convolution structure thus becomes trivial. Working with hard, beam, soft, and PDFs in the bare factorization theorem we have

    σ^(0,1)/ Q^2  Y Τ   = σ̂^  H^_ab(Q^2, ϵ) δ(Τ) f_a(x_a,ϵ) f_b(x_b,ϵ)
     
       
       + σ̂^∑_a'∫_x_a^1  z_a/z_a  Q  ^_a a'(Q Τ, z_a, ϵ) f_a'(x_a/z_a ,ϵ) f_b(x_b,ϵ)
     
       
       + σ̂^∑_b'∫_x_b^1  z_b/z_b  Q  ^_b b'(Q Τ, z_b, ϵ) f_a(x_a,ϵ) f_b'(x_b/z_b ,ϵ)
     
       
       + σ̂^  f_a(x_a,ϵ) f_b(x_b,ϵ) S^(Τ,ϵ)
     .

Note the extra factor of Q in the beam contributions, arising from t_a,b and Τ
having different mass dimensions.
We have written factorization_NLO in a form similar to our master formulas,
such that we can easily read off the one-loop beam function kernels and soft function. The ϵ arguments in factorization_NLO all refer to ultraviolet divergences and can be removed by SCET counterterms to obtain the renormalized factorization theorem. To obtain this it is important to include virtual graphs in the various sectors as well as zero-bin subtractions for the beam functions. 





 §.§ Leading-Power Expansion of Matrix Elements



The leading-power behavior of real emission matrix elements in the soft and collinear limits
is universal, see e.g. <cit.>, and has already been used in universality.
Here, we briefly review the relevant formulas, and give the relevant one-loop expressions.

Given the Born process

    κ_a(q_a) + κ_b(q_b) → L(q_a + q_b)
     ,

where the incoming momenta are given by

    q_a^μ = Q e^Y n^μ/2 ,  
     q_b^μ = Q e^-Y^μ/2 ,

and we write the one-emission process as

    κ'_a(q'_a) + κ'_b(q'_b) → L(q'_a + q'_b - k) + κ_1(k)
     .

In the soft limit k^μ≪ q_a^μ, q_b^μ, the squared matrix element obeys the LP relation

    ^(0)_a' b' → L k(Q,Y; {k})    
     = 16 π^2𝐂/k^+ k^-×δ_a a' δ_b b' ^_a b → L(Q,Y)
     ,

where we made explicit that a soft emission can not change the incoming flavors.
𝐂 = C_F, C_A is the Casimir constant for ab=qq̅, gg.

In the LP n-collinear limit, the particle k arises from the splitting κ'_a →κ_a + κ_1.
If this splitting is allowed, at LP we have (in the notation of master_formula)
q'_a = q_a / z_a and q'_b = q_b, and the LP limit of the matrix element is given by

    ^(0)_a' b' → L k(Q, Y, {k})
     = 8π^2/Q e^Y k^+ P_a a'(z_a, ) ×δ_b b' _a b → L^(Q, Y)
     .

Similarly, in the -collinear limit arising from κ'_b →κ_b + κ_1,
at LP we get q'_a = q_a, q'_b = q_b / z_b,

    ^(0)_a' b' → L k(Q, Y, {k})
     = 8π^2/Q e^-Y k^- P_b b'(z_b, ) ×δ_a a' _a b → L^(Q, Y)
     .

The one-loop splitting functions in d=4-2 dimensions are given by <cit.>

    P_qq(z,)    = C_F [ 1+z^2/1-z - (1-z)]
     ,
    
     P_gq(z,)    = C_F [ 1+(1-z)^2/z -  z]
     ,
    
     P_qg(z,)    = T_F [ 1 - 2z(1-z)/1-]
     ,
    
     P_gg(z,)    = 2 C_A [ z/1-z + 1-z/z + z(1-z) ]
     .

Note that we flipped the notation of qg and gq relative to <cit.>,
following the standard convention.





 §.§ NLO Soft Function




The NLO LP soft function follows from combining sigmaLP with M2_soft_LP_app
using the same steps as in softmaster,

    σ_s^(0,1)/ Q^2  Y Τ   
     = ∫_0^∞ k^+  k^-/(k^+ k^-)^f_a(x_a) f_b(x_b)/(4π)^2 2 x_a x_b ^4(4π)^/Γ(1-)^(0)(Q, Y; {k})  δ[Τ - Τ̂({k})]
    
       
     = ^(Q,Y)/2 x_a x_b ^4𝐂/π f_a(x_a) f_b(x_b)
       ×e^γ_Eμ^2/Γ(1-)∫_0^∞ k^+  k^-/(k^+ k^-)^+1δ[Τ - Τ̂({k})]
    
       
     = σ̂^(Q,Y) 𝐂/π f_a(x_a) f_b(x_b)
       e^γ_E/Γ(1-)2/1/μ (Τ/μ)^-1-2
       
     = σ̂^(Q,Y) 𝐂/4 π f_a(x_a) f_b(x_b)
       [ -4/^2δ(Τ) + 8/1/μ_0(Τ/μ)
       
       
               - 16/μ_1(Τ/μ) + π^2/3δ(Τ) + () ]
     .

The _n(x) are the standard one-dimensional plus distributions, see e.g. <cit.> for details.
Note that there, the precise definition of the MS scheme is important.
We use

    μ^2 ≡μ^2_ MS = 4π/e^γ_E^2
     .

If one were to use μ^2 = (4π)^/Γ(1-)^2, one would miss the π^2/3 term.
For the NLP results presented in the main text, both definitions yield identical results.

Taking sigLPsoft and adding the virtual soft diagram, and then comparing to factorization_NLO, the one-loop bare soft function can be read off as

    S^(Τ,ϵ) =
     𝐂/4 π[ -4/^2δ(Τ) + 8/1/μ_0(Τ/μ)
      - 16/μ_1(Τ/μ) + π^2/3δ(Τ) ]
     .

The finite terms precisely yield the renormalized one-loop soft function <cit.>.





 §.§ NLO Beam Function




Applying the LP master formula sigmaLP to the n-collinear case
and following the same steps as in sigma_NLO_NLP_coll gives

    σ_n^(0,1)/ Q^2  Y Τ   
     = ∫_x_a^1  z_a/z_a f_a'(x_a/z_a) f_b(x_b)/2 x_a x_b ^4Q e^Y/ρ(Q Τ e^Y/ρ)^-z_a^/(1-z_a)^(4π)^/(4π)^2 Γ(1-)
       ×(Q,Y; {k})
     .

Using the universal n-collinear limit M2_coll_LP_app, we obtain

    σ_n^(0,1)/ Q^2  Y Τ   
     = σ̂^(Q,Y)
       ∫_x_a^1  z_a/z_a  f_a'(x_a/z_a) f_b(x_b)  /2πe^γ_E/Γ(1-)1/Τ(Q Τ/μ^2e^Y/ρ)^-z_a^ P_a a'(z_a, )/(1-z_a)^ ,

where P_a a'(z) is the standard ϵ-dependent splitting kernel at NLO.
Comparing to factorization_NLO, we can read off a result that will enable us to obtain the real radiation bare NLO beam function kernel,

    ^(1) real_ij(t, z, ϵ) =
     /2π[ e^Y/ρe^γ_E/Γ(1-)1/μ^2(t/μ^2e^Y/ρ)^-1-z^ P_i j(z, )/(1-z)^]
     .

The splitting function P_ij(z) may contain divergences as z → 1,
which are regulated by the overall (1-z)^-.
All divergences thus arise from the two expansions

    1/μ_ρ^2(t/μ_ρ^2)^--1   =
     - δ(t)/ + 1/μ_ρ^2_0(t / μ_ρ^2)
     - 1/μ_ρ^2_1(t/μ_ρ^2) + ()^2
     ,
    
     (1-z)^-1-   = - 1/δ(1-z) + _0(1-z) - _1(1-z)  + ()^2
     ,

where we defined μ_ρ^2 = μ^2 ρ/e^Y for ease of notation.
As written I_ij does not yet contain the corresponding collinear virtual and zero-bin contributions. 




  
Example: qq Kernel
From Pij, we obtain

    z^ P_qq(z, )/(1-z)^   
     = C_F z^/(1-z)^[ 1+z^2/1-z - (1-z)]
    
       
     =  P_qq(z) + C_F { - 2/δ(1-z) - 3/2δ(1-z)
     
           
       + [ (z-1)  - (1+z^2) _1(1-z) + (1+z^2) ln(z) _0(1-z) ]}
       + (^2)
     ,

where the LO quark splitting function is given by

    P_qq(z) = C_F [Θ(1-z) 1+z^2/1-z]_+ = (1+z^2) _0(1-z) + 3/2δ(1-z)
     .

Adding the corresponding virtual collinear and zero-bin contributions, I_ij yields

    ^(1)_qq(t, z, ϵ) =
      C_F/2π{   2/^2δ(t) δ(1-z)
      + δ(1-z)/[ 3/2δ(t)  - 2 1/μ_ρ^2_0(t/μ_ρ^2) 
     ]
      - δ(t)/P_qq(z)/C_F
       
       + 2 δ(1-z) 1/μ_ρ^2_1(t/μ_ρ^2) +  1/μ_ρ^2_0(t/μ_ρ^2) [P_qq(z) - 3/2δ(1-z) ]
      
       
      + δ(t) [ (1+z^2) _1(1-z) - 1+z^2/1-zln(z)  + (1-z) - π^2/6δ(1-z) ]
     } .

Note that all divergences are proportional to δ(1-z), such that they cancel
after adding the soft,  collinear and the virtual hard contribution from H_ab^ NLO(Q^2,ϵ),
as the latter also has the universal structure (for Drell-Yan)

    σ_virt^(0,1)/ Q^2  Y Τ   = σ̂^(Q,Y) f_a(x_a)  f_b(x_b)δ(Τ)  C_F/π[
       -1/^2 - 1/(3/2 - lnQ^2/μ^2) + () ]
     .

The cancellation of the 1/^2 and the remaining 3/(2) pieces is obvious
from comparing to I_qqS_LP_NLO. The P_qq(z)/ϵ term cancels with the ultraviolet divergence from the bare quark PDF. 
The remaining ln(Q^2/μ^2)/ϵ term cancels when combining the _0(Τ/μ)/ϵ and _0(t/μ^2)/ϵ distribution terms.
The remaining O(ϵ^0) piece in I_qq gives the renormalized beam function
and agrees with the result in Stewart:2010qs.



  
Example: qg Kernel
For the full LP correction to Drell-Yan production, q q̅→ Z,
we also require the quark-gluon kernel.
Here we only need

    z^ P_qg(z, )/(1-z)^   
     = z^/(1-z)^ T_F [ 1 - 2z(1-z)/1-]
    
       
     = P_qg(z) [ 1 + ( lnz/1-z + 1 )] -   T_F + (^2)
     ,

where the finite quark-gluon splitting function is defined as

    P_qg(z) = T_F [(1-z)^2 + z^2 ]
     .

Equation (<ref>) thus yields

    ^(1)_qg(t, z, μ)    =
      T_F/2π{
      - δ(t)/P_qg(z)/T_F  
      + 1/μ_ρ^2_0(t/μ_ρ^2) P_qg(z)
      + δ(t) [ P_qg(z) ( ln1-z/z - 1) + 1 ] } .

Again the P_qg(z)/ϵ divergence cancels against the same mixing term from the bare gluon PDF. 
The remaining O(ϵ^0) terms give the mixing term in the one-loop quark beam function, agreeing with the result in <cit.>.





§ COMPARISON OF NLP CONTRIBUTIONS FOR GG → H G AT NLO




Here, we give the explicit calculation to obtain our full NLP result for hadronic Τ
in the gg → Hg channel in the form of sigma_NLL_comparison.
Our result prior to integration by parts is obtained by inserting kernel_Cgg with ρ = 1 into sigma_NLP_NLO_1,

    σ^(2,1)/Τ^   
     =  C_A/π∫_0^1  x_a  x_b  
       2πδ(x_a x_b ^2 - m_H^2) |_gg→ H^(m_H)|^2/2 x_a x_b ^2∫_x_a^1  z_a/z_a e^Y/m_H
         ×{
        f_g(x_a/z_a) f_g(x_b) [
          (lnΤ^/m_H e^Y + 1 ) δ(1-z_a)
          
       
          - _0(1-z_a) + 1 - z_a + 9z_a^2 - 5 z_a^3 + 7 z_a^4 - 3 z_a^5/2 z_a^2]
       
           
       + x_a/z_a f'_g(x_a/z_a) f_g(x_b) (1-z_a + z_a^2)^2/2 z_a
           
       + f_g(x_a/z_a)   x_b f'_g(x_b) [
          ( - lnΤ^/m_H e^Y - 1  ) δ(1-z_a)
          
       
          + (1+z_a)(1-z_a+z_a^2)^2/2z_a^2 _0(1-z_a)
       ]
       }
       
     + ( a ↔ b  ,  e^Y ↔ e^-Y)
     .

Here, we also separated the pure _0(1-z_a) term from terms regular as z_a → 1.

We now apply integration by parts to the f_g f'_g piece, except for its -δ(1-z_a) term.
In the notation of sigma_NLP_NLO_2, this is achieved by choosing

    D^(2,1)(z_a, z_b, Τ^) = 4 C_A e^Y/m_H[   - lnΤ^/m_H e^Yδ(1-z_a)
     
       
     + (1+z_a)(1-z_a+z_a^2)^2/2z_a^2 _0(1-z_a) ]
     δ(1-z_b)
     .

Here, we only consider being inclusive in Y, so we do not write down the boundary term.
Eq. (<ref>) becomes

    σ^(2,1)/Τ^   
     =  C_A/π∫_0^1  x_a  x_b  
       2πδ(x_a x_b ^2 - m_H^2) |_gg→ H^(m_H)|^2/2 x_a x_b ^2∫_x_a^1  z_a/z_a e^Y/m_H
         ×{
        f_g(x_a/z_a) f_g(x_b) [ 2 δ(1-z_a)
          + ((1+z_a)(1-z_a+z_a^2)^2/2z_a^2 - 1 ) _0(1-z_a)
          
       
          + 1 - z_a + 9z_a^2 - 5 z_a^3 + 7 z_a^4 - 3 z_a^5/2 z_a^2]
       
           
       + x_a/z_a f'_g(x_a/z_a) f_g(x_b) [
          - δ(1-z_a) lnΤ^/m_H e^Y
          + (1-z_a+z_a^2)^2/z_a^2 _0(1-z_a)
          
       
          - (1-z_a)(1-z_a+z_a^2)^2/2z_a^2]
       
           
       - f_g(x_a/z_a)   x_b f'_g(x_b)  δ(1-z_a)
       }
         
     + ( a ↔ b  ,  e^Y ↔ e^-Y)
     ,

where as usual, e^Y = √(x_a/x_b).
Next, we apply the following integration by parts:

    ∫_x_a^1  z_a/z_ax_a/z_a f'_g(x_a/z_a) f_g(x_b)
       [ -(1-z_a)(1-z_a+z_a^2)^2/2z_a^2]
    
    
     =    ∫_x_a^1  z_a/z_a  f_g(x_a/z_a) f_g(x_b)
      2 - 3 z_a + 5 z_a^3 - 6 z_a^4 + 3 z_a^5/2 z_a^2 .

Putting this back into sigma_NLL_comparison_2, we can rewrite it in a form close to sigma_ggHg_NLL_lit,

    σ^(2,1)/Τ^   
     =  C_A/π∫_0^1  x_a  x_b  
       2πδ(x_a x_b ^2 - m_H^2) |_gg→ H^(m_H)|^2/2 x_a x_b ^2∫_x_a^1  z_a/z_a e^Y/m_H
         ×{
        f_g(x_a/z_a) f_g(x_b) [
          2 δ(1-z_a) + ( (1-z_a+z_a^2)^2/z_a^2 - 1 ) _0(1-z_a)
          
       
          + 3 z_a^2 + 1 - z_a + z_a^3/z_a^2]
       
           
       + x_a/z_a f'_g(x_a/z_a) f_g(x_b) [
          - δ(1-z_a) lnΤ^/m_H e^Y
          + (1-z_a+z_a^2)^2/z_a^2 _0(1-z_a)
       ]
       
           
       - f_g(x_a/z_a)   x_b f'_g(x_b)  δ(1-z_a)
       }
         
     + ( a ↔ b  ,  e^Y ↔ e^-Y)
     .

To compare this result to sigma_ggHg_NLL_lit, use the relations

    e^Y/m_H = x_a /m_H^2 ,  lnΤ^/m_H e^Y = lnΤ^/x_a  .



jhep
