Light quark mediated Higgs boson threshold production in the next-to-leading logarithmic approximation

We study the amplitude of the Higgs boson production in gluon fusion mediated by a light quark loop and evaluate the logarithmically enhanced radiative corrections to the next-to-leading logarithmic approximation which sums up the terms of the form αsnln2n−1mH/mq\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\alpha}_s^n{\ln}^{2n-1}\left({m}_H/{m}_q\right) $$\end{document} to all orders in the strong coupling constant. This result is used for the calculation of the process cross section near the production threshold and gives a quantitative estimate of the three and four-loop bottom quark contribution to the Higgs boson production at the Large Hadron Collider.


Introduction
Accurate theoretical predictions for the (inclusive and differential) Higgs gluon fusion cross section are indispensable for the determination with high precision of the Higgs boson couplings [1]. A dominant component of the gluon fusion process originates from Feynman diagrams with a virtual top quark inside the loop. Due to the hierarchy of the top quark and Higgs boson masses, this component can be accurately determined by expanding around the heavy top quark limit [2][3][4]. In this approach, where top quark loops are reduced to effective point-like vertices, gluon fusion cross sections are now known precisely at very high orders in perurbative QCD [5][6][7][8][9][10][11][12][13][14].
With the achieved precision of a few percent, contributions of lighter quarks of a suppressed Higgs Yukawa coupling cannot be ignored. 1 For light quarks, the top quark effective field theory calculations are inapplicable. The relevant Higgs production probability amplitudes need to be computed with their exact quark mass dependence or, alternatively, by means of a systematic expansion around the antithetic asymptotic limit in which the quark mass is vanishing. The exact quark mass dependence for the gg → H amplitude is only known through two loops [15][16][17][18][19][20][21][22]. The two-loop amplitudes for the top-bottom interference in the next-to-leading order cross section [23][24][25] for the production of a Higgs boson in association with a jet have been computed by means of a small quark mass expansion.
In the small quark mass limit the radiative corrections are enhanced by a power of the logarithm ln(m H /m q ) of the Higgs boson to a light quark mass ratio. For the physical values of the bottom and charm quark masses the numerical value of the logarithm is quite large and it is important to control the size of the logarithmic corrections beyond two JHEP07(2020)195 loops. For the gg → H amplitude the leading (double) logarithmic corrections have been evaluated to all orders in strong coupling constant α s in refs. [26,27]. The abelian part of the double-logarithmic corrections for the gg → Hg amplitude of Higgs plus jet production has been obtained in ref. [28]. Though the leading logarithmic approximation gives a qualitative estimate of the QCD corrections beyond two loops, subleading logarithmic corrections can be numerically important. Their computation is necessary to quantifying the theoretical uncertainty estimate. In particular, in the leading logarithmic approximation the numerical predictions vary significantly with values of light quark masses being taken in different renormalization schemes. The logarithmic terms sensitive to ultraviolet renormalization, which cancel the dependence of the amplitude on the renormalization scheme, are formally subleading. Extending the analysis beyond the leading logarithmic approximation is therefore highly desired. In this paper, we present the analysis of the light quark loop mediated gg → H amplitude in the next-to-leading logarithmic approximation which sums up the terms of the form α n s ln 2n−1 (m H /m q ) to all orders in perturbation theory and use this result for the evaluation of the bottom quark effect on the Higgs boson production cross section near the threshold.
The paper is organized as follows. In the next section we discuss the general structure of the leading and next-to-leading logarithmic approximation, derive the factorization and perform the resummation of the next-to-leading logarithmic corrections to a model masssuppressed amplitude of quark scattering by a scalar gauge field operator. In section 3 we apply the method to the analysis of the gg → Hg amplitude mediated by a light quark loop. The result is used in section 4 for the calculation of the bottom quark contribution to the cross section of the Higgs boson production in gluon fusion in the threshold approximation. Section 5 is our conclusion.
2 General structure of the leading and next-to-leading logarithms The logarithmically enhanced contributions under consideration appear in the high-energy or small-mass limit of the on-shell gauge theory amplitudes. To the leading order of the small-mass expansion these are renowned "Sudakov" logarithms which have been extensively studied since the pionering work [29]. The structure of the Sudakov logarithms in the theories with massive fermions and gauge bosons is by now well understood [30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48]. The characteristic feature of the gg → H amplitude however is that it vanishes in the limit of the massless quark m q → 0 i.e. is power-suppressed. The asymptotic behavior of the power suppressed amplitudes may be significantly different from the Sudakov case and attract a lot of attention in many various contexts (see e.g. [26][27][28]). At the same time the logarithmically enhanced corrections to the on-shell (or almost on-shell) amplitudes in the high energy limit are universally associated with the emission of the virtual particles which are soft and/or collinear to the large external momenta. It is instructive first to review the origin and the structure of such corrections to the leading-power amplitudes. In the next section we discuss the asymptotic behavior of the quark form factor to the next-to-leading logarithmic approximation, which will be used for the analysis of the Higgs production in section 3.

