Two-loop bottom mass effects on the Higgs transverse momentum spectrum in top-induced gluon fusion

We compute bottom mass ($m_b$) corrections to the transverse momentum ($q_T$) spectrum of Higgs bosons produced by gluon fusion in the regime $q_T \sim m_b \ll m_H$ at leading power in $m_b/m_H$ and $q_T/m_H$, where the gluons couple to the Higgs via a top loop. To this end we calculate the quark mass dependence of the transverse momentum dependent gluon beam functions (aka gluon TMDPDFs) at two loops in the framework of SCET. These functions represent the collinear matrix elements in the factorized gluon-fusion cross section for small $q_T$. We discuss in detail technical subtleties regarding rapidity regulators and zero-bin subtractions in the calculation of the virtual corrections present for massive quarks. Combined with the known soft function for $m_b \neq 0$ our results allow to determine the resummed Higgs $q_T$ distribution in the top-induced gluon fusion channel at NNLL$^\prime$ (and eventually N$^3$LL) with full dependence on $m_b/q_T$. We perform a first phenomenological analysis at fixed order, where the new corrections to the massless approximation lead to percent-level effects in the peak region of the Higgs $q_T$ spectrum. Upon resummation they may thus be relevant for state-of-the-art precision predictions for the LHC.


Introduction
The transverse momentum (q T ) spectrum of the Higgs boson is one of the most important observables at the LHC. High-luminosity measurements by the ATLAS and CMS detectors promise q T -differential Higgs cross section data with relative experimental uncertainties of eventually only a few percent, see e.g. ref. [1]. This precision can be exploited to discover potential new physics effects from modified (effective) Higgs couplings in the spectrum at low and moderate q T [2][3][4][5] (i.e. q T ≪ m H and q T ≲ m H with m H the Higgs mass, respectively). The shape of the spectrum in the low-q T region, where the peak of the distribution is located, for instance, is sensitive to modifications of the Higgs-bottom Yukawa coupling [2]. To tap the full potential of such new physics analyses the theory uncertainty of the Standard Model (SM) background should be at the same level (or smaller) than the experimental one.
The state-of-the-art theoretical predictions of the Higgs q T spectrum in gluon fusion, which is by far the dominant Higgs production channel at the LHC, already come close to the fewpercent precision goal. They have reached the fiducial N 3 LL ′ + N 3 LO level 1 , i.e. they include fixed-order inclusive as well as fiducial corrections and resummation of large logarithms ∝ ln(q 2 T /m 2 H ) up to third order in QCD [6][7][8][9][10][11]. Recently, even some of the ingredients required to achieve N 4 LL resummation became available [12][13][14]. Small-q T resummation is not only necessary to properly describe the shape of the spectrum, but also improves the prediction for the total fiducial cross section of Higgs production measured by the LHC experiments [6]. This in turn allows e.g. to probe the effective Higgs coupling to gluons. The currently available high-order resummed results for the Higgs q T spectrum [6][7][8][9][10][11] are obtained in the heavy-top limit (m t → ∞) 2 , while all other SM quarks are treated as massless.
Finite top mass effects become relevant for large q T (≳ m t ). The full top mass dependence of the Higgs q T distribution is known at NNLO (= NLO 1 in Higgs + 1jet production) [15]. Regarding bottom mass effects there exists quite some literature, see e.g. refs. [16][17][18][19][20][21][22][23][24][25][26]. Concentrating on the q T ≲ m H region and working in the heavy-top limit we distinguish bottom mass corrections proportional to (at least one power of) the bottom Yukawa coupling y b ∼ m b /m H and those at leading power in 1/m H and thus O(y 0 b ). To the best of our knowledge all available literature is concerned with the former type of corrections.
In gluon fusion these arise from a bottom (instead of a top) loop connecting the produced Higgs boson and the two incoming gluons at the amplitude level. Virtual corrections to the bottom quark mediated gluon-gluon-Higgs (ggH) form factor give rise to large (double) logarithms ∝ ln(m 2 b /m 2 H ), which partly compensate the y b m b /m H ∼ m 2 b /m 2 H suppression. Their systematic resummation has been achieved recently [27]. 3 For the class of corrections where a real gluon is attached to the bottom loop inducing the ggH interaction it is currently unknown how to consistently resum potentially large logarithms of the type ln(m 2 b /q 2 T ) or ln(m 2 b /s) [20,21]. In the intermediate-q T region, where (formally) m b ≪ q T ≲ m H ∼ √ s with s the partonic center of mass energy, the bottom mass corrections proportional to one power of y b (top-bottom interference) are dominant. All of them were computed at NLO in ref. [19] and at NNLO in ref. [18], where only the leading terms of an expansion in m 2 b /m 2 H , m 2 b /q 2 T , m 2 b /s are kept in the relevant two-loop amplitudes [29]. Recently, also the full NNLO result for the top-bottom interference contribution became available [16]. In ref. [17] different heuristic prescriptions to supplement these results with partial or ambiguous resummation of large logarithms at NNLL were studied. The bottom mass effects on the Higgs spectrum at moderate q T were found to be of O(5%) (and negative), while the remaining uncertainties due to unknown higher-order (logarithmic) terms ∝ y b were estimated to be at the one-percent level. The leading electroweak corrections to the spectrum in the q T ≲ m H range are somewhat smaller in size [26].
For a determination of y b from the shape of the spectrum at low q T also the bottom annihilation channel to Higgs production (∝ y 2 b ) plays an important role. This contribution is known to NNLL + NNLO [30] in the "five-flavor scheme", where bottom quarks are included in the parton distribution functions (PDFs) and which effectively corresponds to the m b → 0 limit with fixed y b . Beyond this approximation, i.e. consistently assuming Λ QCD ≪ m b ≪ m H , the resummed q T spectrum in bottom quark annihilation can be calculated in full analogy to the "primary mass effects" in the Z-boson q T spectrum (using four-flavor PDFs) following ref. [31].
In the present paper we consider bottom mass corrections to the Higgs q T spectrum at leading power in 1/m H ∼ 1/s. In contrast to the corrections discussed above these are independent of the bottom Yukawa coupling y b ∼ m b /m H and insensitive to the structure of the ggH interaction mediated by a top loop, i.e. we can safely work in the heavy-top limit. The q T -dependent contributions of this type can be written as a series of (m b /q T ) 2n terms (with n ∈ N) and first appear at NNLO in the spectrum. Hence, they come with an additional factor of the strong coupling α s compared to the NLO contributions ∝ y b . Near the peak of the Higgs q T distribution at q T ≈ 2m b ≈ 10 GeV, both types of effects may therefore be of similar size. The aim of this work is to compute the leading O(y 0 b ) bottom mass corrections in the regime Λ QCD ≪ m b ∼ q T ≪ m H and to provide a first analysis of their numerical impact on the NNLO (= NLO 1 ) spectrum.
Cross sections for sufficiently inclusive measurements in high-energy processes that involve largely different (energy) scales can often be shown to factorize to a good approximation. This means that the physics at the various scales can be described independently by separate factorization functions, which typically simplifies the calculation substantially. The factorization of the q T -differential cross section of color-singlet production in the presence of a massive quark flavor was worked out in ref. [31] in the context of the Drell-Yan process. For the (four) relevant hierarchies between the hard scale Q set by the invariant mass of the color singlet, its transverse momentum q T ≪ Q, and the quark mass, generically denoted by m, the appropriate factorization theorems were formulated there using soft-collinear effective theory (SCET) [32][33][34][35][36][37]. 4 This factorization framework also applies to Higgs production in gluon fusion with m = m b and allows to systematically resum all types of large logarithms at leading power in the small scale ratios. The factorized cross section for the regime q T ∼ m ≪ Q takes a special role in this approach, because the factorization theorems for the adjacent regions, i.e. q T ≪ m ≪ Q and m ≪ q T ≪ Q, represent its large and small mass limits. The latter can therefore be derived together with the corresponding power corrections by an expansion of the former in q T /m and m/q T , respectively.
The only missing (and arguably most complex) ingredients to compute the bottom mass effects on the Higgs q T spectrum we are interested in are the transverse momentum dependent (TMD) gluon beam functions for q T ∼ m ≪ Q (= m H ). These factorization functions describe the initial state radiation collinear to the proton beams in the gluon fusion process at the LHC. 5 It is the main purpose of this work to calculate the quark mass corrections to the gluon TMD beam functions to NNLO, while their massless version is already known to N 3 LO [42,43]. This will enable the small-q T resummation in the spectrum with full bottom mass dependence at NNLL ′ and N 3 LL level (and leading power in 1/Q). The beam functions are independent of the hard scattering process. Our results therefore not only contribute to the Higgs q T distribution at hadron colliders, but also to many other TMD cross sections. In particular, all analytic expressions in this paper directly carry over to the transverse momentum spectrum of any color singlet final state produced by gluon fusion.
The paper is structured as follows: In sec. 2 we discuss the factorization theorem for the q Tdifferential color-singlet production cross section and its ingredients in the regime q T ∼ m ≪ Q. We also give some details on the renormalization group (RG) evolution of the TMD beam functions. The two-loop calculation of the quark mass corrections to the gluon beam function is presented in sec. 3. The renormalized results are derived and summarized in sec. 4. In sec. 5 we cross check our results containing the full dependence on m/q T with known expressions in the small and large mass limits, m ≪ q T and m ≫ q T . In sec. 6 we analyze the numerical impact of the computed NNLO bottom mass corrections to the Higgs q T spectrum. We conclude in sec. 7.

