The Higgs-gluon form factor at three loops in QCD with three mass scales

We report on the complete three-loop corrections to the Higgs-gluon form factor in QCD. While previous calculations are based on QCD with a single heavy quark of arbitrary mass, we extend the study to QCD involving two different massive quark flavors. Thereby, the full set of possible Feynman diagrams at three-loop order is taken into account. Employing differential equations for the relevant master integrals, we determine the form factor in terms of analytic expansions. Outside the radii of convergence, we compute high-precision numerical samples over the two-dimensional physical parameter space. Our new findings will enter as virtual corrections the computation of the top-bottom interference in hadronic Higgs-boson production at next-to-next-to-leading order (NNLO) in QCD.


Introduction
After the discovery of the Higgs-boson at the Large Hadron Collider (LHC) in 2012 by the collaborations ATLAS and CMS [1,2], high energy physics has entered the phase of precision measurements of properties of the Higgs-boson.In this context, precise theory predictions for Higgs observables play an indispensable role in the search for footprints of physics beyond the Standard Model (SM).Recent experimental measurements of Higgsboson couplings [3,4] with SM particles of masses ranging over multiple orders of magnitude are in agreement with the SM expectations.In view of the growth of available data as well as the high-luminosity phase of the LHC, the associated statistical uncertainties will decrease, requiring a reduction of theory uncertainties in order to keep up with the experimental precision.
The predominant mode for the production of a Higgs-boson at the LHC is gluon fusion [5].Its total cross section exceeds the cross section of any other production channel by at least one order of magnitude at the typical center-of-mass energies.The effects of massive quarks beyond next-to-leading order (NLO) account for roughly one third of the total theory uncertainty imposed on the gluon fusion cross section.Being a loop-induced process, it is sensitive to physics beyond the SM.Therefore, a firm understanding of the Higgs production cross section is mandatory for current and future precision physics.
One of the building blocks in the calculation of the total cross section is the Higgs-gluon form factor parametrizing the amplitude for the scattering of two on-shell gluons with a possibly off-shell Higgs-boson in QCD.While at the two-loop order, the form factor is known for quarks of arbitrary mass running in the loops, both in numerical and analytical form [6][7][8][9], addressing the exact effects of massive quarks at three loops becomes more challenging for reasons concerning the reduction of symbolic expressions as well as the emergence of new mathematical structures of the underlying integrals.
To overcome these difficulties, the form factor at three loops was first approximated as an expansion for heavy quarks [10,11], which is sufficient for the characterization of top quark mass effects.We note that even at four loops, the expansion of the form factor in the limit of a heavy quark mass in known [12].In order to include the impact of the lighter quark flavors at three loops, the expansion around the large mass limit was matched with the non-analytic terms at threshold, calculated in Ref. [13], and convergence was improved using Padé approximants [14].Subsequently, the full exact result for a single massive quark flavor of arbitrary mass became available by numerically solving the relevant integrals [15].Up to now, the corresponding analytical result has only been feasible for certain subsets of Feynman diagrams [16,17].
In the paper at hand, we are concerned with a new class of Feynman diagrams entering the Higgs-gluon form factor at three loops in QCD.It is important to emphasize that the exact calculation in Ref. [15] takes into account the effect of a single massive quark flavor.However, starting from three loops, Feynman diagrams can develop more than just one closed fermion loop.Hence, in order to completely describe the form factor at three loops with general quark masses, a treatment of the subset of Feynman diagrams incorporating two different massive quark flavors as shown in Fig. 1 is obligatory.In Figure 1: The complete set of three-loop Feynman diagrams relevant to the Higgs-gluon form factor featuring two quark loops.While the quark coupling to the external Higgsboson must be massive, the quark running in the other loop can either be massive or massless.In the present work, we focus on the general case involving two heavy quark flavors with possibly different masses in the loops.
the following we present their computation consisting of a combination of numerical and analytical methods, both of which employ differential equations for the relevant integrals.The paper is structured as follows.We introduce our notation and definitions relevant for the understanding in Sect. 2. Subsequently, the methods used throughout the calculation are reviewed in Sect.3. Our results are discussed in Sect. 4 and we draw conclusions in Sect. 5.