Sudakov form factor
We consider the limit of the large Euclidean momentum transfer Q 2 = −(p 2 − p 1 ) 2 and slightly off-shell external Euclidean quark momenta m 2 q ≪ −p 2 i ≪ Q 2 . In this limit the "Sudakov" radiative corrections enhanced by the logarithm of the small ratio p 2 i /Q 2 are known to exponentiate [29][30][31][32][33][34] and to the next-to-leading logarithmic approximation the Dirac form factor F 1 of the quark scattering in an external abelian vector field reads is the strong coupling constant at the renormalization scale µ, γ (1) q = 3C F /2 is the one-loop quark collinear anomalous dimension, and I DL (I SL ) is the double (single) logarithmic one-loop virtual momentum integral. Let us consider the evaluation of the above integrals in more detail. The double-logarithmic contribution is generated by the scalar three-point integral where l is the gluon momentum. For m 2 q = 0 the above integral can be computed by using the expansion by regions method [79][80][81] with the result where the ellipsis stand for the nonlogarithmic terms, ε = (d − 4)/2 is the parameter of dimensional regularization and the contributions of the ultrasoft, collinear and hard regions are given separately (for application of the expansion by regions to the form factor analysis see [38]). Note that the logarithmic contribution can be read off the singularity structure of the ultrasoft and collinear regions where the quark propagators may be taken in the eikonal approximation. As we see the integral I DL generates pure double-logarithmic contribution associated with the overlapping soft and collinear divergences and does not generate any single logarithmic term. On the other hand the double-logarithmic term can be obtained JHEP07(2020)195 directly by means of the original Sudakov method [29]. In this case the propagators in eq. (2.2) are approximated as follows where we introduce the standard Sudakov parametrization of the soft gluon momentum l = up 1 + vp 2 + l ⊥ with Euclidean l 2 ⊥ ≥ 0. The logarithmic scaling of the integrand requires −p 2 1 /Q 2 < |v| < 1, −p 2 2 /Q 2 < |u| < 1 and the additional kinematical constraints uv > 0 has to be imposed to ensure that the soft gluon propagator can go on the mass shell. After integrating eq. (2.2) over l ⊥ we get (2.5) Note that the lower integration limits in eq. (2.5) are determined in the leading logarithmic approximation only up to a constant factor which we choose in such a way that eq. (2.5) reproduces the result of the explicit evaluation eq. (2.2) to the next-to-leading logarithmic accuracy. As it follows from eq. (2.4) the logarithmic contribution to I DL is saturated by the "soft" virtual momentum 2 with l 2 ≈ 0 corresponding to the gluon propagator pole position with the quark propagators carrying the large external momenta being eikonal.
The single-logarithmic one-loop term gets contribution from the part of the vertex diagram linear and quadratic in the virtual gluon momentum as well as from the quark self-energy diagram and can be reduced to the sum of two scalar integrals (2.6) The first integral develops the logarithmic contribution when the virtual momentum becomes collinear to p 1 . In the light-cone coordinates where p 1 ≈ p − 1 and p 2 ≈ p + 2 it takes the following form i 2π where we neglected m 2 q and p 2 i . The integral over l − can then be performed by taking the residue of the quark propagator with the external momentum p 2 which gives with a condition 0 ≤ l + ≤ p + 2 on the second light-cone component of the virtual momentum. Then eq. (2.7) takes the following form (2.9)

JHEP07(2020)195
Since l + ∼ p + the integral over the transversal component of the virtual momentum is logarithmic for l ⊥ ≪ Q and is regulated by m 2 q or −p 2 2 at small l ⊥ . Thus for m 2 q ≪ −p 2 1 with the logarithmic accuracy we get and The β 0 term in the exponent eq. (2.1) sets the scale of the strong coupling constant in the double-logarithmic contribution and is a geometric average of the hard scale Q 2 and the ultrasoft scale p 4 i /Q 2 , as follows from the evolution equations analysis [37,38]. For the analysis of the Higgs boson production we need a generalization of the above result to the quark scalar form factor F S . The structure of the Sudakov logarithms does not depend on the Lorenz structure of the amplitude. However, in contrast to the vector case the scalar form factor has a nonvanishing anomalous dimension. The physical scale for the Yukawa coupling to an external scalar field is given by the momentum transfer which results in an additional dependence of the form factor on Q. Thus to the next-to-leading logarithmic accuracy we have where the exponential factor is identical to eq. (2.1), in the renormalization group factor γ (1) m = 3C F is the quark mass anomalous dimension, and ν is the renormalization scale of the Yukawa coupling.
The light quark mediated gg → H amplitude is suppressed by the quark mass and its high-energy asymptotic behavior is significantly different from the leading-power Sudakov form factor considered above. To determine the structure of the logarithmically enhanced corrections in this case in the next section we consider an auxiliary amplitude of a massive quark scattering by an abelian gauge field operator. This rather artificial amplitude is a perfect toy model which reveals the main features of the problem in the most illustrative way with minimal technical complications and the corresponding analysis can be easily generalized to the other mass-suppressed amplitudes and nonabelian gauge groups.