Factorization with massive quarks
As the prototype of a TMD observable we consider in this work the q T spectrum of a color singlet state X with invariant mass Q produced in proton-proton collisions. Using this process as an example we discuss in the following factorization in SCET with an active heavy quark flavor of mass m and n l massless quark flavors in the regime where Λ QCD ≪ q T ∼ m ≪ Q. The relevant effective field theory (EFT) modes in this kinematic region are n a -collinear, n bcollinear, and soft. They are defined by the scaling of their typical momenta: using the (light-cone) notation with opposite light-like beam directions n a , n b andn a ≡ n b ,n b ≡ n a , n 2 a =n 2 a = 0, n a ·n a = 2. In addition to the perturbative modes in eq. (2.1) we also define nonperturbative ultra-collinear modes with the momentum scaling (Λ 2 QCD /Q, Q, Λ QCD ) and (Q, Λ 2 QCD /Q, Λ QCD ) to describe the incoming protons (inside the two beams with radius ∼ 1/Λ QCD ) and their constituents. The associated ultra-collinear fields are part of a SCET with n l massless quark flavors, where the heavy quark field has been integrated out. The matching between the two SCET versions with and without massive quarks yields the beam function matching coefficients we compute in this paper. For a detailed account on the EFT framework for the factorization including other possible kinematic regimes with different hierarchies between the scales q T , m, Q as well as the connections among them we refer to ref. [31].
The mode setup in eq. (2.1) is of SCET II type, because soft and collinear degrees of freedom have parametrically the same invariant mass ∼ m 2 ∼ q 2 T and are only separated in rapidity. As a consequence quantum corrections to soft and collinear operators will in general generate rapidity divergences, which require renormalization and eventually manifest themselves as large rapidity logarithms ∼ ln(m/Q) or ∼ ln(q T /Q) in fixed-order predictions of physical observables. We will employ the MS-type approach to rapidity renormalization devised in refs. [40,44]. It allows to systematically resum the rapidity logarithms by means of rapidity renormalization group equations (RRGEs). The corresponding rapidity renormalization scale is denoted by ν, whereas µ represents the standard virtuality-type MS renormalization scale.
The factorization theorem for the gluon fusion process we are concerned with is given in analogy to the (quark-antiquark initiated) Drell-Yan process studied in ref. [31] by where and Y is the rapidity of the color-singlet state. Here and in the following the superscripts {n l + 1} and {n l } on the m-independent factorization functions indicate whether the associated operators belong to SCET with n l + 1 or n l active quark flavors. Their renormalization group (RG) evolution (w.r.t. µ) is performed in the same flavor scheme above and below their characteristic (matching) scale. Inside these functions the running QCD coupling α s (µ) must be consistently evaluated in the respective {n l +1} or {n l } flavor scheme in order to avoid large logarithms when µ is of the order of the characteristic scale. In eq. (2.3) this applies to the hard function H {n l +1} gg governed by the hard scale Q and the parton distribution functions (PDFs) f {n l } k governed by the hadronization scale Λ QCD . The hard function is process-dependent but independent of the observable (here in particular q T ). At leading power in m/Q it corresponds to the squared coefficient of a SCET current operator resulting from the matching between QCD and SCET carried out with n l + 1 massless active quark flavors at µ ∼ Q. Explicit expressions of the hard function for gluon fusion Higgs production (gg → H) can be found up to NNLO e.g. in ref. [45] and at N 3 LO (in the large top mass limit) in ref. [46].
The m-dependent factorization functions in eq. (2.3) are the TMD gluon beam function kernels I µν gk and the gluon-fusion TMD soft function S gg . They can be regarded as matching coefficients at leading power in Λ QCD /m and are located at the flavor threshold where the massive quark is integrated out, i.e. their matching scale is ∼ m. Accordingly, the RG evolution above and below µ ∼ m is performed with n l + 1 active flavors (concerning the running of soft and beam functions) and n l active flavors (concerning the PDF running), respectively. The renormalized m-dependent functions I µν gk and S gg can be expressed in terms of α s (µ) in either the {n l + 1} or the {n l } flavor scheme without introducing large logarithms for µ ∼ m. For this reason we did not assign flavor scheme superscripts to I µν gk and S gg in eq. (2.3). For concreteness we will, however, by default use α {n l +1} s (µ) in explicit expressions of these functions and indicate this by superscripts as I µν{n l +1} gk and S {n l +1} gg when necessary. The massdependent TMD soft function for gluon fusion processes S gg is up to three-loop order related to the corresponding soft function S qq in the quark-antiquark channel by Casimir scaling. At NNLO it can therefore be directly obtained from the result for S qq computed in ref. [31] for Drell-Yan processes by replacing the quadratic Casimir coefficient C F → C A . We give the explicit expression in app. C.5.
In this work we are mainly concerned with the matching coefficients I µν gk of the TMD gluon beam functions B µν {n l +1} g onto the standard collinear PDFs f {n l } k . The leading power matching relation was used to formulate eq. (2.3) and reads (following refs. [47][48][49]) with µ ∼ q T ∼ m and ν ∼ Q representing the virtuality and rapidity matching scales, respectively. Here and in the following we use the symbol ⊗ z for the (Mellin-) convolution (2.6) The parton indices q andq in the sums in eqs. (2.3) and (2.5), stand for all massless quark and corresponding antiquark flavors: q = u, d, s, . . . We will use Q as the index for the massive quark flavor in the following and g represents the gluon. The universal perturbative matching kernels I µν gk describe the n a,b -collinear initial-state radiation characterized by the beam function matching scales for the case of a gluon entering the hard scattering process.
The resummation of logarithms ln(q 2 T /Q 2 ) and ln(m 2 /Q 2 ) is accomplished by performing the (R)RG evolution of the different factorization functions in eq. (2.3) from their characteristic scales to common renormalization scales µ, ν. The resummation kernels for each factorization function (not shown in eq. (2.3) for compactness) contain the resummed logarithms and are obtained by solving the corresponding (R)RGEs. The µ evolution of the TMD gluon beam function is determined by the RGE Similar RGEs hold for the hard and soft functions in eq. (2.3). RG consistency in eq. (2.5) implies that the dependence on the matching scale µ of the coefficients I µν gk is subject to gg (⃗ p T , m, µ, ν) . (2.10) The symbol ⊗ ⊥ denotes the convolution 6 in two-dimensional transverse momentum space. The ν-independence of the cross section in eq. (2.3) implies the RG consistency condition The m-dependence of the TMD beam function matching coefficients I µν gk is currently unknown. Their quark mass dependent contributions are the only missing pieces in eq. (2.3) at NNLO and will be computed in the present work to this order, i.e. O(α 2 s ). The tensor structure of the I µν gk can be decomposed as Note that for unpolarized proton beams the beam function matching kernels I gk and J gk depend on ⃗ p T only via ⃗ p 2 T due to rotation symmetry. Only I gk acquires a non-zero contribution from tree-level matching: 7 (2.14) Here and in the following the superscript (n) with n = 0, 1, 2, ... indicates an n-th order contribution in the perturbative (loop) expansion, i.e. F (n) ∼ α n s for (fixed-order) functions F = B g , I gk , S gg , etc. and γ (n) ∼ α n+1 s for anomalous dimensions γ = γ Bg , γ ν,Bg , γ f,ij , etc. The tensor structure of J gk is orthogonal to g µν ⊥ (upon contraction of all Lorentz indices). Thus, for cross sections that are insensitive to the gluon polarizations like eq. (2.3), only the two-loop matching coefficients I (2) gk are required at NNLO, or N 2 LL ′ and N 3 LL when including resummation. Moreover, there are no quark mass corrections to the one-loop coefficient J (1) gk . Hence, at NNLO (N 2 LL ′ , N 3 LL) accuracy, we only have to consider mass effects on the "unpolarized" gluon beam function B g ≡ g µν B µν g with matching kernel I gk = g µν I µν gk . We present their calculation in the next section.