Notation and definitions
We consider corrections to the Higgs-gluon form factor parametrizing the amplitude for the production of a Higgs-boson in the scattering of two on-shell gluons carrying momenta p 1,2 , helicities λ 1,2 , and colors a 1,2 .Since the present work is intended to conclude the study of the form factor performed in Ref. [15], we closely follow the definitions thereof and adopt the notation to a large extent.Hence, we define the bare form factor C accordingly: Since in the SM, the amplitude under consideration is loop-induced, an overall factor of the bare strong coupling constant αs has been factored out.Likewise, the vacuum expectation value, v, originating from the coupling of the Higgs-boson to massive quarks is isolated.
The gluon polarization vectors are normalized in the usual way.
The kinematic configuration of the underlying scattering process gives rise to a single kinematic variable s = 2 p 1 • p 2 coinciding with the virtuality of the potentially off-shell Higgs-boson, which is in the following denoted by m 2 H .The unrenormalized form factor C admits a perturbative expansion in the bare strong coupling constant: Owing to the fact that massive particles are required to mediate the interaction between the massless initial-state gluons and the Higgs-boson, the Higgs-gluon form factor naturally becomes a function of mass parameters.In the paper at hand, we focus on to the calculation of the predominant QCD corrections at three loops retaining the mass of two different quark flavors.In particular, our analysis rests on QCD with n f = n t +n b +n l quark flavors and we refer to the massive flavors as the top quark and the bottom quark, respectively, while the remaining quarks are treated as massless particles.However, in the subsequent discussion, the masses of the top and bottom quark should be understood as variable parameters.In order to keep track of contributions to the form factor stemming from different classes of Feynman diagrams, we refrain from setting the labels n t and n b to their physical values.
With the aid of the introduced flavor-labels indicating the number of top (n t ), bottom (n b ) and massless (n l ) quark loops in the underlying Feynman diagrams, we decompose the three-loop coefficient of the form factor accordingly: Contributions originating from Feynman diagrams with a single closed top or bottom quark loop are captured with the coefficients C (2) t and C (2) b , respectively.As opposed to the corrections to the form factor at lower orders in perturbation theory, the three-loop corrections develop for the first time contributions of the type depicted in Fig. 1, namely, Feynman diagrams with two closed fermion loops.In QCD with a single massive quark flavor, the only possible configurations are given by either assigning the massive flavor to both fermion loops, or allowing for a massless flavor in the second fermion loop while the Higgs must couple to the massive quark.In the decomposition of the form factor in Eq. (2.3) valid for two massive quark flavors, the former contribution is covered by the coefficients C (2) tt and C (2) bb , while the latter one is accounted for in terms of C (2) tl and C (2) bl , respectively.Here, the first letter in the subscript refers to the quark flavor attached to the Higgs-boson.
We point out that the subset of Feynman diagrams featuring one massless fermion loop in addition to the massive quark loop has been addressed in Ref. [16] by means of analytical techniques, such that C (2) tl and C (2) bl are known in terms of harmonic polylogarithms.Furthermore, the complete set of three-loop Feynman diagrams in QCD with a single massive quark flavor has been subject to the exact computation in Ref. [15] based on semi-analytical methods.Hence, all coefficients that appear in the first and second line on the right-hand side of Eq. (2.3) are known either in terms of analytic functions or expansions supplemented with high-precision numerical samples over a one-dimensional parameter space.However, it is important to stress that a complete knowledge of the three-loop QCD corrections, C (2) , in the most general case requires the calculation of terms in the third line of Eq. (2.3) in addition to the aforementioned constituents.The missing contribution arises in QCD with two massive quark flavors.More precisely, a class of Feynman diagrams with two distinct massive quarks running in the loops as illustrated in Fig. 1 emerges from the underlying theory.Thus, the new contribution being proportional to n t n b consists of two building blocks denoted by the coefficients C bt , where in analogy to the other coefficients, the first letter in the subscript indicates the quark flavor that couples directly to the Higgs-boson, whereas the second letter refers to the quark flavor running in the second loop.The calculation of these two coefficients is at the core of the present paper.
First, we note that in contrast to the known coefficients in Eq. (2.3) capturing Feynman diagrams with only one massive quark flavor, the new contributions, C tb and C (2) bt , require the introduction of an additional dimensionless scale.Consequently, the Higgs-gluon form factor becomes a function of two dimensionless variables, which we define as the individual quark masses rescaled with respect to the mass of the Higgs-boson squared: Inspired by the fact that three mass scales enter the parametrization, we dub the coefficients under consideration, C tb and C bt , three-scale coefficients in the remainder of the text.The three-scale coefficients depend on x and y as arguments: bt ≡ C It is obvious, that the set of underlying Feynman diagrams contributing to C (2) bt is a copy of diagrams accounted for in terms of C (2) tb with flavors of the quarks swapped.Therefore, the three-scale coefficients are related by symmetry: bt (y, x) . (2.6) In our calculation, we profit from the fact that knowing one of the two coefficients in the entire parameter space spanned by x and y is sufficient for a conclusive determinations of the Higgs-gluon form factor at three loops.In accordance with Eq. (2.6), the coefficients are related to each other by simultaneously exchanging flavors and variables.
In order to yield a physically meaningful result, the bare form factor must be renormalized by redefinition of couplings and masses.The strong coupling constant is renormalized in the MS scheme with n f = n t + n b + n l flavors.Additionally, we perform the decoupling of the massive quark flavors using the decoupling constant derived in Refs.[18,19].Consequently, the strong coupling constant α s ≡ α (n l ) s (µ) exhibits a dependence on the renormalization scale µ dictated by the β-function for n l massless quarks.The gluon wave function as well as the masses of the heavy quarks are renormalized on-shell.To this end, we recomputed the relevant renormalization constants up to two loops to sufficiently high powers in the dimensional regulator ε.In principle, the on-shell gluon wave function renormalization constant can be taken from Ref. [20] since at two loops, the contribution originating from a second massive flavor adds up linearly to the result with one heavy flavor.Likewise, Z m has been computed in the on-shell scheme in Ref. [21], however, in the presence of a second heavy quark flavor only up to the finite part in the dimensional regulator.The integrals needed for the construction of a mass-counterterm in the present calculation have been solved in Ref. [22] in closed form in ε and the masses.We point out, that the corrections to the Higgs-gluon form factor at one and two loops are provided analytically to sufficiently high orders in ε in Ref. [23], such that the derivation of the UV-counterterm is straightforward.
After UV-renormalization, the form factor still exhibits IR divergences regularized in terms of poles in ε.The IR singularities can be factorized employing the I-operator of Ref. [24], which is explicitly given for the form factor at hand in Ref. [25].Alternatively, the IR singularities may be isolated in the MS factorization scheme [26,27], such that only the poles in ε are absorbed.In Ref. [15], the three-loop corrections to the Higgs-gluon form factor were presented in terms of a finite remainder defined with the aid of the I-operator and the translation to the MS factorization scheme was provided for convenience.
We point out that the finite remainder for the three-scale coefficients considered in this paper is independent of the particular scheme choice.Therefore, we define finite remainders, C tb and C (2) bt , without specifying the scheme choice and implicitly treat the form factor in the MS factorization scheme.Moreover, it is worth to mention that the finite remainders of the corresponding three-scale coefficients are devoid of explicit logarithms of the renormalization scale.