Massive quark scattering by a gauge field operator
We consider quark scattering by a local operator (G µν ) 2 of an abelian gauge field strength tensor for the on-shell initial and final quark momentum p 2 1 = p 2 2 = m 2 q ≪ Q 2 . A detailed discussion of this process in the leading logarithmic approximation can be found in refs. [26,27]. Below we extend the analysis to the next-to-leading logarithmic approximation. The leading order scattering is given by the one-loop diagram in figure 1(a). Due to helicity JHEP07(2020)195 conservation the corresponding amplitude is suppressed at high energy by the quark mass and with the logarithmic accuracy reads where, α q = e 2 q /(4π), e q is the quark abelian charge and the scalar double-logarithmic integral over the virtual quark momentum reads (2.14) The integral can be evaluated through the expansion by regions with the result where L = ln(Q 2 /m 2 q ) and the ellipsis stand for the nonlogarithmic terms. The singlelogarithmic collinear contribution is given by the integral (2.16) and the ultraviolet divergent contribution reads where ν is the corresponding ultraviolet renormalization scale. For ν = m q the linear logarithmic terms in eq. (2.13) cancel and the renormalized amplitude in the next-toleading logarithmic approximation takes a simple form Figure 2. The two-loop Feynman diagrams for the quark scattering by the (G µν ) 2 vertex (black circle) which contribute in the next-to-leading logarithmic approximation. Symmetric diagrams are not shown.
The integral I ′ DL in eq. (2.14) can be evaluated by using the Sudakov parametrization in the same way as it was done in the previous section for I DL . After integrating over l ⊥ we get [27] 1 where the normalized logarithmic variables read η = − ln v/L, ξ = − ln u/L. As in eq. (2.5) we choose the integration limits in such a way that eq. (2.19) reproduces the result of explicit evaluation eq. (2.15) to the next-to-leading logarithmic accuracy. A characteristic feature of the mass-suppressed amplitude is that in contrast to the Sudakov form factor the double-logarithmic contribution is generated by a soft quark exchange between the eikonal gauge boson lines. The additional power of m q originates from the numerator of the virtual quark propagator which effectively becomes scalar and therefore is sufficiently singular at small momentum to provide the double-logarithmic scaling of the one-loop amplitude.
To determine the factorization structure of the next-to-leading logarithms we consider the two-loop radiative corrections and start with the diagrams figure 2(a,b), which include all the double-logarithmic contributions. Following refs. [26,27] we use a sequence of identities graphically represented in figure 3 to move the gauge boson vertex in the diagram figure 2(a) from the virtual quark line to an eikonal photon line. Let us describe this procedure in more detail. We choose the virtual momentum routing as shown in figure 3(a). To the next-to-leading logarithmic approximation the integration over at least one of the virtual momenta has to be double-logarithmic. We start with the case when this is the virtual quark momentum l, which corresponds to the soft quark l 2 ∼ m 2 q . The integral over the gauge boson momentum l g in this case has the same structure as in the one-loop contribution to the vector Sudakov form factor i.e. has the double-logarithmic scaling when the gauge boson momentum l g is soft and the single-logarithmic scaling when l g is collinear to either p 2 or l. As for the vector form factor the ultraviolet-divergent part of figure 3(a) is cancelled against the corresponding ultraviolet-divergent parts of figures 2(c,d). Then we can decompose the lower quark propagator as follows The second term in the above equation can be neglected if l g is collinear to p 2 or soft. In the former case l − g ≪ l + g (cf. eq. (2.8)) while l ⊥ g factor makes the integral over l ⊥ g nonlogarithmic JHEP07(2020)195 (cf. eq. (2.9)). For soft l g one has l ⊥ g 2 ∼ l − g l + g so the second term in eq. (2.20) is proportional to l − g . Moreover one can neglect l 2 g in the denominator of the second quark propagator which becomes S(p 2 + l g ) ∼ 1/l − g and is cancelled by the l − g factor so that the integral over l g depends on l 2 ∼ m 2 q and m 2 q only and does not generate logarithmic corrections. Note that the above approximation is not valid for l g collinear to l, which will be considered separately. In a covariant gauge only A − light-cone component of the photon field can be emitted by the eikonal quark line with the momentum p 2 , while the emission of the A + and transverse components is suppressed. Thus the interaction of the virtual photon which is soft or collinear to p 2 to the quark line in figure 3(a) can be approximated as follows which is equivalent to the QED Ward identity. The right hand side of eq. (2.21) corresponds to the diagram figure 3(b) where the crossed circle on the quark propagator represents the replacement S(l) → S(l) − S(l + l + g ) and the 1/l + g factor is absorbed into the upper eikonal quark propagator. By the momentum shift l → l−l + g in the second term of the above expression the crossed circle can be moved to the upper eikonal photon line which becomes 1 figure 3(c). The opposite eikonal line is not sensitive to this shift since p − 2 ≈ 0. On the final step we use the "inverted Ward identity" on the upper eikonal gauge boson line to transform the diagram figure 3(c) into figure 3(d) with an effective dipole coupling 2e q p µ 1 to the eikonal gauge boson, where e q is the quark charge. Note that we can replace 2p 1 (l + l + g ) by −(p 1 − l − l g ) 2 in the gauge boson propagator as long as l g ≪ Q since p + 1 ≈ 0. Then for the soft virtual momentum l g one can also use the eikonal approximation for the propagators carrying the external momenta p i . By adding the symmetric diagram and the diagram figure 2(b) we get a "ladder" structure characteristic to the standard eikonal factorization for the Sudakov form factor. This factorization, however, requires the summation over all possible insertions of the soft photon vertex along each eikonal line while in the case under consideration the diagram in figure 1(b) with the soft exchange between the photon lines is missing. This diagram can be added to complete the factorization and then subtracted. After adding the diagram figure 1(b) the integral over the soft gauge boson momentum factors out with respect to the leading order amplitude. The additional diagram JHEP07(2020)195 with effective gauge boson interaction accounts for the variation of the charge propagating along the eikonal lines at the point of the soft quark emission and is characteristic to the leading mass-suppressed amplitudes [26,27]. At the same time for the virtual momentum collinear to p 2 the eikonal approximation can be used for the propagators carrying the external momentum p 1 while the integration over l − g is performed by taking the residue of the S(p 2 +l g ) propagator. Thus the integration over the collinear gauge boson momentum completely factors out in the sum of the diagrams figure 2(b) and figure 3(d), as well as in the symmetric diagrams contribution when l g is collinear to p 1 , without any additional subtraction required for the soft virtual momentum. This can be easily understood since in contrast to the soft emission the collinear one is not sensitive to the eikonal charge nonconservation.
In the next-to-leading logarithmic approximation we also need to consider the case when the virtual quark momentum l is either hard or collinear and the corresponding integral is single-logarithmic while the gauge boson momentum l g is soft and the corresponding integral has double-logarithmic scaling. For hard l the integration over l g trivially factorizes and is double-logarithmic only for the diagram figure 2(b). If l is collinear to p 2 the integral over l g in the diagram figure 2(a) is not double-logarithmic since the soft exchange is between two collinear quark lines. At the same time for l being collinear to p 1 the integral over l g in this diagram is double-logarithmic but the above algorithm with the momentum shift is not applicable since the upper line gauge boson propagator has to be on-shell. However in this case the diagrams figure 2(a,b) already give all possible insertions of the soft gauge boson vertex along the eikonal quark line collinear to p 1 so the integration over the soft momentum l g also factorizes [82].
The factorized soft and collinear contributions together with the external quark selfenergy diagrams figure 2(c) add up to the one-loop contribution to the universal factor which describes the next-to-leading Sudakov logarithms for the amplitudes with the quarkantiquark external on-shell lines [36,38,41] where the coupling constant is renormalized at the scale m q , β 0 = −4/3 and γ (1) q = 3/2 are the corresponding beta-function and the collinear quark anomalous dimension, and the "factorization scale" µ f is the mass parameter of the dimensional regularization used to deal with the soft divergence not regulated by the quark mass. Hence the next-to-leading logarithmic approximation for the amplitude can be written in the following form In this equation the function g(−x) of the variable x = − αq 4π L 2 < 0 incorporates the double-logarithmic non-Sudakov contribution of figure 1(b) with an arbitrary number of the effective soft gluon exchanges [27]. It is given by the integral