Calculation of the two-loop beam functions
The SCET operator matrix element defining the bare unpolarized TMD gluon beam function, which accounts for the effects of the n-collinear initial-state radiation (n = n a , n b ), reads where p n (p − ) denotes the incoming spin-averaged proton with lightlike momentum p µ = p − n µ /2 and x ≡ ω/p − . The operator B µc n⊥ ≡ 2 tr[T c B µ n⊥ ] is the gauge-invariant n-collinear gluon field strength in SCET: The SCET label momentum operators ⃗ P n⊥ and P n ≡n · P n [34] act on the n-collinear gluon fields A µ n to their right. For more details on the involved SCET operators and Wilson lines (W n ) we refer to refs. [35,49,51].
The beam function kernels I gk are in practice computed from a perturbative matching calculation of partonic beam functions B g/j obtained by replacing the incoming proton with parton states (with the same momentum) onto corresponding partonic PDFs f k/j [49,51]. The operator between the external states in eq. (3.1) is local in time. We can therefore evaluate the corresponding real-emission Feynman diagrams in fig. 1b and fig. 2 directly as loop diagrams without cutting (or taking a discontinuity or imaginary part). 8 This is in analogy to the SCET calculation of the standard PDFs [45,51], but in contrast to that of the virtuality-dependent beam functions [45,[51][52][53]. Apart from the vertices for the B µc n⊥ operator insertions [45,53] the usual (time-ordered) QCD Feynman rules can be used [33].
The bare partonic beam functions are ultraviolet (UV) and infrared (IR) divergent due to the separation of the collinear regions from the hard and ultra-collinear regions in virtuality. As argued in ref. [54], for the case of the unpolarized gluon beam function UV and IR divergences can be regulated by conventional dimensional regularization upon replacing in (the real emission contribution to) eq. (3.1) with d = 4 − 2ϵ. For the strong coupling α s we employ the MS and for the (bottom) quark mass m the on-shell renormalization scheme. We adopt the notation for the n-loop heavy-flavor correction to the partonic beam function, and analogously for I ν,B , etc. In both parts we let α s evolve with n l + 1 active flavors, i.e. α s ≡ α {n l +1} s (µ) throughout this paper, unless indicated otherwise.

Rapidity regulator
To regulate the rapidity divergences present in the real emission as well as the purely virtual contributions to B g (⃗ p T , m, x) we choose the symmetric Wilson line regulator of refs. [40,44].
The same rapidity regulator has been used in the calculation of the NNLO TMD soft function S gg for m b ̸ = 0 [31] and the massless NNLO TMD beam functions in ref. [50] with which we will combine the massive quark contributions to be computed in the present paper. 9 It may be implemented by modifying the n-collinear Wilson lines as 10 Logarithmic rapidity divergences manifest themselves as 1/η poles in loop (and phase space) integrals. The parameter w obeys the RRGE ν d/dν w = −η w/2 and is set to one after the derivation of the ν anomalous dimension. To obtain the correct anomalous dimensions it is crucial to take the limit η → 0 before ϵ → 0 [40].
Because the relevant rapidity-divergent two-loop diagrams for B (2,h) g/j in fact correspond to one-loop graphs, where either a gluon line is dressed with a massive quark bubble or a triple gluon vertex is replaced by a massive quark triangle subgraph, only single 1/η poles occur in our calculation. This is consistent with the requirement that the rapidity divergences of beam and soft functions cancel in the cross section eq. (2.3) at NNLO. Moreover, since we have at most a single gluon attached to a Wilson line, the n-collinear "group momentum" operator n · P g [40] can be replaced with the standard label momentum operator P n for our purposes.
For technical reasons we will furthermore employ the "δ-regulator" of ref. [63] in the calculation of the virtual diagram in fig. 3a. This corresponds to assigning an offshellness to the involved eikonal (Wilson line) propagators. Although diagram 3a is rapidity-finite, we introduce an auxiliary δ-regulator to avoid rapidity divergences in (intermediate) expressions generated by an integration by parts (IBP) reduction to master integrals, see sec. 3.4. In the context of IBP reduction the δ-regulator proves more efficient than the η-regulator, which modifies the power of eikonal propagators to non-integer values, in the sense that it yields a smaller set of master integrals.
Before applying them in our calculation let us point out an interesting peculiarity of the rapidity regulators of the η-and δ-type. Consider the following rapidity-finite one-loop integral where M ∼ m denotes some mass parameter (and we suppress the causal i0 prescription). An integral of this type for example contributes to the unregulated virtual diagram in fig. 3b, where it arises after integrating the massive quark bubble (leaving a single-parameter integral, see 9 NNLO results for the massless gluon TMD beam functions obtained with different rapidity regulators are found in refs. [55][56][57][58], see also ref. [59]. 10 In general, there are exceptions from this prescription. Consider, for example, the real quark-antiquark cut and the real gluon cut of diagram 2i, which are separately rapidity divergent (as x → 1). In our calculation of Bg(⃗ pT , m, x) the sum of the cuts vanishes exactly, such that the diagram in total does not contribute. For measurements that put different weights on the one-and two-particle final states like the ones in refs. [60][61][62], however, the contributions from fig. 2i must cancel with similar terms from diagrams 2g and 2h related by gauge symmetry (just like the corresponding rapidity divergences in the soft function). This cancellation must be preserved by the rapidity regulator. In practice this means that we have to add (i.e. cancel) these particular terms before regulating the involved diagrams by different powers of |n · Pg| −η . See app. B.1 for a zero-bin calculation where such cancellation is crucial.
sec. 3.2) and canceling a factor of k − ≡n · k in the numerator and the (Wilson line propagator) denominator. Implementing the η-regulator in this diagram according to eq. (3.5) yields Similarly, with the δ-regulator we have 11 In both cases the limit of vanishing regulator is not continuous: While this does not pose a conceptual problem (as long as the rapidity regulator is correctly implemented such that the combined soft and collinear contributions to an observable reproduce the leading-power full-theory result at fixed order), it forces us to consistently regulate also rapidity-finite terms which may complicate their calculation in practice.
The discontinuous behavior of the rapidity-regulated integral in eq. (3.9) is caused by the absence of an external minus momentum component (here p − ) in the denominator of the integrand. For example, the η-regulated rapidity-finite integral ( has a smooth η → 0 limit 12 , and analogously using the δ-regulator: On the other hand, if the scalar loop integral in eq. (3.6) corresponds to a term of an n-collinear SCET diagram like fig. 3b the independence of the large lightcone momentum component p − ∼ Q gives rise to a non-vanishing soft zero-bin (aka soft-bin) [63,64]: Adopting soft scaling of the loop momentum, i.e. k µ ∼ (m, m, m), and expanding the integrand (with M ∼ m) accordingly leaves the integral unchanged. Hence, subtracting the zero-bin exactly removes the contribution of eq. (3.6) from the diagram. This cancellation is not affected by the rapidity regulator. Indeed, we find by explicit calculation that after zero-bin subtraction (see app. B) all rapidity-finite contributions to the bare two-loop beam function have a smooth η → 0 limit. We conjecture that this holds true for any rapidity regulator of η-and δ-type and for any collinear matrix element at any loop order, such that one can always safely drop the rapidity regulator in rapidity-finite terms when consistent zero-bin subtractions are performed.

Dispersion relation for massive bubble diagrams
The two-loop Feynman diagrams in fig. 1b, fig. 2e-i, and fig. 3b-d contain the one-loop off-shell gluon self-energy consisting of a massive quark bubble as subdiagram. For their evaluation in general covariant gauge (ξ = 0 corresponds to Feynman gauge) we conveniently employ the dispersion relation [65,66] Note that the gluon propagators to the left and right of the one-particle irreducible (1PI) gluon self-energy bubble are included in eq. (3.12), color indices are suppressed, and we have introduced the bookkeeping parameter κ ≡ 1 in the second line for convenience of presentation, see below. The first term in the second line represents a weighted integral over a massive gluon propagator (in Landau gauge), the second term is proportional to a massless gluon propagator. The vacuum polarization function due to a virtual massive quark pair is defined by whereμ ≡ µ e γ E /2 (4π) −1/2 and α s ≡ α s (µ) is the running coupling in the MS scheme. Note that the first term in the second line of eq. (3.12) corresponds to the insertion of the onshell renormalized vacuum polarization function and is thus UV finite for given p µ . The UV divergence of the massive quark bubble is contained in the second (massless) term. Using eq. (3.12) we can write each of the two-loop diagrams with an off-shell massive quark bubble as sum of two parts. One part corresponds to the one-loop diagram where the dressed gluon propagator is replaced by a massive gluon propagator with mass M ≥ 2m, which must be integrated over. The other part equals the corresponding massless one-loop diagram times a factor −Π (1) (0, m 2 ). In the virtual diagrams fig. 3b-d the latter contribution vanishes because the loop integral is scaleless. Of course we can also use 13 This leads to one-loop-type integrands including a massive propagator denominator to the power of ϵ with mass m/ y(1 − y), see first line of eq. (3.14). The integration over y, just like the integration over M , is conveniently performed after the loop integration. In our beam function calculation we used both methods and checked that the results agree for all two-loop diagrams with a massive offshell bubble. The approach based on the dispersion relation in eq. (3.12) allows a particularly transparent discussion of the main features of the relevant two-loop diagrams, since their calculation is effectively reduced to a one-loop problem with a massive gluon and integer powers of propagator denominators. The integration over M does not affect important properties of the original two-loop graph like the presence of a non-vanishing zero-bin, a rapidity divergence, or the gauge-dependence. We will therefore mainly refer to this method in the presentation of our beam function calculation.
In refs. [31,[65][66][67] it was argued that the terms (∝ p µ p ν ) labeled by κ in eqs. (3.12) and (3.15) cancel among the two-loop diagrams contributing to gauge-invariant SCET matrix elements such as soft functions or quark jet and beam functions. The statement also holds for the gluon beam function in eq. (3.1). For the real-emission diagrams this can be understood from the analogy to the cancellation of the terms linear in the gauge parameter ξ within the (massless) one-loop calculation of B (1) g/g (or resorting to a Ward identity). Note that (one/two-particle) real-emission and purely virtual contributions are separately gaugeinvariant, i.e. independent of ξ (and thus κ). It is straightforward to explicitly verify that κ drops out separately in the real-emission diagrams 2e and 2f as well as in the sum of diagrams 2g-i already at the integrand level. 14 For the purely virtual diagrams in fig. 3b-d the analogy to the ξ terms of the corresponding massless one-loop graphs is more subtle, because the latter vanish in dimensional regularization. However, the gauge-invariant coefficients I gg can also be obtained from the matching of partonic beam functions and PDFs with offshell external legs or an artificial gluon mass to regulate IR singularities. In this case also the massless virtual diagrams contribute to B (1) g/g and require zero-bin subtractions. The ξ-independence of the result again suggests that also the κ terms from diagrams fig. 3b-d must cancel upon zero-bin subtractions. Indeed, we find by explicit calculation that their total κ term before zero-bin subtraction exactly equals the total virtual zero-bin contribution. Hence, B (2,h) g/g is independent of κ. The crucial role of zero-bin subtractions for the gauge invariance of SCET matrix elements involving massive gauge bosons was already pointed out in ref. [63].
The zero-bin contributions from real-emission graphs vanish, see app. B.1. We thus conclude that dropping the κ terms from the start removes all zero-bin contributions. At the same time this eliminates all terms (∝ I 0 ) that cause the issues with the rapidity regulator discussed in sec. 3.1. We will therefore mostly exclude the κ terms in the following presentation of our beam function calculation. Instead we will treat them separately and explicitly demonstrate that they exactly cancel the non-vanishing virtual zero-bins in app. B.2.

Real emission diagrams
The relevant (one-and two-particle) real emission diagrams for the computation of B (a) where the splitting function P gq is given in eq. (C.2) and we defined for (later) conveniencê As noted above it is easy to see that the contributions ∝ κ from the diagrams for B (2,h) g/g with a massive quark bubble cancel. In particular, using the dispersion relation in eq. (3.12), the calculation of diagrams fig. 2g,h resembles the one of the real emission graphs contributing to B (2,h) q/g in ref. [31]. The graphs in fig. 2 with a coupling to a collinear Wilson line are rapidity divergent. (Diagram 2i vanishes upon integration.) Implementing the η rapidity regulator, according to eq. (3.5), for these diagrams amounts to an overall factor of 15 This factor regulates the 1/(1 − x) poles associated with rapidity divergences by translating them to 1/η poles via the expansion in terms of the plus distributions L n (y) ≡ θ(y) ln n y y The rapidity divergences of diagrams 2d-f (and their mirror graphs) cancel exactly. The rapidity divergence of diagram 2h instead cancels with its soft analog in fig. 6a within the cross section in eq. (2.3) at fixed order. The one-real-gluon cuts of the diagrams in fig. 2 (if present) give rise to terms singular in ⃗ p 2 T , i.e. proportional to δ (2) (⃗ p T ) or the plus distributions 16 via the expansion The cuts through two massive quark lines result in regular (non-singular) functions of ⃗ p 2 T . The reason is that the limit ⃗ p 2 T → 0 is effectively tied to the limit m 2 → ∞, where massive quarks in the final state are kinematically not allowed. The massive cuts therefore yield terms proportional to 1/m 2 rather than to 1/⃗ p 2 T (in d = 4 dimensions) and therefore need no regularization in terms of distributions.
The total zero-bin contribution to B (2,h) g/g associated with real emission diagrams is scaleless and vanishes similar to that from the massless diagrams. Details on the corresponding zero-bin calculation are presented in app. B.1.

Virtual diagrams
The one-loop contribution to B g/g due to a massive quark flavor is given by the diagram in fig. 1a. It corresponds to a wavefunction-type correction to the tree-level result B Here and in the following we use the shorthand notation The two-loop virtual diagrams contributing to B (2,h) g/g are displayed in fig. 3. As stated above and shown by explicit calculation in app. B.2 the total virtual zero-bin contribution to B (2,h) g/g exactly equals and thus cancels the terms proportional to the bookkeeping parameter κ in the unsubtracted expression for B (2,h) g/g . The calculation of the virtual diagrams is somewhat more involved than the one of the real emission graphs in fig. 2, because the loop integrations are not constrained by the measurement δ-functions in eq. (3.1). In fact, up to an overall δ (2) (⃗ p T ) the virtual contribution equals the one of the massive PDF matching coefficient M (2) gg 15 In the following we will suppress the parameter w and only restore it implicitly when required by the RRG formalism [40], i.e. for the derivation of the ν anomalous dimension of the beam function in eq. (4.17). 16 In refs. [40,50] the notation L T n (⃗ pT , µ) ≡ (−1) n computed in direct QCD [69]. Given the complete (one-and two-particle) real-emission contribution we could even extract the virtual part from the known M (2) gg by taking the small mass (m ≪ ⃗ p T ) limit of the TMD beam function. It is nevertheless instructive to calculate the virtual diagrams (and their zero-bins) in our SCET setup. In sec. 5.1 we will then use the small mass limit as a strong cross check of our explicit calculation. In the following we briefly discuss the evaluation and features of the different virtual diagrams in fig. 3: Diagram 3a has a massive quark triangle subgraph and is arguably the most difficult to compute. As there is no corresponding soft diagram with a massive triangle, all zero-bins are power-suppressed and do not contribute to B (2,h) g/g . For the same reason the diagram is rapidity-finite and does a priori not require any rapidity regulator. To simplify the involved two-loop Feynman integrals we want to use automated integration by parts (IBP) reduction to a minimal set of master integrals, a standard tool of modern multi-loop calculations. However, when eikonal (Wilson line) and massive propagators are present at the same time a naive implementation of IBP reduction often fails, see e.g. the discussion in ref. [67]. For the unregulated rapidity-finite two-loop diagram 3a the IBP program FIRE5 [70] for example leads to a rapidity-divergent expression due to unregulated rapidity-divergent master integrals. In the present case there is a pragmatic solution to this problem: We introduce an auxiliary rapidity regulator which ensures well-defined integrals in the result and at intermediate steps of the IBP reduction. After evaluation of the master integrals the dependence of the IBP reduced expression of the diagram on the rapidity regulator cancels out. For practical reasons we choose the δ-regulator for this calculation. The IBP reduction with FIRE5 [70] yields four master integrals, which we solved by direct integration in Feynman parameter representation as functions of m and δ/p − (in the limit δ → 0). The final result for diagram 3a is δ-independent, but depends on the gauge parameter ξ. The ξ-dependent terms exactly cancel the ones of diagram 3b. As a check we computed the ϵ poles of the diagram (for ξ = 0) also by direct loop integration without IBP reduction and rapidity regulation.
Diagram 3b is rapidity divergent and ξ-dependent. Implementing the η-regulator according to eq. (3.5) it entails rapidity-finite terms (∝ κ) that are however discontinuous in the η → 0 limit as discussed in sec. 3.1. These terms cancel exactly after zero-bin subtraction, see app. B.2, leaving a smooth result for the rapidity-finite part as η → 0 (which therefore can also be computed without regulator). Like for all relevant diagrams with a massive quark bubble a particular convenient and transparent way to calculate this diagram is via the dispersion relation as described in sec. 3.2. Note that for the virtual diagrams the Π(0, m 2 ) term in eq. (3.12) only gives rise to a vanishing contribution proportional to the corresponding virtual massless one-loop diagram, and can thus be dropped.
Diagram 3c is proportional to the bookkeeping parameter κ when deploying eq. (3.12). To derive the integrand one must carefully implement the Feynman rule for the triple gluon field strength vertex given in app. D of ref. [53]. One way to do this is to assign an auxiliary offshellness δ to the Wilson line propagators and perform the necessary color index contractions. The resulting Feynman integral then takes the form of I δ (1, 1) in eq. (3.8). Setting δ → 0 before the loop integration we end up with a result for diagram 3c that is proportional to I 0 (1, 1) in eq. (3.6) and corresponds to a (soft) Wilson line self-energy correction, cf. app. A. This explains why Diagram 3c exactly vanishes upon zero-bin subtraction, see app. B.2. It can thus be effectively omitted in the calculation of B (2,h) g/g as a whole (independent of whether or not an unnecessary rapidity regulator is applied).
Diagram 3e represents the full (QCD) two-loop wave function correction from 1PI massive quark vacuum polarization diagrams. One specific diagram included in this correction is diagram 3d. It contains the complete κ term of the two-loop wave function renormalization. The complete two-loop 1PI wavefunction contribution from a heavy flavor can be obtained from diagram 3e by inserting the known expression (with κ = 1) for the massive vacuum polarization function Π (2) (0, m 2 ) given in ref. [71]. 17 In order to verify our statement of sec. 3.2 that dropping the κ terms (and thus the zero-bins) also in the purely virtual contribution yields the correct result we must carefully extract the κ term due to wavefunction renormalization. To this end we evaluate where Π (2,κ) denotes the κ Term of the gluon self-energy subgraph in diagram 3d (suppressing color indices). Note that the gluon propagator connecting the subgraph with the O(g 0 ) vertex for the B σc n⊥ B c n⊥σ operator is (∝ g µν ⊥ and) in fact replaced by the derivative ∂/∂p 2 . This is necessary, because g µν ⊥ Π (2,κ) µν /p 2 (unlike the full g µν ⊥ Π (2) µν /p 2 ) is not regular in the on-shell limit p 2 → 0. The limit is equivalent to p + → 0 with p 2 = p − p + and may conveniently be taken after the derivative but before the loop integration in eq. (3.26). The contribution of diagram 3d is subtracted from that of diagram 3e (with κ = 1) to obtain the two-loop 1PI wavefunction correction to B (2,h) g/g without κ term (and zero-bins). Diagrams 3f and 3g represent the one-particle reducible contributions from wavefunction renormalization. They are straightforward to compute using eqs. (3.13) and (3.14). As usual the external leg correction diagrams 3d-g (and their left-right mirror graphs) must be multiplied with the fractional numbers given in the caption of fig. 3 to obtain the correct wavefunction renormalization contribution.
The soft zero-bin contributions of the virtual diagrams are necessary when the virtual κ terms are included. Their calculation is discussed in app. B.2. When expressing the bare result for B (2,h) g/g in terms of the MS renormalized coupling α s ≡ α {n l +1} s (µ) and on-shell renormalized mass m we have to add the counterterm contributions The massless one-loop result B (1,l) g/g is given in eq. (C.4). In the second term of eq. (3.27) the one-loop heavy-flavor contribution in eq. (3.24) is multiplied with the full one-loop coupling counterterm for n l + 1 active flavors, The one-loop mass counterterm in the on-shell scheme is

Two-loop TMD beam function results
We now have the full (real + virtual) two-loop results for the heavy flavor contributions B for n = 1, 2. The latter are fixed by RG consistency, which relates them to known expressions for hard and soft anomalous dimensions [31]. Our beam function calculation provides an explicit confirmation of these results, which serves as an important cross check.
The renormalized matching kernels I gk in eq. (2.5) are related to the bare partonic beam functions via 18 where Z Bg is the MS renormalization factor of the TMD gluon beam function operator and the sum over all massless partons i is understood. Throughout this paper we always express (any contributions to) B g/j , Z Bg , and I gj in terms of α s ≡ α {n l +1} s (µ) as indicated explicitly in eq. (4.1) by the superscript {n l + 1}. We stress that this also applies to the massless n-loop expressions B (n,l) g/j , Z (n,l) Bg , and I (n,l) gi , which arise from n l light quark flavors and gluons only. For compactness of notation we will often drop the {n l + 1} superscript in the following. The (ultra-) collinear PDFs live in n l -flavor QCD, where the heavy flavor has been integrated out, and are therefore naturally expressed in terms of α {n l } s (µ). In the MS scheme we have with the one-loop splitting functions P with L m ≡ ln(m 2 /µ 2 ). Expanding eq. (4.1) to O(α s ) and using With B (1,h) g/g in eq. (3.24) and B (1,h) g/q = 0 we thus obtain , where we have used eq.
In order to write eq. (4.8) in a compact form, we have introduced the auxiliary function The beam function anomalous dimensions are derived from the counterterm Z Bg as follows: At O(α s ) and O(α 2 s ) we thus obtain the heavy-flavor contributions In the derivation of anomalous dimensions within the RRG formalism using the η regulator it is in general important to retain higher order ϵ terms in the 1/η poles of the corresponding counterterms. In this particular case it is crucial to include the O(ϵ/η) term of Z (1,l) Bg , as given in eq. (C.7), in the formulas for the two-loop anomalous dimensions. Moreover, for eq. (4.15) it is necessary to restore the full µ dependence for finite ϵ in the 1/η terms of Z 1,l Bg and Z 2,h Bg , which are ∝ µ 2ϵ and µ 4ϵ , respectively. The results in eqs. (4.14) and (4.15) exactly equal the contributions of a single light flavor, see ref. [50]. The heavy-flavor contribution to the two-loop rapidity anomalous dimension in eq. (4.17) satisfies, according to eq. (2.12), the RG consistency relation 2γ ν,Sg is part of the anomalous dimension of the TMD soft function S gg . It is determined via Casimir rescaling from the result in ref. [31] and given in eq. (C. 16).
The results in this section represent the heavy-flavor contributions to the gluon TMD beam function and its anomalous dimensions through O(α 2 s ). The corresponding massless contributions I The massless part of the soft function S (n,l) gg and its anomalous dimensions are also found in ref. [50], while the soft heavy-flavor contributions are computed in ref. [31] and collected for convenience in app. C.5. We emphasize again that our explicit results imply that both massive and massless contributions must be evaluated with α s ≡ α {n l +1} s (µ). In the coefficients of the α s expansion of the massless contributions from ref. [50] we must however consistently replace n f → n l (e.g. in the expression for β 0 ).

Beam function asymptotics
In this section we study the limiting behavior of the beam function matching coefficients in eqs. (4.8) and (4.9) when m ≪ p T ≪ Q ("small mass limit") and p T ≪ m ≪ Q ("large mass limit"). Following ref. [31] we establish in this way the connection between the cross section in eq. (2.3) and the corresponding factorization formulas in the two limits (with p T ≡ |⃗ p T | ∼ q T ), where also the matching coefficients I gi themselves exhibit a factorized structure. The associated factorization ingredients and thus the asymptotic expressions for the I gi are either available in the literature or can be inferred by consistency from related factorization formulas for other processes. Verifying the expected limiting behavior therefore represents a valuable and strong check of our results for the I gi . In particular the virtual contributions, which are proportional to δ (2) (⃗ p T ) and thus unaffected by the small and large mass expansions, are directly cross-checked.

Small mass limit
In the small mass limit m ≪ q T (for Λ QCD ≪ m, q T ≪ Q) the factorized cross section of gluon-fusion color-singlet production takes the form (see ref. [31] for quark-initiated processes) × i∈{Q,Q,q,q,g} k∈{q,q,g} × j∈{Q,Q,q,q,g} l∈{q,q,g} Compared to eq. (2.3) the mass-dependent beam function matching coefficients are factorized in eq. (5.1) into the corresponding coefficients for n l + 1 massless quarks and the known PDF (flavor-threshold) matching factors M ij : Here the explicit superscript {n l +1}, which has been suppressed in eq. (5.1), indicates that the beam function matching coefficients must be evaluated using α s ≡ α {n l +1} s in order to avoid large logarithms ∼ ln n (m 2 /⃗ p 2 T ). The relevant M ij are collected up to O(α 2 s ) for convenience in app. C.3.
We now show explicitly that our mass-dependent two-loop results for the unpolarized coefficients I gi in sec. 4 indeed satisfy eq. (5.2). At one-loop we trivially have because the only contributing diagram is virtual, see fig. 1 a, and thus survives the expansion in eq. (5.2) as a whole. For the matching coefficient I (2,h) gq in eq. (4.9) we find in the small mass limit in agreement with eq. (5.2). The first term in the last line represents the contribution to I (2) gq due to a single massless quark flavor. The explicit expression is given in eq. (C.9). Note that for p T > 0 the m 2 /⃗ p 2 T expansion leading to eq. (5.4) is straightforward. However, to determine the correct distributional structure and in particular to fix the coefficient of δ (2) (⃗ p T ) on the right-hand side of eq. (5.4) we also had to expand the cumulant, i.e. the ⃗ p T -integral from 0 to an arbitrary ⃗ p cut T (with m ≪ |⃗ p cut T |), of the left-hand side. In the same manner we verify 20 that our result for the matching coefficient I (2,h) gg in eq. (4.8) is in accordance with eq. (5.2): Qg (m, z, µ) , (5.6) given in the literature. They can, however, be inferred from consistency relations between the different formulations of the factorization theorem in ref. [72] for deep inelastic scattering in the (x → 1) endpoint region and a corresponding alternative factorization theorem based on the approach of ref. [31]. 21 Concretely, one can show the relation where S g c is, up to Casimir rescaling, the csoft function in ref. [31] and M g (1 − z, m, µ) = lim z→1 M gg (z, m, µ) is the massive PDF matching coefficient in the threshold limit. The explicit expression for H g c that we extracted up to O(α 2 s ) from eq. (5.10) is given in eq. (C.13). We can now check our results for the heavy-flavor corrections to the gluon TMD beam function coefficients in the large mass limit against eq. (5.9). At one loop we consistently have The first term in the two-loop matching coefficient I in eq. (4.9) vanishes in the decoupling limit and we obtain The last term in eq. (5.12) arises from the flavor threshold matching relation eq.
The asymptotic behavior of the beam function coefficients in eqs. (5.12) and (5.13) agrees with eq. (5.9) and thus confirms our results.

Numerical effect of bottom mass corrections
As a first step to quantify the effect of the quark mass corrections obtained in sec. 4 on physical observables, in particular the Higgs transverse momentum spectrum, we assess here their numerical size at fixed order in α s . A full-fledged analysis including resummation and the appropriate renormalization scale variations as well as the matching to the massless full QCD fixed-order result relevant at large transverse momentum is left for future work.
To remove the dependence on the rapidity renormalization scale ν we consider the symmetrized combinatioñ with S (2,h) gg as given in eq. (C.14). In fig. 4 we show plots of (2πp T times)B as a function of p T and p cut T , respectively. For the plots we set µ = m = m b = 4.8 GeV and x = ω/E cm with ω = Q = m H = 125 GeV, E cm = 13 TeV, and we used MMHT2014 NNLO PDFs [73]. The quark mass corrections can be expressed as an infinite series of the subleading terms ∼ (m/p T ) 2n in the small mass expansion with n ≥ 1. Note that fixing µ does not affect the (m/p T ) 2n corrections we want to visualize here: The difference between the result with the full mass dependence (red curve) and its small mass limit (blue dashed curve) is (unlike the individual curves) µ-independent, because the beam and soft function (and equivalently the hard function) µ anomalous dimension is mass-independent. Note also that this difference is of O(α 2 s ), i.e. the full and the small mass results are equal at O(α s ). We observe that for p T ≫ m the deviation between the full result and the small mass limit performed in sec. 5.1 is indeed small, while for p T ∼ 10 GeV ∼ 2m b the deviations are of O(100%) and the small mass result does not provide a sensible approximation of the O(α 2 s T F ) contributions. In the large mass limit and for µ = m the correctionB (2,h) g is proportional to δ (2) (⃗ p T ) as can be verified from the results in sec. 5.2 and ref. [31]. In the plots of fig. 4 the 22 Recall that q andq stand for the n l massless quark and antiquark flavors, respectively. large mass limit is therefore illustrated by (dotted green) horizontal lines at zero (left panel) and a nonzero value (right panel), which are touched by the curves of the full result at p T = 0 and p cut T = 0, respectively. Finally, we assess the NNLO (=NLO 1 ) quark mass corrections ∼ (m/q T ) 2n to the Higgs q T distribution in top-induced gluon fusion at the LHC. In fig. 5 we show their effect relative to the full NLO (= LO 1 ) result, i.e. relative to the leading order spectrum for q T > 0. Concretely, we plot the cross section ratio 23 as a function of q T > 0 with the same input as for fig. 4. The newly computed quark mass power corrections are nonsingular in the m → 0 limit and composed of terms ∝ (m/q T ) 2n with n ∈ N, which can also contain positive powers of ln(m/q T ). The total mass-nonsingular contribution to eq. (6.4) is depicted in fig. 5 by the dotted violet curve. It corresponds to the difference of the full result (solid red curve) and the small mass limit (blue dashed curve), i.e. the mass-singular contribution, which includes massive PDF matching factors, as described in sec. 5.1. Similar to the plots in fig. 4, this difference (the dotted violet curve) is µ-independent, while the full and small mass results for dσ (2,h) /dσ (1) individually depend on µ, as dσ (1,h) is nonzero for µ ̸ = m.
Although the obligatory resummation (with mass-dependent RRG evolution kernels) in the peak region (i.e. for q T ∼ 10 GeV) may quantitatively change the result, fig. 5 should nevertheless give a reasonable estimate of the potential size of the bottom mass corrections to the Higgs q T distribution in gluon fusion mediated by a top loop. As expected, the corrections (represented by the dotted violet curve) become negligibly small for q T ≫ 10 GeV, i.e away from the peak region. Around the peak they amount to ∼ 1 − 2% and increase for smaller q T . 23 Note that in accordance with eq. (2.3) the qT -independent hard function factors do not contribute for µ = m and qT > 0 at the order of interest.
Despite their apparent small size, these corrections will likely matter for precision analyses of the Higgs q T spectrum at N 3 LL ′ accuracy (and beyond), which already has reached the few percent level [6].

Conclusion
We have calculated the TMD gluon beam functions at NNLO in SCET with one massive and n l massless quark flavors. Our results for the massive quark contributions to the renormalized beam function matching kernels are presented in sec. 4. The relevant two-loop diagrams associated with collinear real emissions are shown in fig. 1b and fig. 2. Their calculation is rather straightforward. The purely virtual two-loop diagrams in fig. 3, however, require a careful subtraction of non-trivial zero-bin contributions and involve terms that, upon rapidity regularization, exhibit a discontinuous behavior in the limit of vanishing rapidity regulator even though they are rapidity-finite. Both types of contributions arise from two-loop diagrams with a massive quark bubble subdiagram on a virtual gluon line. More precisely, they are generated by the part proportional to p µ p ν of the virtual gluon propagator dressed with a massive quark bubble, where p µ is the off-shell four-momentum flowing through the gluon line, see sec. 3.2. The corresponding parts of the two-loop diagrams are referred to as κ terms in this work. We explicitly show that the zero-bin subtraction completely removes all κ terms from the final beam function result. In other words, one obtains the correct NNLO expressions by omitting the κ terms in all diagrams with a massive quark bubble and thus also the non-vanishing zero-bin contributions.
An important application of our beam function results is the computation of bottom mass effects on the gluon-fusion production of Higgs bosons with small transverse momenta (q T ≪ m H ). We have derived the full (m b /q T )-dependence of the Higgs q T distribution at leading order in the QCD coupling, i.e. at relative O(α s ), and at leading power in 1/m H , i.e. at O(y 0 b ) with y b the bottom Yukawa coupling. We have also confirmed the anomalous dimensions relevant for the resummation of logarithms ∼ ln(m 2 b /m 2 H ) ∼ ln(q 2 T /m 2 H ) at NNLL ′ . For N 3 LL resummation only the quark mass corrections to the three-loop rapidity anomalous dimension is yet unknown. Apart from the process-dependent hard function, i.e. an overall factor, our results directly carry over to the transverse momentum distribution of any other color singlet final state produced by gluon fusion. We have performed a first numerical analysis at fixed order and found a few-percent level effect of the new bottom mass corrections on the Higgs q T spectrum in the peak region, where q T ∼ m b . A more sophisticated analysis including the resummation of large logarithms based on the factorization approach of ref. [31] for the different hierarchies between m b , m H , and q T ≪ m H is left for future work.

A Review of soft function calculation
The two-loop heavy-flavor contribution S (2,h) gg to the TMD soft function S gg in the factorization formula eq. (2.3) was derived in ref. [31]. 24 The relevant soft function diagrams are displayed in fig. 6. The calculation was performed using the dispersion relation in eq. (3.12) with the κ terms set to zero as justified by a gauge invariance argument following ref. [63]. In this appendix we explicitly demonstrate the cancellation of κ terms among the diagrams in fig. 6. This requires a careful evaluation of diagram 6d, which is closely related to the zero-bin contribution of the collinear diagram 3d, see app. B.2.
By means of eq. (3.12) the two-loop calculation of S (2,h) gg is split into two one-loop calculations with a massive and a massless gluon exchanged between the soft Wilson lines, respectively. The massless calculation yields the known result [40,50] with an overall factor ∝ Π (1) (0, m 2 ). Here we focus on the one-loop diagrams corresponding to the graphs in fig. 6 where the gluon line with a massive quark bubble is replaced by a gluon with mass M over which we eventually integrate according to eq. (3.12). For the real corrections from diagrams 6a and 6b we have to cut the massive gluon propagator, i.e. effectively replace it with a δ-function that puts the gluon momentum on the mass shell (with positive energy). The corresponding one-loop integrands read up to a common overall factor Applying the η-regulator [40] both integrands in eq. (A.1) are multiplied with the same factor ν η |k − − k + | −η . We see that the κ terms of the two real-emission integrands cancel each other. The (normalized) integrand of the virtual diagram 6c is given by (suppressing the i0- 24 The calculation in ref. [31] was done for incoming quark and antiquark, i.e. with Wilson lines in the fundamental color representation. Their two-loop result can however be directly translated to S (2,h) gg due to Casimir scaling by replacing the overall color factor CF TF → CATF . prescription) .
Again the η-regulator adds a factor ν η |k − − k + | −η . Note, however, that the κ term is rapidityfinite and the result for this term is independent of whether the η → 0 limit is taken before or after the integration (in contrast to the collinear integral in eq. (3.7)). This κ term is exactly canceled by the soft wave function renormalization represented by diagram 6d. The integrand [6d], similar to that of diagram 3d, requires an offshellness δ of the external (Wilson) line with the self energy insertion, and the η-regulator must not be applied. In SCET, a consistent implementation of the offshellness in soft diagrams requires the parametric scaling δ ∼ m 2 /Q ≪ m, because δ is related to the offshellness p 2 ∼ Qδ ∼ m 2 of the underlying full QCD (and corresponding collinear) diagrams, like in eq. (3.26). To comply with the EFT power counting we therefore must expand the integrand of diagram 6d in δ/k ∼ δ/M ∼ m/Q, see also ref. [63]. With the same normalization as for eq. (A.2) and taking into account the factor 1/2 for wave function renormalization we have (with k − ↔ k + for two of the mirror diagrams) The 1/δ pole in eq. (A.3) vanishes upon integration over the soft loop momentum k µ as it is antisymmetric under k µ ↔ −k µ . 25 So, (before integration over M ) the total η-regulated contribution of diagrams 6c, 6d, and their mirror graphs is proportional to and thus κ-independent. The final MS renormalized result for S (2,h) gg is given in eq. (C.14).

B Soft zero-bin contributions
In this appendix we give details on the calculation of the soft zero-bin contributions. These must be subtracted from the sum of the diagrams in figs. 2 and 3 in order to avoid double counting soft contributions already contained in the diagrams of fig. 6. The zero-bin contributions are obtained by expanding the integrand (including the measurement δ-functions) of the corresponding collinear diagrams assuming a loop-momentum scaling such that the momentum passing through at least one of the propagators is soft, while the momenta of the other propagators may be soft or collinear [64]. In this way, for each two-loop diagram in figs. 2 and 3 different zero-bins associated with soft-collinear and soft-soft loop-momentum regions arise. Except for some zero-bins from diagrams with a massive quark bubble, however, we find that all of these are power-suppressed (w.r.t. 1/Q) and/or individually scaleless and thus do not contribute to the (leading-power) beam function kernel I gg . In particular, there is no non-trivial zero-bin contribution from diagrams with massive quark triangle or box subgraphs.

B.2 Virtual zero-bin contributions
The relevant (unregulated) zero-bin integrands of the virtual diagrams in fig. 3 are where ℓ µ ∼ M ∼ m and we suppressed a common factor ∝ δ (2) (⃗ p T ) δ(1−x) Im Π (1) (p 2 , m 2 ) /M (as well as the i0 prescription). Implementing the η-regulator adds a factor ν η |ℓ − | −η to the integrand [3b] 0-bin , so that the integrated zero-bin contribution of diagram 3b vanishes (for η ̸ = 0). As discussed in sec. 3.1, the κ term of [3b] 0-bin exactly cancels the term in the integrand of diagram 3b that gives rise to the discontinuous η → 0 limit due to the integral in eq. (3.7). The rapidity-finite part (including the κ term) of the zero-bin subtracted diagram 3b therefore does not need to be rapidity-regulated. The zero-bin integrand [3c] 0-bin exactly equals that of the unsubtracted diagram 3c. The zero-bin subtracted diagram 3c therefore vanishes (regardless of rapidity regularization).
The zero-bin integrand [3e] 0-bin is derived from the unintegrated eq. (3.26) respecting the scaling p 2 ∼ m 2 and thus p + ∼ m 2 /Q. As a consequence of SCET power counting we must take the p + → 0 limit at the level of the zero-bin integrand. This is consistent with the evaluation of the soft diagram 6d in app. A. The factor 1/2 in [3e] 0-bin is due to wavefunction renormalization. The zero-bin contribution of diagram 3e effectively cancels the soft wave function contribution represented by diagram 6d once the corrections from all four external legs (of the soft function and the two partonic beam functions, respectively) are combined in the factorized cross section in eq. (2.3). This is expected because the unsubtracted collinear wavefunction renormalization exactly equals that of full QCD [33].
After zero-bin subtraction the κ terms of diagrams 3b and 3d cancel each other exactly. This resembles the cancellation of the κ dependence from diagrams 6c and 6d within the virtual contribution to the soft function as shown in app. A.

C.1 Splitting functions
The leading-order (one-loop) PDF anomalous dimensions γ ij /π are given by P (0) q i q j (z) = C F θ(z) δ ij P qq (z) , P (0) q i g (z) = P (0) q i g (z) = T F θ(z)P qg (z) , P (0) gg (z) = C A θ(z)P gg (z) + with q i explicitly denoting here the different massless quark flavors, and the usual one-loop (LO) quark and gluon splitting functions

C.3 Massive PDF matching factors
The massive PDF matching coefficients were calculated up to two loops in ref. [69].

C.4 Hard massive threshold correction
The matching correction H g c due to collinear mass modes arising in the limit Q ≫ m ≫ q T can be inferred from the literature via eq. (5.12). Up to two-loop order we find H g c m, µ,