Calculation
The computation of the three-scale coefficients is divided into multiple steps aligned with the workflow of current multi-loop calculations.First, all possible three-loop Feynman diagrams for our process of interest featuring two different heavy quark loops are generated within the framework of DiaGen/IdSolver [28].In doing so, not only the relevant diagrams are generated, but also a graph-theoretical comparison based on the topologies of the corresponding diagrams is performed to allow for the construction of a minimal set of integral families.In parallel, computer code for the matching of scalar expressions with the families of integrals after projection to the form factor is automatically produced.The projection of the sum of Feynman diagrams to the form factor consists of working out the Lorentz and color algebra, which is performed with FORM [29][30][31] and the package color [32], and the matching with the previously defined integral families.Thereby, a linear combination of scalar integrals with coefficients being rational functions in ε, x and y as defined in Eq. (2.4) is obtained.Employing Integration-By-Parts identities [33,34] and the algorithm of Laporta [35], the set of scalar integrals can be reduced to a set of 290 master integrals.We point out that this set comprises all master integrals entering both three-scale coefficients.The reduction to master integrals can efficiently be accomplished with Kira [36][37][38] in combination with FireFly [39,40].We decide to choose a basis of master integrals, such that in all expressions, the dependence on the dimensional regulator factorizes from the dependence on the kinematic variables and masses in denominators of coefficients in front of these master integrals [41,42].
At the heart of the evaluation of the three-scale coefficients lies the efficient computation of the master integrals.In order to tackle the computation thereof, we devise two complementary strategies, one of which is based on deep asymptotic expansions of the relevant integrals and the other one rests on numerical techniques as pursuing a fully analytical approach is out of reach with current technology.Both strategies rely on the method of differential equations [43][44][45][46].We note that in order to set up closed systems of first-order linear differential equations in x and y for the master integrals appearing in the form factor, the basis of master integrals must be extended to 305 integrals.In analogy to the strategy outlined in Ref. [15], we insert truncated ε-expansions for the master integrals into the differential equations and solve for the expansions coefficients, which no longer depend on the dimensional regulator.The truncation order is determined by poles in ε in the form factor and the differential equations, such that the finite part of the three-scale coefficients can be obtained.
The first strategy for solving the master integrals is based on asymptotic expansions and exploits the fact that for physical masses of the heavy quarks and the Higgs-boson, respectively, a hierarchy among the scales is implied.In particular, we consider the full set of master integrals and impose a hierarchy in accordance with Taking advantage of the hierarchy allows us to subsequently expand the master integrals in a diagrammatic language for large m 2 t and the resulting expansion coefficients in small m 2 b using their parametric representation.By successively expanding the master integrals in the double-limit, the dependence on the different scales is factorized, such that only singlescale integrals must be computed.The motivation behind this approach is multi-fold.On the one hand, it significantly reduces the complexity of the calculation in comparison to a fully analytic procedure, which might not even be feasible.On the other hand, there is no need to meticulously consider the analytic continuation of analytic expressions or to rationalize roots, which is usually necessary when internal propagators carry masses.Moreover, the asymptotic expansion attains the form of a power-logarithmic series, such that an efficient evaluation is guaranteed.
In the first step, we apply the asymptotic large mass expansions to the set of master integrals.The corresponding diagrammatic formulation [47][48][49][50] for the asymptotic estimate of an integral F with graphical representation Γ in the limit m 2 t → ∞ reads where γ i denotes an asymptotically irreducible subgraph, {p j } γ is the set of momenta being external with respect to the subgraph, and T γ i refers to a Taylor operator, which expands in all quantities considered small in this context.In this way, the original three-scale integral factorizes into a single-scale vacuum integral and a reduced vertex-function, which is independent of the heavy scale m 2 t .The identification of subgraphs has been implemented in the framework of DiaGen/IdSolver.After application of the Taylor operators on the right-hand side of Eq. (3.1), the expressions can be organized in terms of scalar integrals rendering possible a reduction to master integrals.
It turns out that all master integrals we encounter in the large mass expansions were previously studied in the literature.The vacuum integrals can be at most three-loop integrals, which have been computed in Ref. [51], whereas the most complicated reduced vertex-functions are two-loop integrals investigated analytically, for example, in Refs.[8,52,53].Although the relevant two-loop integrals are, in principle, available in analytic form, we point out that the two crossed vertex-integrals studied in Ref. [53] are expressed in terms of iterated integrals over elliptic kernels.We determined all two-scale integrals by means of differential equations in terms of deep expansions with numerical coefficients in special kinematic points and sampled them over the physical parameter space.
For the sake of presentation of the form factor, we find it more beneficial to further exploit the hierarchy of scales and to analytically expand the vertex-functions that emerge from the asymptotic large mass expansion in the limit of small internal mass.To this end, we employ a variety of different techniques.For most of the integrals, we make use of the method of regions [54][55][56] in accordance with its geometric formulation in Feynman parameter space [57,58].The resulting integrals in different scaling regions are evaluated with the aid of HyperInt [59].Alternatively, the integrals can be treated numerically with Mellin-Barnes methods implemented in the package MB.m [60], where the analytic results are restored via the PSLQ algorithm [61].
However, the sketched method only provides the first few expansion coefficients in the double-limit of large m 2 t and small m 2 b .We would like to extend the depth of the series expansion to high order, rendering possible a precise description of the form factor within a certain radius of convergence.In particular, we seek for a deep asymptotic expansion of the bare three-scale coefficient at x → ∞, with truncation order n and maximum logarithmic power m.Here, the y-dependent expansion coefficients are understood as asymptotic expansions in y → 0 themselves, truncated at k and maximum logarithmic power l.The other three-scale coefficient, C bt , is treated analogously.The individual master integrals admit asymptotic expansions similar to Eq. (3.2).In order to determine the expansion coefficients of the master integrals and in turn the coefficients c nm (y), we exploit the fact that a system of first-order linear differential equations for the set of master integrals provides relations among the coefficients c nm (y) at different values of n and m.Hence, the expansion coefficients can be expressed in terms of a finite set of boundary coefficients, which can be inferred from the leading asymptotic behavior as introduced in Eq. (3.1).Instead of determining the y-dependent expansion A scatter plot displaying generated numerical samples of the three-scale coefficients.The orange line indicates the limit x = y.The regions highlighted in yellow and green can be accessed via expansions as explained in the text.The plot is rotated to match the perspective in Fig. 3. coefficients at the level of master integrals, we found it more advantageous to work out the expansion of the three-scale coefficients multiple times at different rational values for y and reconstruct their analytic form at the level of the form factor coefficients.In practice, we truncate the expansions of C (2) tb and C (2) bt , respectively, at n = 40 and k = 50, which is more than sufficient for the phenomenological application considering m t to take the physical value of the top quark mass and m b to assume the mass of a lighter quark flavor.
It is important to stress that the radii of convergence of the series expansions in Eq. (3.2) and Eq.(3.3) are limited due to the previously imposed hierarchy of scales and the presence of physical thresholds.More precisely, expansions in the limit of large internal mass at x → ∞ are only justified in the region x > 1/4 due to the quark-pair threshold being located at x = 1/4.Moreover, expanding the vertex-integrals encoded in the y-dependent expansion coefficients in Eq. (3.3) implies another restriction related to physical thresholds of the corresponding two-scale remainder in Eq. (3.1).Depending in the particular form factor coefficient, the exact position of the threshold varies, since not all vertex-functions enter both three-scale coefficients.For the coefficient C (2) tb , we observe the appearance of a threshold at y = 1/16 stemming from the aforementioned two crossed vertex-integrals.Thus, the validity of the expansion in y → 0 in restricted to y < 1/16.Contrary to that, the expansion of C (2) bt converges for y < 1/4 due to another quark-pair threshold.The constraint on y can be relaxed by inserting the exact solutions for the reduced vertexfunctions.In order to fully exhaust the potential of the deep asymptotic expansions for x → ∞, we employ our numerical solutions for the relevant vertex functions and sample the form factor coefficients even beyond the restricted radii of convergence in y described before.In Fig. 2, the region covered by expanding C Owing to the previously imposed hierarchy among the mass scales, the deep asymptotic expansions in Eq. (3.2) of the form factor coefficients under consideration have a limited radius of convergence.Although the expansion suffices for the study of mass effects in compliance with the mass hierarchy, another technique must be adopted in order to cover the remaining regions of the two-dimensional parameter space.For this purpose, we decide to utilize the auxiliary mass flow method [62,63] as implemented in the package AMFlow [64] and compute the master integrals at fixed numerical pairs (x, y).Thereby, we not only sample the form factor coefficients in regions of the parameter space beyond the radius of convergence of the aforementioned series expansions, but also validate the quality of the expansions within their radii of convergence.More precisely, the numerical samples generated with AMFlow and analytic methods agree by at least ten digits in the radii of convergence of the expansions.In its default setup, AMFlow dresses the massive propagators with the Feynman causality parameter iη and transports boundary conditions from the limit η → ∞ to η = 0 at fixed values of the pair (x, y).Demanding a precision goal of 30 digits for the individual ε-expansion coefficients of the master integrals, we noticed that more than half of the computation time is spent on the construction of the underlying systems of differential equations in η, hence, on the reduction at fixed (x, y) and symbolic η.We overcome this bottleneck by performing the required reductions only once with symbolic (x, y) and load the constructed differential equations when using AMFlow to compute the master integrals point by point.Thereby the numerical calculation receives a significant speedup of more than 100%.As far as the computation time is concerned, a single phase space point is obtained within approximately 2 hours on a laptop.
In Fig. 2, we visualize the grid of collected numerical samples.In total, the set of master integrals has been probed in more than 2000 different points with AMFlow as indicated by the blue dots.For the construction of the grid, we have compactified the physical parameter space (x, y) ∈ (0, ∞) 2 to the unit square with the following transformation: In particular, pairs in (ρ, σ) ∈ (0, 1) 2 are uniformly distributed.In order to adequately describe possible features close to the thresholds and the limit of small mass, additional samples have been produced in the corresponding regions.