JHEP07(2020)195
where the expression in the exponent is the result of the one-loop integration over the soft gauge boson virtual momentum in figure 1(b). The integral can be solved in terms of generalized hypergeometric function and has the following asymptotic behavior at x → +∞ The function ∆ q (−x) in eq. (2.24) accounts for the all-order non-Sudakov next-to-leading logarithmic corrections. In two loops these corrections are generated by a part of the collinear contribution which does not factor into eq. (2.23) and by the renormalization group logarithms proportional to beta-function. Note that to get the next-to-leading twoloop logarithmic contribution the integration over the virtual quark momentum l should be double-logarithmic and we can use the same approximation as for the calculation of the one-loop amplitude. Let us consider the collinear contribution first. As it has been pointed out the contribution of the virtual momentum l g which is collinear to l in the diagram figure 2(a) cannot be factorized into the external quark lines. Since in the doublelogarithmic approximation the soft quark momentum is close to the mass shell l 2 ≈ m 2 q this contribution can be read off the one-loop result for the Sudakov massive quark form factor with the external on-shell momenta l and p 2 . After adding the symmetric contribution and the self-energy diagram figure 2(d) and integrating over l g we get the factor The renormalizetion group logarithms from the vacuum polarization of the off-shell eikonal photon propagators in figure 2(e) and the symmetric diagram read where we set the renormalizarion scale of the gauge field operator and α q in the leading order amplitude to m q . To get the two-loop cubic logarithms the expressions eqs. (2.28), (2.29) should be inserted into the integral eq. (2.25) for x = 0 corresponding to the leading order amplitude which gives Since the soft emission from an eikonal line of a given charge factorize and exponentiate, the higher order next-to-leading logarithmic corrections associated with eqs. (2.28), (2.29) can be obtained by keeping the exponential factor in the integral eq. (2.25) unexpanded and the corresponding contribution to ∆ q reads and erf(x) is the error function. At x → ∞ eq. (2.32) has the following asymptotic behavior In three loops, however, a new source of the next-to-leading logarithms starts to contribute, which is related to the renormalization group running of the coupling constant in the diagram with the effective soft gauge boson exchange, figure 1(c). This diagram is needed to provide the factorization of the β 0 term in eq. (3.3). Note that the diagram includes only the soft gauge boson exchange. The expression for the two-loop subdiagram with the vacuum polarization insertion to the soft gauge boson propagator is given up to the overall sign by the result for the two-loop logarithmic corrections proportional to β 0 to the off-shell Sudakov form factor eq. (2.1) and reads where L µ = ln(Q 2 /µ 2 ) and µ is the renormalization scale of α q in the double-logarithmic variable x. As for eqs. (2.28), (2.29), the corresponding higher order next-to-leading logarithmic corrections are obtained by inserting eq. (2.34) into the integral eq. (2.25) which gives where we introduced the function with the following asymptotic behavior at x → ∞ Note that the leading term of the expansion in 1/x, eq. (2.37), vanishes at the normalization scale µ = Qm q . The complete next-to-leading logarithmic contribution reads We confirm eq. (2.38) in two-loop approximation through the explicit evaluation of the corresponding Feynman integrals as functions of the quark mass with subsequent expansion of the result in m q . A detailed discussion of such a calculation will be published elsewhere [83].

JHEP07(2020)195
Finally we should note that according to eq. (2.1) the collinear logarithm eq. (2.28) and the renormalization group logarithm eq. (2.34) exponentiate while the higher order renormalization group logarithms of the form eq. (2.29) sum up to the usual running coupling constant factors (2.39) Thus we can rewrite eq. (2.24) in a more orthodox form which would naturally appear within an effective field theory analysis and includes also the terms α n q L m with m < 2n − 1. We however prefer to work with the strict logarithmic expansion, eqs. (2.24), (2.38). These equations reflect the general structure of the non-Sudakov next-to-leading logarithmic corrections and will be generalized to the Higgs boson production in the next section.

gg → H amplitude mediated by a light quark
The leading order amplitude of the gluon fusion into the Higgs boson is given by the oneloop diagram in figure 4(a). The dominant contribution to the process comes from the top quark loop and in the formal limit of the large top quark mass m t ≫ m H is proportional to the square of the Higgs boson mass m H . By contrast for a light quark with m q ≪ m H running inside the loop the amplitude is proportional to m 2 q . Indeed, the Higgs boson coupling to the quark is proportional to m q . Then the scalar interaction of the Higgs boson results in a helicity flip at the interaction vertex and helicity conservation requires the amplitude to vanish in the limit m q → 0 even if the Higgs coupling to the light quark is kept fixed. By using the explicit one-loop result the light quark mediated amplitude can be written in such a way that its power suppression and the logarithmic enhancement is manifest where L = ln(−s/m 2 q ), s ≈ m 2 H is the total energy of colliding gluons, and the result is given in terms of the heavy top quark mediated amplitude M t (0) gg→H , which corresponds to a local gluon-gluon-Higgs interaction vertex and has one independent helicity component.
Though apparently completely different, the gg → H amplitude shares several crucial properties with the quark scattering amplitude considered in the previous section: it is suppressed by the leading power of the quark mass due to the helicity flip, the doublelogarithmic contribution is induced by the soft quark exchange and the (color) charge is not conserved along the eikonal lines. Moreover, since the eikonal (Wilson) lines are characterized by the momentum and color charge but not the spin, the factorization structure of soft and collinear logarithmic corrections described in the section 2.2 directly applies to JHEP07(2020)195 the case under consideration. In particular the non-Sudakov double-logarithmic corrections are determined by the diagram figure 4(b) with the effective soft gluon exchange and are described by the same function g(x) with x = (C A − C F ) αs 4π L 2 > 0, where the color factor accounts for the eikonal color charge variation and the opposite sign of the argument is dictated by the direction of the color flow. 3 Thus the next-to-leading factorization formula for the amplitude takes the following form (cf. eq. (2.24)) where the Sudakov factor for a gluon scattering is (see e.g. [84]) and we renormalize the strong coupling constant in the leading order amplitude at the infrared factorization scale µ f = µ. The extra renormalization group factor appears in eq. (3.2) since the leading order amplitude is defined in terms of the quark pole mass while the physical renormalization scale of the quark Yukawa coupling is m H . As for the quark scattering amplitude, the next-to-leading logarithmic term ∆ g (x) gets contributions from the collinear gluon exchange which does not factorize into eq.  figure 4(c,d). In these diagrams the integral over the soft quark momentum l should be evaluated in the double-logarithmic approximation which effectively put the soft quark on its mass shell l 2 ≈ m 2 q . In the diagram figure 4(c) the second quark line carrying the external momentum p 2 is off-shell by the amount (lp 2 ) so that the gluon vertex correction can depend only on the ratio m 2 q /(lp 2 ). However this is the physical gluon vertex with a smooth massless quark limit m q → 0 which means that JHEP07(2020)195 it does not produce the next-to-leading logarithmic contribution. Note that this argument does not work beyond the next-to-leading logarithmic approximation when the soft quark can go off-shell. At the same time the single collinear logarithms associated with the Higgs boson vertex corrections, figure 4(d), do not vanish and can be read off the expansion of the Sudakov scalar form factor eq. (2.12) Note that the renormalization group running of the scalar coupling has alredy been taken into account in the factorization formula eq. (3.2). To get the corresponding all-order corrections the above factor should be inserted into the integral representation of the function g(z) which gives the following contribution 4 to ∆ g 2γ (1) Thus the total result for the next-to-leading logarithmic contribution to the amplitude reads with the following perturbative expansion where L µ = ln(−s/µ 2 ) and the functions g γ (x) and g β (x) are defined by eq. (2.32) and eq. (2.36), respectively. In the above equation the three-loop contribution proportional to the beta-function vanishes when the strong coupling constant in the double-logarithmic variable x is renormalized at the scale µ = m H (m q /m H ) 2/5 , which can be considered as the "optimal" renormalization scale in this case. Abelian part of the series expansion of eq. (3.6) agrees with the result [51] for the Higgs boson two-photon decay amplitude.