Results
From the numerical samples of the master integrals obtained via AMFlow and the deep asymptotic expansions, we construct the bare three-scale coefficients.Their renormalization and subtraction of IR divergences yielding finite remainders is carried out in accordance with the prescriptions defined in Sect. 2. The observation that all poles in ε of the bare results cancel against those of the counterterms serves as a first consistency check of our calculation.Since the three-scale coefficients are related by the symmetry relation tb is parametrized in terms of x and y as defined in Eq. (2.4).The left plot features the real part and the right plot shows the corresponding imaginary part.In order to visualize various features of the illustrated coefficient, we confine the vertical axis to the indicated ranges. in Eq. (2.6), we show results only for the finite remainder coefficient C (2) tb .The real and imaginary part, respectively, are depicted in Fig. 3. To give insight into different kinematic limits, we sliced the three-dimensional illustrations by fixing one of the two variables to a constant value as presented in several subfigures of Fig. 4.
In the following, we want to discuss the finite remainder in different regions of the physical parameter space and point out some details.To facilitate a firm understanding of the limits accessible via expansions as explained in Sect.3, we provide the asymptotic estimate of the finite remainder C (2) tb including the leading and subleading behavior.For convenience, the exact expansion coefficients are replaced by the corresponding numerical values.In the double-limit x → ∞ and y → 0, the finite remainder is given by with L x = log(x) and L y = log(y).We note that the leading asymptotic behavior starts contributing at constant order in x and y.Furthermore, it is worth to mention that beyond the subleading power, square roots in x and y start to enter the expansion.Their appearance can be traced back to the renormalization of the quark masses, since the expansion of the bare form factors as well as the other renormalization constants do not involve roots.Contrary to that limit, the asymptotic estimate for y → ∞ and x → 0 takes the form tb (x, y) = (0.559207 + 0.288409i) + (0.225378 − 0.239183i) L x − (0.0396296 Here, the finite remainder is power-suppressed in both variables.Hence, the finite remainder vanishes in this limit as expected for a massless quark flavor coupling to the Higgs-boson.As opposed to the previous limit, contributions proportional to square roots in x or y due to renormalization of the quark masses cancel out completely. To the best of our knowledge, the contribution to the Higgs-gluon form factor under consideration has not been addressed in the literature so far.Thus, performing independent cross-checks of our results, for example, in special kinematic limits is possible only in a single case.The only result to compare with is given by the limit of equal quark masses being covered in Ref. [15].Starting from general values in the two-dimensional parameter space, the corresponding limit x = y can be smoothly reached, where we find perfect agreement with previous findings.However, it must be mentioned that this evidence does not provide the strongest check for the correctness of our results, since only the limit but not the neighborhood is probed.In Fig. 4g, the corresponding limit of the finite remainder is depicted.It is not difficult to convince ourselves that the finite remainder vanishes in x = y = 0 due to the well-known subleading power behavior.Furthermore, we note that the imaginary part disappears for large x and y, while the real part approaches the constant −77/1728.Between these boundaries, the finite remainder is a continuously differentiable function and varies smoothly at the physical thresholds.Beyond the limit of equal quark masses, the finite remainder is most conveniently analyzed in the regions of large scale hierarchy, where the expansions in Eq. (4.1) for x ≫ y and Eq.(4.2) for y ≫ x, respectively, are valid.At first sight, it is striking that the finite remainder as presented in Fig. 3 develops a divergent asymptotic behavior for small values of y below the threshold located at y = 1/4 independent of x > 0, whereas on the other side of the threshold for y > 1/4, the finite remainder is relatively flat.According to the asymptotic expansion in Eq. (4.1), the divergent behavior of the finite remainder is of logarithmic origin, while in the opposite limit, the flat landscape can be attributed to the power-suppressed character of the form factor as formulated in Eq. (4.2).As far as the power-suppression is concerned, it is clear that the finite remainder approaches zero for x → 0, being in compliance with the expectation that the form factor vanishes when the Higgs-boson couples to a massless quark.Similarly, the finite remainder goes to zero for y → ∞ and x < ∞.We note that the region for y → ∞ becomes flat only after renormalization.In fact, the bare three-scale coefficient exhibits a logarithmic divergence at the finite order in y being completely canceled via renormalization.
Special care must be taken when considering the transition of the real part for y → ∞ from x ̸ = y to the limit where both masses become infinitely large at the same time.The finite remainder vanishes in the addressed regions, which is underlined with the slices in Figs.4c , 4d and 4e.The Figs. 4a and 4c consistently showcase the vanishing of the imaginary part in agreement with the known limit of equal mass, where x and y approach infinity at the same rate.The corresponding real part converges much slower to the nonzero constant as stated above and seems to experience a sharp turn when one parameter is kept finite while approaching infinity with the other one.The intricate asymptotics of the real part in the neighborhood of the boundary point 1/x = 1/y = 0 are stressed in Figs.4a and 4c.More precisely, in Fig. 4a, we observe that the real part converges to zero for y → ∞ and large but finite x, whereas in Fig. 4c, the real part diverges at large fixed y for x → ∞.With the methods at our disposal, it is not possible to further analyze the precise transition from x = y at infinitely large values to large but finite x or y independently, due to a lack of scale hierarchy in that region.Likewise, the same difficulty applies to the neighborhood of the limit where x and y tend to zero at potentially different rates.The Figs. 4b and 4d highlight the problem at hand.While in Fig. 4b, real and imaginary part uniformly converge to zero for x → 0 as demanded by the limit of equal mass, Fig. 4d displays singular asymptotics for finite x as y → 0.
In contrast to the flat region of the finite remainder, the logarithmically enhanced region seems to induce sizable corrections as far as the three-scale coefficient is considered as a virtual correction in hadronic Higgs production.For small values of y, the corresponding quark flavor running in the loop which is not directly coupled to the Higgs-boson becomes effectively massless.It is important to recall that if we had assigned a massless quark flavor to the second fermion loop from the outset, its effect would have been regularized in terms of poles in the dimensional regulator.Since the singularity is now screened by the mass of the quark, a logarithmic divergence is developed.We note that the logarithms in the mass of the second quark flavor cannot be absorbed by a redefinition of the strong coupling constant.Hence, in the phenomenological application a consistent treatment within the 4-flavor scheme is inevitable, since the logarithms in the mass are supposed to cancel against logarithms stemming from gluon splittings into massive final-state quarks in real corrections by virtue of the KLN theorem [65,66].Clearly, the impact of the logarithmic divergence is maximized for a large separation of scales and for the heavier quark directly coupled to the Higgs-boson to prevent a suppression due to the Yukawa coupling.Logarithmic grows of the finite remainder in the aforementioned regions are visible in the slices in Figs.4a , 4b , 4e and 4f.In this context, it is worth to note that for fixed y, only the real part is subject to the logarithmic growth in the limit of large x.Contrary to the real part, the imaginary part assumes a constant value in the limit as clearly visible in Figs.4b and 4f.This property can also be deduced from Eq. (4.1), in particular, the leading asymptotic behavior.

Conclusions and outlook
We completed the calculation of the Higgs-gluon form factor at three loops in QCD with multiple massive quark flavors.In particular, three-loop corrections to the form factor give rise to a class of Feynman diagrams involving two closed loops of possibly different massive quark flavors.Since the exact three-loop computation of Ref. [15] relies on QCD with a single massive quark flavor, the aforementioned subset of Feynman diagrams was discarded.With the findings presented in this paper, Ref. [15] is supplemented with the missing contribution such that the Higgs-gluon form factor is known for arbitrary configurations of mass parameters.
The contribution originating from Feynman diagrams with two different massive quark flavors is parametrized in terms of two dimensionless variables.In order to describe the contribution to the form factor, we employed two complementary methods both of which rest on the method of differential equations.On the one hand, we constructed deep asymptotic expansions of the underlying master integrals in limits relevant for the phenomenological application.On the other hand, the regions of the parameter space beyond the radii of convergence of the asymptotic expansions were sampled by numerically solving the relevant master integrals with AMFlow.Thereby, the physical region of the parameter space is covered by a combination of analytic expansions and numerical samples.We provide an ancillary file with all results obtained in the course of this work as a supplement to this paper.Additionally, we agree to generate more numerical samples of the form factor upon request.Since our contribution to the form factor is presented as a finite remainder, we refer to the analytic lower-order results known to higher orders in the dimensional regulator [23], which may be used to reconstruct the IR poles or to change the renormalization scheme for the mass parameters.
In combination with the results of Ref. [15], our findings constitute a main building block for the evaluation of the top-bottom interference effects in Higgs production via gluon fusion at NNLO in QCD.In light of the phenomenological application, we want to stress that for a consistent treatment in 4-flavor scheme, the new contribution to the form factor is indispensable as it cancels against large logarithmic corrections due to gluon splittings into massive bottom quarks in real radiation.
Finally, we point out the possibility to apply the methods employed in the course of our work to the calculation of similar form factors, for example, the Higgs-photon form factor or massless quark form factors at the tree-loop level in QCD with two different massive quark flavors.Since the Higgs-photon form factor as well as the massless quark form factors share large subsets of relevant integrals with those encountered in the Higgs-gluon form factor, their evaluation is feasible upon minor modifications.We dedicate the computation to future work.
Figure2: A scatter plot displaying generated numerical samples of the three-scale coefficients.The orange line indicates the limit x = y.The regions highlighted in yellow and green can be accessed via expansions as explained in the text.The plot is rotated to match the perspective in Fig.3.

( 2 )
tb (x, y) at x → ∞ is highlighted in yellow, whereas the corresponding region covered by expanding C(2) bt (x, y) at x → ∞ and utilizing the symmetry relation in Eq. (2.6) is highlighted in green.

Figure 3 :
Figure 3: The coefficient C

Figure 4 :
Figure 4: A set of slices through the three-dimensional visualization of the coefficient C (2) tb at fixed values for x or y as indicated in the individual subfigures.The blue and orange curves represent the real and imaginary part, respectively.