Higgs boson threshold production
The perturbative analysis of the total cross section of the Higgs boson production requires the inclusion of the real emission contribution, which is not yet available for the light quark loop mediated process with arbitrary energy of the emitted partons beyond the next-toleading order approximation. However, near the production threshold where z = m 2 H /s → 1 only the soft real emission contributes to the inclusive cross section. For 1 − z ≪ m b /m H the energy of an emitted gluon E g is much less than m b . Up to the corrections suppressed by E g /m b ≪ 1 such emission does not resolve the bottom quark loop and has the same structure as in the top quark loop mediated process, where it is known through the next-tonext-to-next-to-leading order. Thus we can apply our result for the analysis of the higher order corrections to the Higgs boson threshold production.

Partonic cross section near the threshold
The expansion of eqs. (3.2), (3.6) gives where C b = 1 + ∞ n=1 c n and up to four loops in the next-to-leading logarithmic approximation we get 5 The coefficient c 1 agrees with the small-mass expansion of the exact two-loop result (see e.g. [19]). The higher-order leading logarithmic terms have been obtained in refs. [26,27] while the next-to-leading logarithmic terms are new. The n l part of the β 0 term in the coefficient c 2 agrees with the small-mass expansion of the analytical three-loop result [86]. The three-loop leading logarithmic term has been recently confirmed through numerical calculation [87]. Eq. (4.2) corresponds to the coefficient (C A − C F ) 2 /5760 of the L 6 s /z term in eq. (C.1) of [87], which agrees with the numerical value 0.0004822530864 given in this paper. The three-loop next-to-leading logarithmic term however depends on the ultraviolet renormalization and infrared subtraction scheme. We keep the bottom quark as an active flavor since it contributes to the running of the strong coupling constant in the relevant scale interval m b < µ < m H . At the same time in ref. [87] a massive quark is treated separately and is decoupled even for m q ≪ m H . This changes the running of the effective coupling constant as well as the form of the infrared divergent factor eq. (3.3) and makes the comparison of the results not quite straightforward beyond the leading logarithmic approximation. As far as we can conclude the conversion of eq. (4.2) to the scheme used in [87] gives the coefficient (C A − C F )(11C A /9 − 3C F /2 − 2T F )/640 of the L 5 s /z term in eq. (C.1), which is in agreement with the numerical value 0.001736111111 obtained in [87]. The mass dependence of the three-loop amplitude has been also studied by numerical method based on conformal mapping and Padé approximants in ref. [88] but the accuracy of this method in the small-mass limit is not sufficient for a comparison with our result.
The dominant contribution to the cross section is due to the interference of M b gg→H with the top quark loop mediated amplitude. We can define the nonlogarithmic part of the gluon Sudakov factor Z 2 g in such a way that it coincides with the gluon form factor which accounts for the virtual corrections to the gg → H amplitude in the heavy top quark effective theory. After combining it with the soft real emission we get the known effective JHEP07(2020)195 theory expression for the radiative corrections to the top quark loop mediated threshold cross section. Then we can write the above top-bottom interference contribution in the factorized form In eq. where α s is the MS coupling constant renormalized at the scale µ, and the effective theory threshold cross section has the following perturbative expansion where v ≈ 246 GeV is the vacuum expectation value of the Higgs field. The coefficients of the expansion eq. (4.5) are the same as for the top quark loop mediated threshold cross section. The corresponding expressions in terms of delta-and plus-distributions up to n = 2 for N c = 3 and the factorization scale µ f = µ read [92,93] σ 0 = δ(1 − z) , where ζ n = ζ(n) is a value of Riemann zeta-function. Though eq. (4.6) includes terms beyond the next-to-leading logarithmic accuracy we prefer to keep them since the logarithmic expansions for c n and σ n are of quite different nature. Note that in the effective theory cross section eq. (4.5) we can take the massless bottom quark limit and in eq. (4.6) the number of active flavors is n l = 5. By re-expanding eq. (4.3) in α s we get the next-toleading logarithmic result for the leading order (LO), next-to-leading oder (NLO), and the JHEP07(2020)195 next-to-next-to-leading order (NNLO) contributions where L ν = ln(m 2 H /ν 2 ), the normalization factor is is the bottom quark Yukawa coupling renormalized at the scale ν, and we use the expansion The expression for the next-to-next-to-next-to-leading order (N 3 LO) contribution is rather cumbersome and can be obtained from eq. (4.3) by using our result for c 3 and the N 3 LO threshold cross-section from [94].

Hadronic cross section in the threshold approximation
With the partonic result from the previous section we can estimate the bottom quark loop contribution to the hadronic cross section given by where f g (x i ) is the gluon distribution function and s denotes the square of the partonic center-of-mass energy. In addition to the expressions of the last subsection, as numerical input, we use the hadronic cross section coefficients for top quark contributions in the infinite mass limit evaluated in the threshold approximation through N 3 LO as obtained by ihixs2 [95]. The numerical results for the top-bottom interference contribution to the cross section in different orders in α s are presented in table 1 for the following values of input parameters: α s (M Z ) = 0.118, n l = 5, ν = µ f = µ = m H /2, m b (m b ) = 4.18 GeV, m H = 125 GeV. The above choice of µ ensures a good convergence of the series eq. (4.5) for σ eff gg→H+X [6]. At the same time σ eff gg→H+X and C b in eq. (4.3) are separately renormalization group invariant and in general one can use a different value of µ in the series for C b , eq. (4.2). However, the corresponding optimal value µ = m H (m b /m H ) 2/5 ≈ m H /3 is quite close to the one for σ eff gg→H+X and therefore we use the same renormalization scale in both series. Note that in the NLL approximation there is no difference between the pole mass m b and the MS mass m b (m b ) and we use the latter as the bottom quark mass parameter for its better perturbative properties.
In table 1 we also present the result obtained with the threshold partonic cross section retaining full dependence on m b , which is available up to NLO, and the leading logarithmic δσ pp→H+X -1.023 -2.000 Table 1. The bottom quark loop corrections in picobarns to the Higgs boson production cross section of different orders in α s given in the leading logarithmic approximation (LL), the nextto-leading logarithmic approximation (NLL) and with full dependence on m b . All the results are obtained with the threshold partonic cross section at center of mass energy of 13 TeV and renormalization/factorization scales set equal to half the Higgs boson mass. Following the conventions of ref. result obtained with the same σ n coefficients but all the subleading terms in eqs. (4.2), (4.9) being neglected. As we see, both perturbative and logarithmic expansions have a reasonable convergence. In LO and NLO, where the full mass dependence is known, we find that the NLL cross section is within 42% and 3% of the exact result, respectively. The inclusion of the NLL terms is crucial for reducing the scale dependence as it determines the scales of the bottom quark mass, Yukawa coupling and the strong coupling constant in the LL result. For example, the renormalization group invariant product of the Yukawa factor eq. (4.9) and the coefficient C b in NNLO is decreased by about 17% with the scale variation from m H /3 to 2m H in the LL approximation and only by about 5% in the NLL one. The scale dependence of the different orders of perturbative expansion for the threshold cross section in the NLL approximation is shown in figure 5.
We should emphasize that the results in table 1 and figure 5 are obtained in the threshold limit. As discussed, for example, in ref. [57], threshold corrections are not uniquely defined for the hadronic cross section integral. Diverse definitions lead to important numerical differences in the estimate of the cross section. In this paper we have adopted the simplest choice of a flux for the partonic cross section which is given by eq. (4.10). For top quark contribution only, with the same flux choice, the threshold N 3 LO cross section in the infinite top quark mass limit constitutes ∼ 65% of the full cross section. While the threshold contribution may not be adequate for precise estimate of the cross section, it does constitute a physical quantity (in contrast to infrared divergent amplitudes) and can therefore be used to detect whether the large logarithms pose any challenges for the convergence of the perturbative expansion. From the above numerical results we conclude that while in NLO the subleading logarithms are sizable, beyond NLO the logarithmically enhanced corrections are modest and under control.
In figure 6 we plot the ratios of the NNLO and N 3 LO top-bottom interference contribution to the threshold cross section to the NLO result in the NLL approximation. For a wide range of scales the K-factors are in the interval from 1.03 to 1.04. To get an estimate of the total NNLO correction to the bottom quark contribution we can apply the JHEP07(2020)195  where we use the numerical values from table 1. Similar procedure gives the N 3 LO correction of −0.02 pb.
The accuracy of the NLL approximation in LO and NLO is rather good. For a rough estimate of its accuracy in NNLO we can use the sum of the (highly scheme dependent!) subleading L n s /z three-loop terms with n = 0, . . . , 4 in eq. (C.1) of [87]. By applying the corresponding K-factor to the LO result δσ pp→H+X in table 1 we get a correction of 0.049 pb,

JHEP07(2020)195
which constitutes approximately −40% of the NLL correction eq. (4.11), very much like for the LO terms. This tempts us to assign in general a 40% uncertainty to the NLL approximation. However, the analysis of the electroweak Sudakov logarithms [39,41] quite similar to the mass logarithms discussed in this paper suggests that in NNLO the next-to-next-toleading logarithms can be numerically equal to the LL and the NLL terms. This gives us a conservative estimate of 100% uncertainty of the result eq. (4.11). Assuming also 20% uncertainty due to the N 3 LO and higher order corrections together with the 50% uncertainty of the threshold approximation discussed above and adding up the errors linearly we obtain a rough estimate of the bottom quark mediated contribution to the total cross section of Higgs boson production in gluon fusion beyond NLO to be in the range from −0.32 to 0.08 pb. This falls within a more conservative estimate of ±0.40 pb given in [6] on the basis of the K-factor for the top quark mediated cross section and the scheme dependence of the result. The above interval can be further reduced by evaluating the next-to-next-to-leading logarithmic contribution and getting an approximation valid beyond the threshold region. The latter however requires the analysis of the logarithmically enhanced corrections to the hard real emission which currently is not available even in the LL approximation.

Conclusion
We have derived the all-order next-to-leading logarithmic approximation for the light quark loop mediated amplitude of Higgs boson production in gluon fusion. To our knowledge this is the first example of the subleading logarithms resummation for a power-suppressed QCD amplitude. By using this result an estimate of the high-order bottom quark contribution to the Higgs boson production cross section has been obtained in threshold approximation. Despite a large value of the effective expansion parameter L 2 α s ≈ 40α s the corresponding perturbative series does converge. In NLO the next-to-leading logarithmic approximation is in a quite good agreement with the known complete result. For the yet unknown NNLO and N 3 LO corrections we have obtained −0.12 pb and −0.02 pb, respectively. With a rather conservative assessment of accuracy of the next-to-leading logarithmic and the threshold approximations we give a rough estimate of the bottom quark mediated contribution to the total cross section of Higgs boson production in gluon fusion beyond NLO to be in the range from −0.32 to 0.08 pb.
The actual accuracy of the logarithmic and threshold approximations however is difficult to estimate and an exact computation of quark mass effects is therefore expected to be important in consolidating the theoretical precision of the top-bottom interference contribution to the inclusive Higgs cross section. With the computation of the complete three-loop gg → H amplitude in ref. [87] and recent advances for two-loop pp → H + jet amplitudes [96] an exact NNLO result is within reach.