Resummation and Matching of $b$-quark Mass Effects in $b\bar{b}H$ Production

We use a systematic effective field theory setup to derive the $b\bar{b}H$ production cross section. Our result combines the merits of both fixed 4-flavor and 5-flavor schemes. It contains the full 4-flavor result, including the exact dependence on the $b$-quark mass, and improves it with a resummation of collinear logarithms of $m_b/m_H$. In the massless limit, it corresponds to a reorganized 5-flavor result. While we focus on $b\bar{b}H$ production, our method applies to generic heavy-quark initiated processes at hadron colliders. Our setup resembles the variable flavor number schemes known from heavy-flavor production in deep-inelastic scattering, but also differs in some key aspects. Most importantly, the effective $b$-quark PDF appears as part of the perturbative expansion of the final result where it effectively counts as an $O(\alpha_s)$ object. The transition between the fixed-order (4-flavor) and resummation (5-flavor) regimes is governed by the low matching scale at which the $b$-quark is integrated out. Varying this scale provides a systematic way to assess the perturbative uncertainties associated with the resummation and matching procedure and reduces by going to higher orders. We discuss the practical implementation and present numerical results for the $b\bar{b}H$ production cross section at NLO+NLL. We also provide a comparison to the corresponding predictions in the fixed 4-flavor and 5-flavor results and the Santander matching prescription. Compared to the latter, we find a slightly reduced uncertainty and a larger central value, with its central value lying at the lower edge of our uncertainty band.


Introduction
The formulation of reliable predictions for heavy-quark initiated processes has been the subject of much study over many years, in particular for the determination of parton distribution functions (PDFs) in the context of deep-inelastic scattering (DIS) [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. At the LHC, important examples of heavy-quark initiated processes are Higgs or vector-boson production in association with heavy quarks. In this paper, we are interested in Higgs production in association with b quarks, i.e., the inclusive bbH-induced cross section. In a typical hard-scattering process with protons in the initial state, there are at least two parametrically separate scales. First, the hard scale µ H ∼ Q, where Q denotes the physical quantity that determines the momentum transfer in the hard interaction, e.g., Q = −q 2 in DIS or Q = m H for the case of Higgs production that we will be interested in. Second, the low scale µ Λ ∼ Λ QCD , which separates the perturbative and nonperturbative regimes and is typically taken to be of order the proton mass, µ Λ ∼ 1 GeV. In the limit µ Λ Q we can apply the standard QCD factorization theorem [18][19][20] to compute the hadronic cross section in terms of the partonic cross section convolved with PDFs.
For heavy-quark initiated processes, the mass m of the heavy quark introduces another physical scale. Depending on its value, we can distinguish two parametrically different cases, shown in figure 1: 1 (a) m ∼ Q: There is a single parametric scale µ H ∼ m ∼ Q in addition to µ Λ .
(b) m Q: There are two parametric scales µ H ∼ Q and µ m ∼ m in addition to µ Λ .
When working in the limit m ∼ Q, the heavy quark never appears in the initial state of the hard partonic process. Instead, it is produced as part of the hard interaction at µ H by an incoming gluon splitting into a pair of heavy quarks. The partonic calculation contains the exact dependence on m, including the correct m-dependent phase space. The gluon splitting into a heavy-quark pair contains a collinear singularity, which is regulated by m, and as a result produces logarithms ln(m/Q). For m ∼ Q, these collinear logarithms are counted as small and are included at fixed order in the α s expansion. When working in the limit m Q, the heavy quark explicitly appears in the initial state of the hard partonic process, and the collinear logarithms are resummed to all orders in α s into an effective heavy-quark PDF. The quark mass m only appears in the boundary condition of the PDF's DGLAP evolution, which starts at the scale µ m ∼ m. The hard process itself is computed in the m → 0 limit. That is, finite-mass effects of O(m/Q), including the exact phase space of the gluon splitting into a massive quark pair, are power corrections and are neglected.
Predictions obtained in strictly one of the above two limits are usually referred to as obtained in a fixed-flavor number scheme. Which of these limits is more appropriate in Figure 1. The two parametric scale hierarchies for inclusive cross sections for heavy-quark initiated processes.
practice depends on the process and the numerical size of the corrections. For b-initiated processes at hadron colliders, the relative importance of ln(m b /Q) corrections has been discussed for example in [21][22][23].
To obtain the best possible theoretical predictions, it is often desirable to have a complete description that incorporates the results from both limits. In this way, the final result is valid in each limit as well as in the transition region in between, and hence one can be agnostic about which parametric regime is the more appropriate one.
For bbH, predictions exist in the 4-flavor scheme (4FS) [24,25], which works in the limit m b ∼ m H , and in the 5-flavor scheme (5FS) [26][27][28][29], which works in the limit m b m H . Currently, both predictions are combined using the pragmatic "Santander Matching" prescription [30], which is a weighted average of the 4FS and 5FS predictions, where the relative weighting depends on the numerical size of ln(m b /m H ).
There are various methods available in the literature, referred to as variable-flavor number schemes (VFNS), which aim to combine the virtues of both limits in a more systematic fashion. That is, they include the full m dependence in the limit m ∼ Q and the resummation of collinear logarithms ln(m/Q) in the limit m Q. There are a number of such schemes available, namely the ACOT scheme [1,2] (and its simplified variants S-ACOT [6], S-ACOT-χ [7,13], and the more recent m-ACOT [31] for hadron-hadron collisions), the TR scheme [3,8], and the FONLL scheme [11,32]. The differences between the schemes essentially amount to how the two limits are combined.
Effective field theories (EFTs) are the standard tool to describe processes with parametrically separated scales, allowing to systematically resum the logarithms of ratios of these scales. In this paper, we discuss the EFT formulation of heavy-quark initiated processes for the case of inclusive cross sections. All the basic ingredients are actually well known in this case. Nevertheless, we find it worthwhile to discuss the EFT formulation in detail, as it provides a conceptually clear field-theoretic derivation, including the transition between the two parametric regimes and a way to assess the associated theoretical uncertainties. This setup can also be extended to more differential cross sections, which we leave for future work. A similar setup has also been used to incorporate quark-mass effects for final-state jets in Refs. [33][34][35].
Our final result for DIS resembles the aforementioned schemes in several ways, but also differs in some key aspects. Most importantly, the b-quark PDF is not treated as an external O(1) quantity. Rather, it contributes as part of the perturbative series of the final result, where it effectively counts as an O(α s ) object. (In this work, we follow the assumption made in all available PDF sets that there is no intrinsic bottom in the proton, such that the effective bottom-quark PDF is generated purely perturbatively.) The application of our method to hadron-hadron collisions is completely straightforward. Our final result for bbH encompasses the merits of both 4F and 5F schemes. It contains the full 4FS result at NLO, including the exact m b dependence and phase space. In addition, it improves the 4FS result with the all-order resummation of collinear logarithms up to NLL order. In the m b → 0 limit, our result corresponds to a reorganized 5FS result, where the perturbative series is expanded to NLO with the b-quark PDF counted as O(α s ).
In the next section, we discuss the general setup in detail, focussing on DIS to be specific. In section 3 we briefly discuss the similarities and differences with respect to other heavy-flavor schemes in the literature. Then in section 4, we apply this framework to bbH production. We discuss in detail the perturbative uncertainties and present our final numerical results at LO+LL and NLO+NLL. We also compare to the predictions in the 4F and 5F schemes using a consistent set of inputs. We conclude in section 5.

EFT formulation of heavy-quark initiated processes
In this section, we discuss the EFT formulation in detail. For simplicity and to be specific we frame the discussion in the context of heavy-quark production in DIS, where we have to deal with only one strongly-interacting initial state. In this case we associate Q = −q 2 . The extension to hadron-hadron collisions is straightforward and will be discussed in section 4. For definiteness we consider the heavy quark to be the b quark, 2 and treat the four lighter quarks as massless. We take Q < m t , so we can essentially ignore the top quark (i.e., we either integrate it out at the scale µ H ∼ Q or it has already been integrated out at a higher scale).
In section 2.1, we review the case µ Λ m b ∼ Q, corresponding to figure 1(a), where a single matching step at the hard scale µ H ∼ Q ∼ m b is required. We will refer to this as the fixed-order region or limit. This also serves to introduce our notation and language. In section 2.2, we discuss the case µ Λ m b Q, corresponding to figure 1(b), where two separate matching steps, at µ H ∼ Q and µ m ∼ m b , are performed. We will refer to this as the resummation region or limit. In section 2.3, we discuss the appropriate perturbative counting in our result, and in section 2.4, we combine the results in both limits to yield our final predictions valid in both limits and anywhere in between.
Throughout this paper, we use roman indices i, j, k to denote the light flavors, i.e., the four light quarks and the gluon. We also use the convention that any repeated indices are implicitly summed over (also a repeated index b implies a sum over b andb). For clarity, we will focus on the dependence on the relevant physical and renormalization scales, but suppress all other kinematic dependences. In particular, we will not write out the dependence on the momentum fractions and the Mellin-type convolutions in them. We will denote the number n f of light active flavors as superscripts for quantities where the distinction is relevant, e.g., α [4] s vs. α [5] s .

m b ∼ Q: Fixed order
In this case, shown in figure 1(a), the b-quark mass is treated parametrically as of the same size as Q. At the scale µ H ∼ Q ∼ m b , all degrees of freedom with virtualities ∼ Q 2 ∼ m 2 b , including the heavy b quark, are integrated out. We match full QCD onto a theory of collinear gluons and collinear light quarks with typical virtuality Λ 2 QCD . 3 This matching step is precisely equivalent to the standard operator product expansion (OPE) in DIS [36][37][38][39][40][41][42][43][44][45], which we briefly review now.
We define the DIS operator O DIS (Q, m b ), whose proton matrix element determines the DIS cross section (or equivalently the hadronic tensor or DIS structure functions), At the scale µ H , it is matched onto a sum of nonlocal PDF operators where a sum over light quarks and gluons i = u, d, s, c, g is understood, and "⊗" denotes the Mellin-type convolutions in the momentum fractions. The Wilson coefficients D i (Q, m b , µ H ) are also called coefficient functions. The O [4] i (µ) are the standard MSrenormalized quark and gluon PDF operators [45] 4 , whose proton matrix elements define the nonperturbative PDFs, f [4] i (µ) = p O [4] i (µ) p .
Since the b quark is being integrated out and not present in the theory below µ H , there is also no O b operator and no D b coefficient on the right-hand side of eq. (2.2). As indicated, the right-hand side of eq. (2.2) is the leading term in an expansion in Λ 2 QCD /Q 2 , where the 3 In SCET, this is the purely n-collinear quark and gluon sector, which is equivalent to a boosted version of QCD, where n µ = (1, n) and n is the direction of the incoming proton. In lightcone coordinates, the momentum of the collinear modes scales as pc ∼ (Q, Λ 2 QCD /Q, ΛQCD). In principle, there could also be soft modes with momentum scaling ps ∼ (ΛQCD, ΛQCD, ΛQCD), and also Glauber modes. Since their contributions cancel in the inclusive cross section [20], they are not needed here. 4 For corresponding operator definitions in SCET and a discussion of their equivalence see e.g. refs. [46][47][48].
low scale Λ 2 QCD is set by the external proton state we are eventually interested in. For ease of notation, we will not indicate these power corrections in the rest of this section.
In the above, α s (µ) ≡ α [4] s (µ) and O [4] i (µ) are renormalized with n f = 4 active quark flavors. That is, we use MS with dimensional regularization with respect to the four light quark flavors, while b-quark loops are renormalized in the decoupling scheme, such that the b-quark decouples from the theory below µ H (see Appendix A).
Since O DIS determines the full-theory cross section, it does not have an explicit dependence on µ H , i.e., it does not receive additional operator renormalization. It only has an implicit dependence on µ H through the renormalization of α s (µ H ), which cancels order by order in perturbation theory. On the other hand, the coefficients D i are explicitly µ dependent, and their µ dependence cancels against the explicit µ dependence of the operators The full dependence on the physical scales Q and m b , which are treated as hard scales, resides in the Wilson coefficients D i (Q, m b , µ H ). The coefficients D i at some scale µ contain logarithms ln(µ/Q) ∼ ln(µ/m b ). Therefore, they are computed by a perturbative matching calculation (see below) at the hard scale µ H ∼ Q ∼ m b , where they contain no large logarithms.
The PDFs f i (µ) at some scale µ contain logarithms ln(µ/Λ QCD ). Hence, the input PDFs that are determined from the experimental data are defined at a low scale µ Λ Λ QCD , which should still be large enough for perturbation theory to be valid. All contributions from lower scales, including the nonperturbative regime, are absorbed into the input PDFs f i (µ Λ ). The renormalization of the PDF operators leads to their renormalization group equation (RGE) where γ ij are the PDF anomalous dimensions, which are given in terms of the standard Altarelli-Parisi splitting functions [38,40,43]. The solution of this RGE yields the standard DGLAP evolution [43,49,50] relating the operators (and PDFs) at a scale µ 0 to the operators (and PDFs) at the scale µ, O [4] i (µ) = U [4] ij (µ, µ 0 ) ⊗ O [4] j (µ 0 ) . (2.5) As denoted, the anomalous dimensions and evolution factors involve n f = 4 light quark flavors. By definition, the coefficients and operators in eq. (2.2) must be evaluated at the same scale, which means the operators on the right-hand side give PDFs at µ H containing large logarithms ln(µ H /µ Λ ). These logarithms are resummed by using eq. (2.5) to evolve the PDFs from the low scale µ Λ up to µ H , f [4] i (µ H ) = U [4] ij (µ H , µ Λ ) ⊗ f [4] j (µ Λ ) . (2.6) Equivalently, we can perform the resummation for the Wilson coefficients. The coefficients and operators obey inverse RGEs, since their scale dependences must cancel each Figure 2. Schematic leading-order matching for heavy-quark production in DIS for m b ∼ Q. Figure 3. Schematic NLO matching for heavy-quark production in DIS for m b ∼ Q.
other. After performing the matching at the scale µ H , the coefficients are evolved from µ H down to µ Λ , The evolution factor is precisely the same as in eq. (2.6). Evolving the coefficients down corresponds to successively integrating out virtualities between µ 2 H and µ 2 Λ . After evolving down to µ Λ , we can take the proton matrix element to obtain the final DIS cross section (2.8) The full cross section on the left-hand side contains large logarithms ln(Q/Λ QCD ) and ln(m b /Λ QCD ), which on the right-hand side are factorized into logarithms ln(Q/µ H ) and ln(m b /µ H ), which are considered small and reside in the coefficients, large logarithms ln(µ H /µ Λ ), which are resummed into the evolution factor, and logarithms ln(µ Λ /Λ QCD ), which are absorbed into the PDFs. This result for the cross section is precisely the 4FS result. Since the m b dependence is included at fixed order in eq. (2.8), we will refer to it as the fixed-order ("FO") result.
The Wilson coefficients D i (Q, m b , µ H ) are determined in perturbation theory by matching the matrix elements of both sides of eq. (2.2) between the same partonic external states j, The left-hand side corresponds to the full-theory matrix element. The partonic matrix element of the MS-renormalized PDF operators on the right-hand side are the partonic PDFs, f [4] i/j (µ H ) = j O [4] i (µ H ) j ≡ Γ [4] ij . (2.10) They are equivalent to the collinear MS subtractions often denoted as Γ [4] ij . The matching in eq. (2.9) is performed order by order in α s (µ H ) ≡ α [4] s (µ H ), for which we expand each of the pieces as (2.11) The leading-order matching is shown schematically in figure 2. At the lowest order, only the gluon external state contributes. Light quarks in the external state first contribute at NLO. The b-quark does not appear as external state in the matching calculation, since it cannot appear anymore on the right-hand side. Writing out the dependence on the momentum fraction z explicitly, the partonic PDFs at LO are simply given by g is directly given by the LO diagram on the left of figure 2, The schematic matching at NLO is shown in figure 3. There are virtual and real emission corrections to the gluon channel as well as the contribution with light (anti)quarks in the external state. Using pure dimensional regularization to regulate both the UV and IR, the MS-renormalized partonic PDFs at NLO are f [4](1) where β 0 (n f ) = (11C A −4T F n f )/3 (with n f = 4 here) and the one-loop (LO) gluon splitting functions are Together with the LO coefficients from eq. (2.13), the NLO matching coefficients are The 1/ terms in the f (1) i/j in eq. (2.14) are collinear IR divergences. They precisely cancel between the two terms on the right-hand side, such that NLO Wilson coefficients are free from IR divergences. While the same IR regulator must be used in the full and effective theories, the Wilson coefficients are independent of the specific IR regulator. On the other hand, the coefficients do explicitly depend on the UV renormalization scheme of the PDF operators, which is where the standard MS scheme is used. In other words, the fact that eq. (2.14) are a pure pole contribution is just an artifact from using pure dimensional regularization for both UV and IR divergences. 5 Using any other IR regulator, e.g., putting the external states off shell, the f (1) g/j would look different, but the final results for the D (2) i would be exactly the same once the IR regulator is taken to zero.

m b Q: Resummation
In this case, shown in figure 1(b), there is a parametric hierarchy between the b-quark mass m b and Q. The calculation now proceeds via a two-step matching. First, at the scale µ H ∼ Q, all degrees of freedom with virtualities ∼ Q 2 are integrated out and full QCD is matched onto a theory of collinear gluons, collinear light quarks, and in addition collinear massive b-quarks, all with typical virtualities p 2 c ∼ m 2 b . 6 Next, we evolve from µ H down to the intermediate scale µ m ∼ m b . At µ m , all degrees of freedom with virtualities ∼ m 2 b are integrated out, including the massive b-quark, and the theory is matched onto the same theory as in the previous section 2.1 of collinear gluons and collinear light quarks with typical virtuality Λ 2 QCD . The matching at the scale µ H proceeds as before, except that above µ m the bottom quark is still a dynamical degree of freedom. Analogous to eq. (2.2), the DIS operator is matched onto a sum of nonlocal PDF operators, which now includes the bottom-PDF Here, we use the notation C i,b for the Wilson coefficients to distinguish them from the D i coefficients in the previous subsection. The PDF operators, O [5] i and O [5] b , have the same structure as those in eqs. (2.2) and (2.3). The essential difference is that they are now renormalized with n f = 5 active flavors. That is, b-quark loops are now renormalized using MS with dimensional regularization.
Eq. (2.17) corresponds again to the standard OPE in DIS. It is important to note, however, that the expansion performed in eq. (2.17) is by construction an expansion in p 2 /Q 2 , where p 2 is the typical virtuality of the external states in the theory below µ H . 5 Technically, all loop corrections to the bare PDF matrix elements are scaleless and vanish, which means the UV and IR divergences are precisely equal with opposite sign and cancel each other. Adding the UV counterterms then leaves the IR divergences. In this case, the µ dependence in f i/j (µ) is purely through αs(µ) and the coefficients f i/j have no explicit µ dependence as written in eq. (2.11). 6 In SCET this would be a theory containing massive collinear fermions [51]. The collinear modes have momentum scaling pc ∼ (Q, m 2 b /Q, m b ). The corresponding soft modes with momentum scaling ps ∼ (m b , m b , m b ) are again not needed since they cancel. This implies that the production of secondary b-quarks can only arise from the splitting of collinear gluons.
Compared to the fixed-order case in eq. (2.2), where we had p 2 ∼ Λ 2 QCD , we now have p 2 ∼ m 2 b , which not only includes b-quarks but also collinear gluons of that virtuality. Thus, as indicated in eq. (2.17), it is always an expansion in m 2 b /Q 2 . In particular, matrix elements with external b quarks (e.g. in the matching calculation below) are expanded in the m b → 0 limit. The coefficients C i,b (Q, µ H ) contain the full dependence on the physical scale Q but are independent of the low scale m b . On the other hand, the operators do contain an implicit m b dependence since they involve massive b-quark fields.
In principle, one is free to reabsorb some of the neglected O(m 2 b /Q 2 ) corrections in eq. (2.17) into the coefficients. This corresponds to including some subleading power (subleading twist) corrections in the leading-power (leading-twist) result and letting the leading-power resummation act on them. However, we stress that to correctly include the subleading power corrections in the resummation requires extending the factorization in eq. (2.17) to the subleading order. As mentioned before, the power corrections in m 2 b /Q 2 can be important in practice and should be added back such that in the fixed-order limit µ m → µ H we recover the fixed-order result of the previous subsection. This is discussed in detail in section 2.4. In the rest of this section, we do not indicate the power corrections for ease of notation.
After the matching at µ H in eq. (2.17), we want to evolve the theory from µ H down to µ m . The renormalization of the PDF operators again gives rise to their RGE, the solution of which is given by DGLAP evolution, relating the operators at different scales, O [5] i (µ) = U [5] ij (µ, µ 0 ) ⊗ O [5] j (µ 0 ) + U [5] ib (µ, µ 0 ) ⊗ O [5] b (µ 0 ) , O [5] b (µ) = U [5] bj (µ, µ 0 ) ⊗ O [5] j (µ 0 ) + U [5] bb (µ, µ 0 ) ⊗ O [5] b (µ 0 ) . (2.18) The difference to eq. (2.5) is that now n f = 5 and O [5] b contributes to the evolution, which we have written out explicitly. Taking the proton matrix elements on both sides yields the corresponding evolution of the PDFs from µ 0 to µ in the theory above µ m . Equivalently, we can use eq. (2.18) to evolve the Wilson coefficients from µ H down to µ m , Next, at the scale µ m , the operators O [5] i (µ m ) and O [5] b (µ m ) are matched onto the set of operators O [4] i (µ m ), which are precisely the ones appearing in eq. (2.2) and do not include a b-quark operator, By integrating out the b quark, the m b dependence implicit in the O [5] j,b (µ m ) is now fully contained in the matching coefficients M jk (m b , µ m ) and M bk (m b , µ m ). In particular, the O b operator does not exist in the theory below µ m and its effects are moved into the M bj coefficient. In addition, secondary b-quark loops are integrated out, which corresponds to switching the UV renormalization scheme for b quarks at the scale µ m from MS to the decoupling scheme, and the M ij and M bj contain the associated matching (threshold) corrections.
The remaining steps now proceed as in the previous subsection. The operators (or PDFs) at µ m still contain logarithms ln(µ m /Λ QCD ), which are resummed by using eqs. (2.5) and (2.6) to evolve them from µ Λ up to µ m . Equivalently, we can think of evolving the products C x (Q, µ m )M xk (m b , µ m ) from µ m further down to µ Λ (with x = j, b). The final expression for the DIS cross section is then given by The full cross section on the left-hand side contains large logarithms ln(Q/m b ) and ln(m b /Λ QCD ). On the right-hand side these are factorized into logarithms ln(Q/µ H ) and ln(m b /µ m ), which are considered small and reside in the coefficients C i,b (Q, µ H ) and M(m b , µ m ), large logarithms ln(µ H /µ m ) and ln(µ m /µ Λ ), which are resummed into the evolution factors U [5] (µ H , µ m ) and U [4] (µ m , µ Λ ), and finally logarithms ln(µ Λ /Λ QCD ), which are absorbed into the PDFs at µ Λ . We will refer to eq. (2.22) as the resummed ("resum") result, since it has all logarithms ln(Q/m b ) resummed.
In the traditional 5F scheme, the resummed result in eq. (2.22) is written as where the combinations are interpreted as the evolved 5F PDFs including a PDF for the bottom quark f [5] b . To all orders in α s , eqs. (2.23) and (2.24) are simply a different way to write eq. (2.22). In practice, however, the evolution and matching corrections are always carried out to a certain finite order, where the different interpretations lead to different perturbative countings yielding different results. This is discussed in detail in section 2.3. Basically, in eq. (2.23) the 5F PDFs are traditionallly regarded as external O(1) inputs, and the perturbative order counting in α s is only applied to the coefficients C i and C b . In contrast, in eq. (2.22), we only regard the f [4] l (µ Λ ) as external O(1) quantities, while the perturbative order counting is applied to all terms in curly brackets. As we will see, one advantage of doing so is that this renders the order counting consistent between the resummed and fixed-order results, which facilitates their combination, as discussed in detail in section 2.4. Figure 4. Schematic leading-order matching at µ H for heavy-quark production in DIS for m b Q. Figure 5. Schematic NLO matching at µ H for heavy-quark production in DIS for m b Q.
We also note that the publicly available 5F PDF sets are constructed as in eq. (2.24), with the notable difference that the matching scale µ m is commonly identified with and fixed to the heavy-quark mass, µ m ≡ m b . However, it is clear from our discussion that µ m is a (in principle arbitrary) perturbative matching scale and it is important to keep it conceptually distinct from the parametric m b dependence. In our results, we will utilize the µ m dependence to estimate the intrinsic resummation uncertainties.

Matching at µ H ∼ Q
The Wilson coefficients C i (Q, µ H ) and C b (Q, µ H ) are computed in perturbation theory by taking partonic matrix elements of both sides of eq. (2.17), The calculation proceeds analogous to section 2.1.1. The essential difference is that now b quarks are present in the theory below µ H and so we also have to consider external b-quark states to determine the C b matching coefficient. As discussed earlier, the fulltheory matrix elements on the left-hand side are expanded to leading order in m 2 b /Q 2 . The partonic matrix elements on the right-hand side now lead to partonic PDFs similar to those of eq. (2.10), now including also b-quarks, To perform the matching, we expand both sides of eqs. (2.25) and (2.26) in powers of α s (µ H ) ≡ α [5] s (µ H ), where all the pieces are expanded analogously to eq. (2.11). The leading-order matching is illustrated in figure 4. At LO, the partonic bottom PDF is The limit m b → 0 explicitly highlights that the full-theory matrix element is expanded in The NLO matching is illustrated schematically in figure 5. At this order, there are 1-loop and real-emission corrections to the bottom-quark LO contribution as well as a contribution from a gluon channel. The partonic PDFs at NLO for finite m b are (see e.g. ref. [52]) Together with the LO b-quark coefficient in eq. (2.29), the NLO matching coefficients are The logarithms of m b inside the matrix elements of the DIS operator precisely match those in the partonic PDFs in eq. (2.30), such that the m b → 0 limit is finite. The reason is that for the matching at µ H , the finite bottom mass is nothing but an IR regulator for the collinear divergences associated with bottom quarks, which cancels in the matching.
Since the matching coefficients are independent of the IR regulator, we can also take the m b → 0 limit at the beginning, as long as we use another IR regulator, such as dimensional regularization. In this case, the computation of the coefficients C i,b becomes much simpler since there is one less scale involved. The partonic PDFs are then the usual ones in pure dimensional regularization, completely analogous to eq. (2.14), and P qg (z) as in eq. (2.31). The NLO coefficients are then given by 35) and are precisely the same as in eq. (2.32). Figure 6. NLO matching at µ m .

Matching at µ m ∼ m b
To compute the matching coefficients M ij at the low scale µ m , we calculate matrix elements of both sides of eqs. (2.20) and (2.21) with the same external partonic states. Using the definitions in eq. (2.10) and eq. (2.27), we have (2.36) Now, the b quark cannot appear anymore as an external state, since it is integrated out on the right-hand side. (Hence, there are no equivalent matching equations for f [5] b/b or f [5] i/b .) At the same time, m b is now the hard scale which appears in the matching coefficients (and cannot be set to zero). The matching coefficients M ij are known fully to O(α 2 s ) [53]. (They are also known partially to O(α 3 s ), see e.g. refs. [54,55] and references therein.) Expanding eq. (2.36), the LO matching is simply bg , which are illustrated in figure 6, Note that the precise number of flavors used in α s here is an O(α 2 s ) effect, which will then also generate nontrivial matching conditions for M (2) gg and M (2) bg .

Perturbative expansion and order counting
We now discuss the perturbative counting for the cross section for the two scale hierarchies in figure 1. Since the gluon and light quark PDFs at the scale µ Λ are nonperturbative objects fitted from data, we make the standard assumption and count them as external O(1) quantities, To determine the cross section to a certain perturbative accuracy, a perturbative counting should be applied to all remaining terms in the cross section that are computed in perturbation theory. In the fixed-order case m b ∼ Q, this implies that the standard perturbative counting in terms of powers of α s appearing in the hard matching coefficients applies. On the other hand, for m b Q, the perturbative counting should be applied also to the matching coefficients at µ m and the evolution factors between µ H and µ m . We will argue that for phenomenologically relevant hard scales this implies that an appropriate perturbative counting takes the effective bottom PDF to be an O(α s ) object.

Fixed order
First, recall that the DGLAP evolution factors U ij (µ 1 , µ 2 ) resum single logarithms of the ratio µ 1 /µ 2 to all orders in α s . For this purpose, one performs a logarithmic counting where one expands in powers of α s while counting α s ln(µ 1 /µ 2 ) ≡ α s L ∼ 1. That is, one formally counts L ∼ 1/α s . We can then write 40) where U N k LL are functions of α s L to all orders in α s . Combining this with eq. (2.39), we can also count the evolved 4F PDFs as O(1) quantities The PDFs evolved at N k LL are then usually called N k LO PDFs. Note that while the evolution mixes the PDFs, it does not induce a parametric difference between the lightparton PDFs. Also, in the limit µ 1 → µ 2 we have U [4] ij → δ ij . Hence, we can generically treat f [4] i (µ) as external O(1) quantities for any µ regardless of how large the logarithms ln(µ 1 /µ 2 ) actually are, and this is the standard praxis.
For the fixed-order (4F) cross section in eq. (2.8), the perturbative counting in α s is then directly applied to the D i coefficients, so it has the perturbative expansion [with a s ≡ α [4] s (µ H )/(4π) as in eq. (2.11)] (2.42) Taking into account eq. (2.40), obtaining the cross section at an accuracy of order α k s (N k LO) then requires the N k LO matching coefficient with the N k LL evolution (i.e. N k LO PDFs). Here the order is counted relative to the lowest nonvanishing order, which for heavy-quark production in DIS is O(α s ). When combining the expansion in eq. (2.40) with the α s expansion of D i (Q, m b , µ H ), one could in principle reexpand the product of the two series. In practice, this is usually not done, since the f [4] i (µ H ) are treated as external O(1) inputs as mentioned above.

Resummation
We now discuss the perturbative counting for the resummed cross section in eq. (2.22). First, we can use the same arguments as in eq. (2.41) to treat the 4F PDFs at µ m as O(1) inputs. We then have to consider the perturbative counting for the terms in curly brackets in eq. (2.22). The situation is more subtle now due to the presence of the additional scale µ m ∼ m b . Depending on the hierarchy between µ H and µ m , there are two different options of how to count the evolution factors U [5] ij (µ H , µ m ).
• For very large hierarchies, we can use a strict logarithmic counting, in which case α s L ∼ 1 and eq. (2.40) generically applies to all the evolution kernels so U [5] ij ∼ U [5] bj ∼ U [5] ib ∼ U [5] bb ∼ 1 . (2.43) • For intermediate hierarchies µ m µ H , the resummation can still be important, so we still use eq. (2.40) to organize the logarithmic order of the resummation. However, we should also take into account that in the limit µ m → µ H the off-diagonal mixing evolution kernels vanish U [5] bg (µ m → µ H , µ H ) → 0 and similarly for U [5] gb . This is because their fixed-order expansion starts at order α s L rather than 1, so they are suppressed by an overall factor of α s L relative to the diagonal U [5] bb and U [5] gg . Therefore we count The counting in eq. (2.43) corresponds to the traditional 5F scheme. With this counting and using eqs. (2.37) and (2.38), the evolved 5F PDFs in eq. (2.24) have the perturbative expansion Hence, they are treated as external O(1) quantities. The resummed result is then written as in eq. (2.23) and has the perturbative expansion [with a H = α [5] (2.46) The N k LO cross section then requires using the N k LO 5F PDFs, which are given by the expansion of eq. (2.45) to O(α k s ) together with the N k LL evolution factors. A rough numerical estimate shows that for µ m ∼ m b ∼ 5 GeV the first case eq. (2.43) applies for hard scales µ H 1 TeV. Thus, the second case in eq. (2.44) is more appropriate for our purposes. This is also confirmed by the fact that for µ H ∼ O(100 GeV) and standard PDF sets one finds that numerically f [5] b (µ H ) f [5] g (µ H ). Adopting this counting, the resummed result in eq. (2.22) has the perturbative expansion bg + a m U [5] bb M (1) bg bg + a m U [5] bb M (1) bg Here, a H ≡ α [5] s (µ H )/(4π) and a m ≡ α [5] s (µ m )/(4π), and for notational simplicity we have suppressed the convolution symbols and all arguments (which are as in eq. (2.22)). In the contributions proportional to f [4] q we have also counted U gq ∼ α s and U bq ∼ α 2 s . Note that, in the region where this counting applies, a H and a m can be regarded as being parametrically (and practically) of the same size.
From eq. (2.47) we see that the counting in eq. (2.44) leads us to include the matching terms C (1) g and M (1) bg , which provide the boundary conditions for the RGE, already at the lowest order, i.e. one order lower compared to the 5F. Furthermore, any cross terms in eq. (2.47) from the matching at µ H and µ m are expanded against each other. In other words, compared to eq. (2.46), we do not have overall 5F PDFs, but rather the contributions ∼ U [5] ij (µ H , µ m ) ⊗ M jk (µ m ) making up the 5F PDFs in eq. (2.45) are expanded together with the hard matching coefficients. As we will see in the next subsection, these features enable us to have an easy and smooth transition to the fixed-order result. Note that this is quite similar to how the primed resummation orders N k LL are implemented in the resummation for differential spectra, see e.g. refs. [56][57][58][59], where this facilitates a clean and smooth transition to the fixed-order result.
We can of course collect the terms proportional to C b and C i in eq. (2.47) into effective PDFs, which we denoted asf b andf i to distinguish them from the standard 5F PDFs in eq. (2.45). With the counting in eq. (2.44) we then havẽ Thus, in the region of scales we consider, the effective b-quark PDF should be treated as an O(α s ) object, while the gluon PDF still starts at O(1). Note though that the O(α s ) terms inf g are different from those in f [5] g . For example, the U [5] gb ⊗ M (1) bg term in eq. (2.45) counts as O(α s ) in f [5] g while it only appears at O(α 2 s ) inf g . Since the light-to-light M gg and the light-to-heavy M bg matching functions are needed at different relative orders iñ f b (m b , µ H ), this definition of the effective bottom PDF differs with respect to the usual f [5] b (µ H ) in eq. (2.45). Hence, in our numerical implementation we cannot use the b-quark PDF from the standard 5F PDF sets. Instead, we need to constructf b (m b , µ H ) ourselves. We do so by creating PDF grids that have the matching coefficients at the required order but the same order in the evolution factors. The technical details are discussed in Appendix B.
Denoting withf Here, we still consistently drop any higher-order O(α 3 s ) cross terms in the product of coefficients and effective PDFs by keeping the effective PDFs to different orders in the different terms. As already mentioned, this is important to ensure a smooth transition to the fixed-order result in the limit µ m ∼ µ H .
On the other hand, for the bbH hadron collider process we are eventually interested in in section 4, we will have two PDFs and the practical implementation of the strict expansion gets quite involved. Therefore, as long as we are only interested in the phenomenologically relevant region µ m ∼ m b µ H ∼ m H , we can also keep the higher-order cross terms to simplify the practical implementation. We then have (2.50) Once we allow keeping higher-order terms, we can also further simplify the practical implementation by replacing the effective PDFsf i,b above by standard 5F PDFs f [5] i,b . These must then be of sufficiently high order such that they include all necessary matching corrections as required by our perturbative counting. However, we note that whenever one keeps higher-order terms for practical convenience, one should check that this does not have a large numerical influence on the results in the kinematic region of interest. We will come back to this in section 4.
We stress, that even when keeping higher-order cross terms, the perturbative counting is still performed for both C i,b andf i,b withf b counted as O(α s ). So even though the leading term in C s ), the resummed result starts at O(α s ). Comparing to eq. (2.42), the resummed result has a perturbative counting consistent with the fixed-order result. It precisely corresponds to a resummed version of the fixed-order result in the m b → 0 limit. This organization and implementation of the resummation is one of the main ways in which our approach differs with other approaches. This will be discussed further in section 3.

Combination of resummation and fixed order
In sections 2.1 and 2.2 we have derived results for the heavy-quark production cross section in DIS that are relevant for two different parametric scale hierarchies. The fixed-order result dσ FO in eq. (2.8) is relevant for m b ∼ Q, as it keeps the exact m b dependence to a given fixed order in α s , but does not include the all-order resummation of logarithms ln(m b /Q). The resummed result in eq. (2.22) is relevant for m b Q, as it resums the logarithms ln(m b /Q) to all orders in α s , but neglects any m b /Q power corrections that vanish for m b → 0. These two results represent two ways of computing the same cross section. In this section, we combine these two results and obtain our final result accurate for any value of We follow the usual approach for combining a higher-order resummation with its corresponding fixed-order result. We write the full result for the cross section as dσ = dσ resum + dσ nons . (2.51) Here, the nonsingular cross section dσ nons contains all contributions that are suppressed by O(m b /Q) relative to dσ resum and vanishes in the limit m b → 0. With this condition, dσ automatically contains the correct resummation in the m b → 0 limit. Furthermore, we require that the fixed-order expansion of eq. (2.51) reproduces the correct FO result, including the full m b dependence. Therefore, where the singular contributions dσ sing are obtained from the fixed-order expansion of the resummed result to the desired order in α s . For dσ nons to indeed be nonsingular and vanish for m b → 0, dσ sing must contain all singular contributions in dσ FO , i.e. all terms that do not vanish as m b → 0. This in turn requires that the resummation to a given order fully incorporates all these fixed-order singular terms. In this sense, the resummed result should be consistent with the fixed-order result. This condition is precisely satisfied by our resummed result with the perturbative counting used in eqs. (2.47) and (2.50), for which the (N)LL result contains the full (N)LO singular terms, as we will see below.
To explicitly identify the nonsingular terms, we need a meaningful and consistent comparison between dσ FO and dσ sing , which means we have to write both in terms of the same external 4F PDFs and expand both in terms of the same α s . For this purpose, it is most convenient to use f [4] i (µ H ) as in eq. (2.42) but perform the expansion in terms of α [5] s (µ H ) as in the resummed result eqs. (2.47) and (2.50). First, for dσ FO , we can simply change the b-quark renormalization scheme for α s used to computed the D i matching coefficients in eq. (2.42) from the decoupling scheme to the MS scheme. This leads to modified Wilson coefficients which are now expanded in terms of the same α [5] s (µ H ) as is used in C i (Q, µ H ), C b (Q, µ H ). From the point of view of the fixed-order calculation, this is actually the more appropriate expansion for m b < µ H . Next, the singular cross section can be easily obtained by evaluating the resummed result in eq. (2.47) or eq. (2.49) at µ m = µ H , Finally, the fixed-order nonsingular cross section is given by At each order in α s , all singular terms in D MS i are exactly cancelled by the corresponding singular terms from the resummed result, such that dσ nons is free of collinear logarithms and vanishes as m b → 0.
We stress that the statement dσ resum | FO = dσ resum | µm=µ H utilized above is quite nontrivial and crucially relies on the fact that with our perturbative counting in the resummed result all the matching corrections M ij are always included to sufficiently high order (basically to the same order in α s to which we have to expand the evolution kernels) such that the µ m dependence precisely cancels in dσ sing to the given order in α s to which we expand. Once we know that this is the case, we can pick any µ m we like to perform the fixed-order expansion of dσ resum . The choice µ m = µ H is then the most convenient, since all the evolution kernels become trivial. For example, at LL we have and therefore  − · · · ] ⊗ f i . Therefore, we can naturally associate them with the gluon and light-quark PDFsf i . We can then absorb the nonsingular corrections into the light-parton coefficient functions by taking while keeping the C b coefficient unchanged. Equivalently, we can replace C i →C i everywhere and impose the condition such that eq. (2.54) vanishes.
The above shows that we can choose to absorb the nonsingular contributions into the resummed result by modifiying the matching coefficients at µ H . The condition in eq. (2.60) implies that the light-parton coefficientsC i (Q, m b , µ H ) can be obtained from the matching at µ H in section 2.2.1 without taking the m b → 0 limit in the light-parton full-theory matrix elements, while for all bottom contributions and coefficients the m b → 0 limit is still taken.
We can now write the final result for the cross section as where the choice of expanding the cross terms or not is kept implicit and will determine to what order the PDFs are kept in each term. We emphasise that in this form the result is very convenient to implement, since it essentially only requires the fixed-order result (after changing the b-quark renormalization scheme for α s ) and the massless resummed result. In section 4 we will apply this strategy to the bbH cross section and provide further details on the construction of the coefficient functions. By choosing to absorb dσ nons into the matching coefficients in the resummed result, we effectively let the leading-power resummation also act on the nonsingular corrections. This introduces power-suppressed higher-order logarithmic terms, which however are beyond the order we are working at. In particular, this does not include the correct resummation of power-suppressed logarithmic terms. (This would require the extension of dσ resum to subleading order in m b /Q, which is well beyond the scope of this work, and also very likely irrelevant at the current precision.) Fundamentally, we only have control over the nonsingular corrections at the level of their fixed-order expansion. The above procedure to include the nonsingular contributions is not unique, and while physically motivated, is ultimately driven by practical convenience. We would like to underline that alternative choices are in principle possible, provided they do not change the resummation in the m b → 0 limit and reproduce the correct fixed-order expansion, in which case they will effectively differ by power-suppressed higher-order logarithmic terms. 7 We will come back to this point in section 3.
From the discussion so far, it is clear that transition between dσ resum and dσ FO is controlled by the scale µ m . To provide a smooth transition between the resummation and fixed-order regions, this scale is promoted to a m b -dependent profile scale µ m → µ m (m b , µ H ). It has the properties that in the resummation region for m b Q it has the canonical resummation scaling µ m ∼ m b , while in the fixed-order region m b ∼ Q it approaches µ m → µ H , such that the resummation is turned off there and the fixed-order result is recovered, with a smooth transition in between. The fact that it is possible to control this transition between limits with a single scale, makes our predictions in the transition region robust and, moreover, variation of this scale and of its functional form provides a solid handle on the associated theoretical uncertainties. The precise definition and variations of the profile function are discussed in detail in section 4.3 for the case of bbH production.

Comparison to existing approaches
In the previous section we used a systematic field-theory analysis to derive a result for the heavy-quark production cross section in DIS accurate for all possible scale hierarchies from m Q to m ∼ Q. Various approaches to the same problem are available in the literature, which go under the name of variable flavor number schemes (VFNSs). In this section, we briefly compare the existing schemes to our result. In what follows, we do not attempt to give an in-depth review of the different schemes but rather focus on similarities and differences with respect to our EFT result. For reviews of the different schemes in the literature we refer to refs. [60,61] and sect. 22 of ref. [62].
VFNSs can in principle differ in various aspects. The first is the general construction, namely how resummation of collinear logarithms is achieved and how nonsingular power corrections are included. Secondly, they can differ in how the perturbative counting is performed, that is, which of the various perturbative ingredients are included at a given order. Third, they can differ in how the heavy-quark threshold is implemented, which in our language corresponds to the exact choice of the low matching scale µ m . The first aspect is the one that primarily distinguishes the different schemes, while the remaining two aspects are more related to choices made within each scheme. Here, we compare to the choices often used in the literature. We stress though that these choices correspond to how a particular scheme has been used or implemented in practice, but (in most cases) they do not necessarily represent restrictions of a particular scheme itself.

Construction
We start by discussing the differences in the basic construction of the cross sections. For m b Q ("below threshold") all schemes use the same fixed-order 4F result in eq. (2.8). For m b Q ("above threshold") the various schemes construct their cross sections as follows: • Zero-Mass (ZM). In this approach, the massless resummed result in eq. (2.23) is used for m b Q, while nonsingular power corrections are neglected at any order in α s . Hence, this scheme is only expected to be accurate for m b Q. Since powersuppressed contributions are not included, it is not accurate close to the heavy-quark threshold and does not reproduce the full fixed-order result. For this reason, we do not discuss it further.
• ACOT [1,2,5]. The ACOT scheme is based on the idea that the power corrections can be fully included in DIS at the level of the matching at the hard scale µ H , eq. (2.17), by generalizing it such that power corrections in m b /Q are included in the definition of the Wilson coefficients. The heavy quark is considered as an active flavor and the quark mass dependence is retained at each matching step, yielding incorporate the nonsingular contributions and reduce to the original C i,b (Q, µ H ) in the m b → 0 limit. In contrast to eq. (2.61), the heavy-quark contributions in eq. (3.1) are computed with a massive on-shell heavy quark in the initial state. To account for the massive kinematics including the presence of a massive (unresolved) heavy quark in the final state, the heavy-quark Bjorken-x can be rescaled, leading to a variant of this scheme called ACOT-χ [7,63]. While the validity of ACOT in DIS can be based on including heavy-quark masses in the hard-scattering factorization [5], its extension to the case of two incoming hadrons is problematic due to the massive kinematics, see Appendix C.
• S-ACOT [6]. The fact that f [5] b is not independent of f [5] i (see eq. (2.24)) allows one to move power corrections betweenC b andC i without spoiling the formal accuracy of eq. (3.1) [5,6]. This was used to construct a simplified variant of ACOT, in which the heavy-quark Wilson coefficients are computed in the massless limit,C b (Q, m b , µ H ) → C b (Q, µ H ), while the full mass dependence is retained in the (modified) light-parton coefficients. This is evidently equivalent to how we include the nonsingular corrections in eqs. (2.60) and (2.61) for practical purposes. To account for the massive heavyquark kinematics, a χ-rescaling is also applied, leading to S-ACOT-χ [7,13]. 8 In ref. [31], a modification of ACOT, dubbed m-ACOT, is used for the case of two incoming hadrons, where the massless limit is applied only to channels with two incoming heavy quarks, while the mass dependence is kept in heavy-light and lightlight channels.
• TR [3,8]. The TR scheme is defined by requiring that the fixed-order result, after being expressed in terms of 5F PDFs, corresponds to the resummed result up to power-suppressed contributions. This requirement fixes the singular contributions. However, there is still freedom for the treatment of nonsingular terms, and this is fixed by making a choice such that the coefficient functions obey a sensible threshold limit. The result is hence different from both ACOT and S-ACOT. Due to the choice of perturbative counting, a discontinuity exists at threshold which is removed by adding a Q-independent contribution to the result above threshold. Though this contribution is formally higher order, it can be sizeable, even far from threshold. The presence of the constant terms complicates the generalization of this scheme to higher-orders and to hadron-hadron collisions.
• FONLL [11,32]. This scheme is constructed by adding the massless resummed result to the full fixed-order result and consistently subtracting the double counting orderby-order in α s . The fixed-order contribution is rewritten in terms of 5F PDFs with the resulting ambiguity fixed through the choice that only light channels contribute, as we have also done in section 2.4. The double-counting terms are equivalent to the singular terms in our notation, and remove from the fixed-order result its massless limit, i.e. all its terms that do not vanish in the m b → 0 limit. The FONLL procedure is thus equivalent to adding the dσ nons to the resummed result. This also makes the FONLL construction formally equivalent to S-ACOT. Finally, a damping factor, which performs the same function as the χ rescaling in S-ACOT, is used to suppress higher-order spurious contributions and guarantee continuity at threshold.
From the point of view of the all-order resummation, all these schemes are equivalent, as they all include the same resummation. As discussed in section 2.4, the minimal and formally correct result above threshold is given by eq. (2.51) as dσ = dσ resum + dσ nons . This result is formally correct in the sense that it correctly resums the collinear massive logarithms and correctly includes the full mass dependence and kinematics at fixed order. It is minimal in the sense that the nonsingular corrections dσ nons are unambiguous and unique when written in terms of 4F PDFs as in eq. (2.54), and are strictly included at fixed order, while the resummation strictly only includes leading-power terms.
As discussed in section 2.4, there is an ambiguity when one tries to partially or fully absorb the nonsingular contribution into the resummed result, which amounts to expressing them in terms of 5F PDFs. The primary perturbative ingredients of ACOT, S-ACOT, TR, and FONLL are the same and they only differ in the way by which they fix this ambiguity. This ambiguity corresponds to power-suppressed higher-order logarithmic terms. Hence, these schemes can be regarded as formally equivalent up to such terms, which are beyond the considered formal accuracy.

Combination of the ingredients
We now move to the second source of scheme differences, namely how the perturbative counting is performed. By construction, the coefficient functions of (S-)ACOT, TR, and FONLL differ from each other and to those in our EFT result by formally higher-order contributions. Therefore, the largest differences between the approaches arise from the perturbative counting. In all practical implementations we are aware of, the perturbative counting used by each scheme is as follows: • ACOT-like schemes, used in the CTEQ family of PDF fits, construct perturbative expansions in the usual way by counting explicit powers of α s in the coefficient functions. As a result, for DIS at LO (α 0 s ) the result below threshold is zero, while above threshold it is nonzero due to the heavy-quark initiated contributions (C • The TR scheme, used in MSTW and HERAPDF fits, is somewhat different as it combines the orders such that the lowest nonvanishing order below and above threshold appear at the same time. This means that at LO the result below threshold is the O(α s ) gluon-initiated contribution, while above threshold it is the O(α 0 s ) heavy-quark initiated contribution. The additional Q-independent term added above threshold is formally of higher order and does not affect the counting.
• FONLL, used in NNPDF fits, also adopts the standard perturbative counting. The NLO and NNLO results are called FONLL-A and FONLL-C respectively. There is an intermediate result, FONLL-B, where the fixed-order terms are computed to order α 2 s (NNLO) but the massless contribution is only included at order α s (NLO).
None of the schemes discussed above adopts a perturbative counting which is directly comparable to our approach of performing the counting on the full perturbative part of the cross section including both evolution and matching. In particular, it implies that the Only representative diagrams at a given order and for a given channel are shown. Notice that, while the diagrams appearing in the Q < µ m boxes contain collinear logarithms due to the heavy quark, the latter are subtracted in the diagrams appearing in the Q > µ m boxes. c 0,1 are the Q-independent terms present in the TR-scheme that ensure continuity at threshold. We also point the reader to a very similar table in ref. [61].
effective heavy-quark PDF should be counted as an O(α s ) object. Since the perturbative counting used by the different VFNSs summarized here does not distinguish between heavy-quark and light-parton initiated contributions, this difference in order counting is a principal difference in our approach. In figure 7 we summarize the perturbative counting adopted in the (S-)ACOT, TR and FONLL schemes as well as the counting we propose. As argued in section 2.3, our order counting is well justified theoretically and appropriate for a wide range of scales (including scales appropriate for DIS experiments in the case of both bottom and charm quarks). It also has several advantages. As we will see in section 4, one of these is that the perturbative convergence tends to be improved with reduced uncertainties from the hard-scale matching. Another advantage, as highlighted in section 2.4, is that it facilitates a smooth transition to the fixed-order result at the heavyquark threshold (provided the counting is strictly applied and higher-order cross terms are neglected), without the need of any rescaling or damping factors.

Matching scale dependence
Finally, the last important difference between the existing schemes and our approach is the position and treatment of the heavy-quark threshold. In all applications we are aware of, the threshold is always set equal to the heavy quark mass, i.e. effectively the resummation scale µ m is fixed to µ m = m b . This is both the scale at which the low-scale matching is performed and also the scale at which one switches from the fixed-order result to the resummed result.
Recently [61,64], it has been suggested to consider an additional switching scale µ S > µ m , at which the computation switches from a fixed-order result to a resummed one, but nevertheless keeping the matching at a different (lower) scale. The effect of this choice is to delay the use of the resummed result, perhaps to a region where mass effects are negligible, though the transition between the resummation and fixed-order regions is not guaranteed to be smooth.
In this work, we exploit the dependence on the matching scale µ m to explicitly control the transition to the fixed-order result as m b → Q, and furthermore to estimate the intrinsic perturbative uncertainty in the resummation and matching procedure. This is in fact the standard practice in resummed calculations involving different resummation scales. This uncertainty should be taken into account as part of the total perturbative uncertainty in the result, which is typically not the case in existing approaches.

Higgs production in association with b quarks
In this section, we extend the framework presented in section 2 to hadron-hadron collisions and apply it to the bbH process, i.e. Higgs-boson production in association with b quarks. Specifically, this process can be defined as Higgs production via the bottom Yukawa coupling Y b , with all other Yukawa couplings set to zero. (As discussed in section 4.2, we do not include the b-quark loop contributions that are usually included in the gluon-fusion process. There we also comment on the inclusion of Y t Y b interference terms that are usually regarded as part of the bbH process.) The bbH process makes up only a tiny fraction, ∼ 1%, of the total Higgs production cross section in the Standard Model (SM). It is nevertheless an interesting process within the SM, since the total bbH cross section is comparable to the total pp → ttH cross section for LHC energies, and because it provides direct access to the bottom Yukawa coupling. Furthermore, this process may be sensitive to new physics effects, since in many BSM scenarios, such as two-Higgs-doublet models with large tan β, the Higgs coupling to bottom quarks can be enhanced.
In the SM, the cross section in the massless 5FS is known at NLO [26,27] and NNLO [28,29], and in the 4FS at NLO [24,25]. NLO predictions matched to parton showers for bbH production have been studied in both 4F and 5F schemes [65]. The 4FS and 5FS calculations can lead to very different results, with cross sections differing by as much as an order of magnitude. For appropriate choices of the factorization scale, the difference can be reduced significantly, leading to more compatible results within the perturbative uncertainties, see the discussions in Refs. [21, 23-25, 28, 66]. As discussed already, the 5FS and 4FS possess different merits, and predictions that combine the advantages of both are highly desirable. The current combined values by the LHC Higgs Cross Section Working Group [67,68] are obtained using the Santander matching prescription [30], which amounts to a weighted average of the cross sections obtained in the two schemes. In contrast, our predictions here are derived from a consistent field-theory setup, and can thus be regarded as a definite improvement over the currently used prescription.
We start in section 4.1 by extending the EFT result of section 2 to the case of two incoming protons. In section 4.2, we give details about the practical setup of our results. In section 4.3, we discuss our procedure to obtain robust estimates of the perturbative uncertainties from separate variations of the µ H and µ m matching scales. In particular, we discuss the profile scales and variations for the matching scale µ m . In section 4.4, we present our result for the bbH cross section as a function of the b-quark mass. This serves as a validation of our matching procedure, confirming that our approach satisfies all the required properties. There, we also discuss the size of nonsingular power corrections suppressed by m b /Q. Finally, in section 4.5, we present our final results at the physical b-quark mass for several Higgs masses and compare to the existing results obtained in the 4FS, 5FS, and the Santander prescription.

Extension of the EFT approach to hadron-hadron colliders
The simplicity of the EFT framework presented in section 2 for DIS makes it possible to straightforwardly extend the setup to the case of two incoming protons. This is certainly not the case for any of the schemes discussed in section 3, whose consistent generalization to hadron-hadron collisions can be highly nontrivial. We first point out that the evolution of the quark operators and the matching at µ m are identical and therefore the evolved PDFs of eq. (2.24) are the same. Of course, the matching at the hard scale is different.
i (µ H )O [5] j (µ H ) + C bk (m H , µ H ) O [5] b (µ H )O [5] k (µ H ) + O [5] b (µ H )O [5] k (µ H ) + O [5] k (µ H )O [5] b (µ H ) + O [5] k (µ H )O [5] b (µ H ) where we have introduced a bbH operator O bbH , and m H is the Higgs mass. Note that we have used the identities C bk = Cb k = C kb = C kb , C bb = Cb b and C bb = Cbb. These two results are the straightforward extensions of the results in eqs. (2.2) and (2.17), where, as discussed in section 2.4, we have made the choice to absorb power corrections into the coefficients for the light channels.
As in the DIS case, in our order counting we take the bottom PDF to be an object of order α s . As discussed in section 2.3, this is appropriate for a hard scale of the order of the Higgs mass, µ H ∼ m H (more generically for µ H 1 TeV). Therefore, up to NLO, the fixed-order result our cross section matches into for µ m → µ H is given by (omitting the arguments for simplicity) For µ m < µ H the resummed and matched cross section is written as 10 where as in eq. (2.62) we have left implicit the strict expansion of the products of effective PDFs and coefficient functions. We notice that in both cases, using the perturbative order counting introduced in section 2.3, LO(+LL) is in fact order α 2 s and NLO(+NLL) includes the order α 3 s corrections. In section 4.4, we discuss the implementation and results of a strict expansion of eq. (4.4) as well as a more practical implementation keeping higher-order cross terms and using standard 5F PDFs.
The 4FS result corresponds to the result of eq. (4.3) used for all scale hierarchies and where the decoupling scheme is used for the b-quark renormalization of α s . The 5FS result on the other hand corresponds to the massless limit of eq. (4.4), replacingf i,b with f [5] i,b and with the perturbative order counting performed only on the coefficient functions (i.e. assuming the bottom PDF of order 1). This has the expansion (1) bb f [5] b f [5] b + a H 4C (1) bg f [5] b f [5] g NNLO (5F) + a 2 H 2C (2) bb + 2C (2) bb f [5] b f [5] b + a 2 H 4C (2) bk f [5] b f [5] k + a 2 H C (2) ij f [5] i f where C (2) ij are the massless coefficients (namely the massless limit ofC (2) ij ). In figure 8 we illustrate the different countings diagrammatically. This highlights that one can regard our results as a resummation-improved 4FS result.
The massless coefficients C bg and C (2) bk required to reaching NLO+NLL accuracy in our result, eq. (4.4), are the same as those of the massless 5FS computation and can be found explicitly in ref. [28]. Trivially extending eq. (2.60) to the case of two light-light light-heavy heavy-heavy Figure 8. Sample diagrams appearing in the computation of the Higgs production cross section in association with b quarks. Diagrams are grouped according to the different countings adopted in our resummed result, eq. (4.4), and in the 5FS result, eq. (4.5). The 4FS counting coincides with the resummed counting in the fixed-order limit where only the diagrams in the first column are considered.
initial-state legs, the matching coefficientsC (2) ij andC (3) ij can be written as, The MS massive coefficients can be obtained from the decoupling-scheme coefficients as described in Appendix A: ij . (4.7) We have implemented analytic expressions for the coefficients D (2) ij in an in-house code and extract the numerical result for D (3) ij from Madgraph5 aMC@NLO [69], after generating the process pp → bbH at NLO. We have explicitly checked that our implementations, including pole scheme to MS scheme changes for Y b in D ij , exactly reproduce the inclusive results of the bbh@nnlo code [28] and of the recent bbH studies of ref. [65].

Setup
Here we summarize the set of input parameters we use to produce the results of sections 4.3 and 4.4. Unless indicated otherwise we always use the setup detailed below.
Collider energy We provide predictions for the LHC at √ s = 8 TeV.
PDFs We have created PDF sets using a modified version of APFEL [70] for the evolution from a fixed low scale where the parametrization of a known PDF set (MSTW2008) has been used. The main reasons for our modifications were the implementation of a general value for the threshold matching scale µ m as well as the generation of the effective PDFsf {k} required in a strict expansion of eq. (4.4). Further details are given in Appendix B.
Higgs mass We use m H = 125 GeV as default.
Bottom mass For all results where the bottom mass is fixed to its physical value, we use a pole mass of m b = 4.75 GeV for the kinematic mass scale that enters in the 4F matrix elements and in the low-scale matching coefficients M ij . For the Yukawa coupling we use the MS mass m b (m b ) = 4.16 GeV as input, see also below. The use of different bottom masses for the Yukawa coupling and in the matrix elements is not unsualthis has been the setup of the 4FS calculations of Ref. [24,25]. What is different to previously used setups (and to the LHCHXSWG) is that the two values we use are not related to each other via a one-loop conversion. This is not a problem, since the two perturbative series in which they enter are unrelated. 11 What is relevant in our case is that we consistently use common values in both the resummation and fixed-order parts of the calculation. The numerical values above are chosen to have reasonable physical values and to enable an as consistent as possible comparison with the default 4FS and 5FS results.
In our results where we vary m b to study the dependence on the bottom mass, m b and m b (m b ) are varied consistently, with the conversion between the two at one loop as required for our NLO calculation.
Yukawa couplings All the Yukawa couplings are set to zero except the bottom quark Yukawa, Y b . The bottom Yukawa is renormalized in the MS scheme and its running is set to 4 loops. In our numerical studies we always evaluate it at the hard scale µ H , which is the appropriate scale for the resummation of large logarithms ln(µ H /m b ) associated with the bbH vertex and hence leads to better perturbative convergence [71].
Bottom loops Contributions to Higgs production where the Higgs couples to a closed b-quark loop are usually included in the gluon-fusion cross section, since their most important effect is due to the interference of the bottom loop with the top loop. As usual, we exclude these contributions from our bbH computation, such that the result has no double counting with the gluon-fusion cross section. Our result still includes bottom loop contributions, but only in diagrams with two bottoms in the final state, (not included in the gluon-fusion cross section) as part of the NLO correction to the gg channel in our result.
Y b · Y t interference In our results we neglect the interference contribution proportional to Y t · Y b by setting Y t = 0. In the SM, this correction is known to be important and reduces the inclusive 4FS NLO cross section by roughly 10% at the LHC for m H = 125 GeV [24,25,65], while in BSM scenarios with large tan β its relative contribution can be much smaller. This interference has been computed in the 4FS where it first enters at NLO via diagrams containing a top-quark loop, whilst in the 5FS up to NNLO this interference does not contribute [28]. For comparisons between 4FS and 5FS predictions it is often preferred that the interference terms are dropped [30,72] since the latter are not present in the 5FS. To better compare with the results in the literature we also make this choice here. However, we emphasize that the Y t · Y b terms can be straightforwardly and consistently included as an additional nonsingular fixed-order piece in our result. To do so, we can simply allow for a nonzero top Yukawa in the fixed-order coefficients D ij . No changes to the resummed part of our result are required at the order we are working.

Scale dependence and theory uncertainties
In this subsection, we discuss in detail the perturbative uncertainties in our results. We begin by looking at the hard scale dependence. We fix m b to its physical value, set µ m = m b , and plot in figure 9 the cross section obtained according to the 4FS (at LO and NLO, obtained using the code of [65]), the 5FS (at LO, NLO, and NNLO, obtained using bbh@nnlo [28]) and our result (at LO+LL and NLO+NLL). As expected, a clear reduction of the scale dependence is observed in all results when moving to higher orders. We also notice that the patterns of scale dependence of the 4FS (green dashed) and 5FS (blue dotted) results are opposite to each other with the former decreasing and the latter increasing with increasing µ H (except at NNLO). This is due to the fact that at LO the scale dependence is dominated by α s for the 4FS result (which clearly increases at small scales), while for the 5FS result it is driven only by the bottom PDF, which vanishes at the bottom threshold and therefore drops rapidly as the scale decreases. Therefore, over a wide range of hard scales the two results differ significantly.
In contrast, the framework we have presented in section 2 leads to cross sections (red solid) that are less sensitive to the choice of the hard scale, even at LO. The reason behind this is a large compensation between the contributions from the bb, bk, and ij channels. This is due to the fact that each b-initiated contribution compensates the collinear subtraction in a gluon (or light quark) initiated contribution and close to the heavy-quark threshold these terms are all of the same order. This leads to a scale dependence in the (N)LO+(N)LL results that has a similar pattern to that of the unresummed 4FS result, however the resummation of collinear logarithms significantly stabilizes the dependence on µ H . As with the 4FS, the 5FS results also have a greater dependence on µ H compared to our resummed results. The reason for this is that the 5FS predictions adopt a standard perturbative counting and thus the compensation observed in the EFT results is not present.
Additionally, figure 9 illustrates that a smaller scale µ H ∼ m H /4 leads to a more stable perturbative expansion for all the results, and also leads to better agreement between the different approaches. The reason for this has been studied in ref. [23] by a careful investigation of the actual size of the logarithms that arise in the 4FS prediction.
Next, we discuss the choice of µ m and its associated perturbative uncertainties. For this purpose, it is important to identify the kinematic region where the resummation is important and where it must be turned off. To this end, in the left plot of figure 10 we show the fixed NLO result and its decomposition into singular eq. (2.53) and nonsingular eq. structure.
In the m b → 0 region, the singular terms clearly dominate, while the nonsingular corrections are suppressed by at least an order of magnitude and tend to zero. This is the resummation region, where the canonical choice µ m = m b is appropriate to resum the large logarithms ln(m b /m H ) in the singular corrections.
With increasing m b the singular contribution starts deviating from the full result, crosses it at around m b ∼ 30 GeV ∼ m H /4, and becomes much larger than the full result in the large-m b region. This large-m b region corresponds to the fixed-order region, which exhibits a delicate balance between singular and nonsingular contributions, with a large cancellation between the two yielding the full result. This means that the distinction into singular and nonsingular is meaningless here. To not spoil this cancellation it is imperative that the resummation is switched off completely, which is done by taking µ m = µ H . The fixed-order region starts at m b m H /4, where the magnitude of both singular and nonsingular is larger than the full result, so there is clearly an O(1) cancellation between them. We have verified that this pattern holds at both LO and NLO and upon variation of the hard scale in the range m H /16 < µ H < m H . We can therefore safely take m b ∼ m H /4 as the point where we should turn off the resummation for any configuration we might consider.
A smooth transition between the canonical value µ m = m b in the resummation region and µ m = µ H in the fixed-order region is achieved by using profile scales [56,57], where the scale µ m is promoted to a function of m b , which smoothly interpolates between these two limits. The use of profile scales is a common practice when performing resummation in EFTs based on RGEs. Following refs. [58,59,73], we choose different sets of profiles that allow us to separately estimate fixed-order and resummation uncertainties, which in the end are added in quadrature.

For our central scales we use
with the profile function where α ∈ [−1, 1] and (4.11) The multiplicative factor f vary (x) tends to 2 in the limit x → 0 and tends to 1 in the limit x → x 3 (as before, we use x 3 = 1). The effect of this factor (when varying α ∈ [−1, 1]) is to vary the arguments of the resummed logarithms in the small m b region by a factor of two, while keeping the hard scale fixed. Hence, we can use these variations to estimate the resummation uncertainty ∆ resum . In the limit x → x 3 (or m b → µ H ) the effect of this variation tends to zero, as it must, and thus the resummation uncertainty vanishes in the fixed-order result as it should. In the transition region between the resummation and fixed-order regions, this variation effectively captures the uncertainty in the transition. In the right-hand plot of figure 10 the yellow band enclosed by the dotted green curves shows the effect of this variation on the central profile µ H = m H /4.
A key advantage of the setup we have discussed above is that it provides a concrete way by which to estimate the theoretical uncertainties. Our result for the cross section is obtained via a two-step matching procedure and the variation of the scales at which this matching is performed is a natural way to arrive at a realistic error estimate. To obtain our estimate of the total theoretical uncertainty, we take ∆ FO and ∆ resum as the maximum variation among each of their respective profile variations. We then obtain ∆ tot by adding the two in quadrature, We emphasize that since all variations we perform amount to variations of scales (albeit more intricate than standard scale variations), the resulting perturbative uncertainties ∆ FO , ∆ resum and ∆ tot decrease when increasing the perturbative order of a calculation.

Cross section and power corrections as a function of m b
In this subsection we study the cross section as a function of m b . The reason for this is to confirm that the result obtained in the framework presented in this paper does indeed smoothly interpolate between resummation and fixed-order regions. It also serves as an important validation of the method we employ to estimate uncertainties. In the left-hand plot of figure 11 we show the LO+LL (dashed dark green) and NLO+NLL (solid navy) cross sections including error bands (green and blue bands respectively). We also plot central values for the associated fixed-order cross sections at LO (dotted dark green) and NLO (dashed navy). The right-hand plot of figure 11 displays the relative size of the total LO+LL uncertainty (green band) and of the NLO+NLL resummation (light blue band) and total (navy band) uncertainties. We emphasise that the results in figure 11 are from an implementation of the strict expansion of the cross section in eq. (4.4).
The first feature to point out is that at large m b both LO+LL and NLO+NLL results tend to their fixed-order counterparts (i.e., tend to the LO and NLO cross sections). This clearly shows that the framework we have introduced indeed fulfills this desired property in the limit of large m b . The fact that this transition occurs smoothly is a natural result of the strict expansion used here (that is, there are no higher-order cross terms in our result that might spoil the full cancellation between resummation pieces). The smooth transition between the low and high m b regions is also a direct consequence of the order counting we adopt. If we were not to regard the b-quark PDF as an order α s object, then the strict expansion we have performed would not be possible and a discontinuity would be present in the region of m b ∼ µ H arising from the higher-order terms. Finally, the smooth transition to the fixed-order results indicate that our method for including the power-suppressed O(m 2 b /m 2 H ) terms to the strict EFT result works perfectly, and that as we have discussed in section 2.4 it is indeed the case that (at least to the order we work to) it is possible to consistently include all power-corrections present in the fixed-order result.
Regarding the estimates of the perturbative uncertainty, figure 11 reveals that the error bands we assign are indeed reasonable and robust over the full range of m b , with the NLO+NLL band fully contained within that of the LO+LL. The right-hand plot of figure 11 indicates that the total uncertainty is dominated by the fixed-order scale uncertainty in the large-m b limit, and the resummation uncertainty vanishes, as it should, in the limit µ m → µ H . However, with decreasing m b we see that the resummation uncertainty becomes nonnegligible, forming an important component of the total error. The total uncertainty is of the order of 12-14% over most of the range of m b considered here, and grows as m b is increased beyond the scale µ H .
Our result also allows us to consistently quantify the size of power corrections of . This relies on the observation that in the small m b limit, due to the vanishing of the nonsingular contributions, our result essentially becomes a re-arranged 5F computation. This means that all terms required to obtain a consistently matched prediction can in principle be extracted from a calculation that sets m b = 0 from the outset. The comparison between such a re-arranged 5F calculation and the result where the power corrections are included allows us to study the size of the latter. To illustrate this, in figure 12 we compare the LO+LL prediction where power corrections in m b have been included (solid blue) to  the same prediction made with strictly massless coefficient functions (dotted gray). 12 It is clear that in the small m b limit the nonsingular O(m 2 b /m 2 H ) terms are unimportant and that the matched result can simply be constructed, with negligible errors due to missing power corrections, from massless coefficient functions. This argument indicates that should S-ACOT or FONLL with standard perturbative counting be applied to the case of bbH, then the resulting cross section will likely be almost the same as that of the 5F prediction.
It is also apparent that these power-corrections do increase in importance, their size exceeding 10%, for m b 10 GeV. Therefore, in such parameter regions including them is vital for a faithful description of the cross section. In figure 12 we have also plotted the 5F LO and NLO results, which deviate visibly from the LO+LL result as m b grows, indicating that in such regions a massless 5F prediction becomes an inadequate description of the process.
Finally, we make some brief comments regarding the difference between the cross section obtained under a strict perturbative expansion of eq. (4.4) (or equivalently eq. (2.61)) compared to that obtained by not expanding the two sets of matching coefficients. As mentioned earlier, not performing a strict expansion will generically lead to a discontinuity in the limit µ m → µ H due to the non-cancellation of spurious higher-order interference terms. This is illustrated in figure 13 where we have plotted the NLO+NLL cross section predictions for three hard scale choices under a strict perturbative expansion (solid) and with no strict expansion (dotted), namely using standard PDFs f [5] i,b . At large m b the solid and dotted curves display significant differences and in particular the latter showing discontinuities for µ m → µ H . In the small m b limit however, it is clear that the differences between expanded and nonexpanded approaches become much smaller. In particular, this means that the total uncertainty band in the region of physical b-quark masses is basically the same in the two approaches. This property can be used to greatly simplify the practical implementation in this region and we have exploited this to produce the results in section 4.5. However, it is also important to point out that in general it is only an expanded result, akin to that of eq. (2.49), that guarantees a consistent and smooth matching between resummation and fixed-order regions. These observations may well be important when considering heavy-quark initiated processes where the value of m/Q is not as small as in the setup we study here.

LHC phenomenology
Here we return to the phenomenologically relevant case of a physical bottom mass and consider the cross section as a function of the Higgs mass. The previous subsection con-firmed that the framework we use leads to a cross section that consistently describes both resummation and fixed-order regions. It is of course of interest to compare in a meaningful manner our (N)LO+(N)LL predictions to other predictions available, namely predictions in the 4F and 5F schemes, as well as to the Santander matching prescription used to combine these two. The latter is a practical formula that combines the 4FS and 5FS predictions for the total inclusive cross section through a weighted average of the two [30], This construction is such that the combined result tends to that of the 4FS when the collinear logarithms, ln(m H /m b ) are small, and to that of the 5FS when the logarithms are large (i.e., in the limit m b /m H → 0). The choice of the weight ω is motivated by the fact that this leads to roughly equal weights being assigned for the 4FS and 5FS numbers around m H ∼ 100 GeV, which is the region of 'best' agreement between the 4FS and 5FS predictions (see for example ref. [66]). Beyond this motivation the choice of ω is arbitrary and there is no strong theoretical argument preventing the choice of alternative weights or different ways of averaging the two cross section predictions. Moreover, the practical formula combines two predictions made using different PDF and m b (m b ) inputs, which is somewhat inconsistent. The estimate of the uncertainty on a Santander matched prediction is given by the error band obtained by applying the formula eq. (4.13) to the upper and lower uncertainty curves of the 4F and 5F predictions.
Regarding the set of inputs we use, we have chosen to stick as closely as possible to those used in the LHCHXSWG [67,68,72] and also those used in recent studies of bbH production [65]. Explicitly, we use the default MSTW2008 PDF sets, at the appropriate order for the LO, NLO and NNLO 5FS predictions, whilst we use the fixed-flavour n f = 4 set for the 4FS predictions. For the 4FS and Santander matched results we explore the effect of using m b (m b ) = 4.16 GeV (i.e., a different m b (m b ) than that used in 5FS predictions), as done by the LHCHXSWG. 13 The central choice of hard scale is µ H = (m H + 2m b )/4, with m b = 4.75 GeV, and we vary this hard scale by a factor of two to obtain the fixedorder uncertainty. The errorbars for the 4FS and 5FS predictions are obtained by setting µ F = µ R = µ H , that is we do not consider µ F = µ R variations here. The bands for the Santander matched cross sections are obtained with the Santander prescription.
In figure 14 we plot the 5FS (blue points) and 4FS (green points) cross sections, the  Santander point has been computed with the 4FS component using this increased value of m b (m b ) and is therefore seen to be higher than the point that uses a consistent MS mass in both 4FS and 5FS components. We have checked that the results of figure 14 are fully consistent with those in the literature.
As seen in the previous subsection the error band of the NLO+NLL predictions is contained within that of the LO+LL. The LO+LL uncertainty band is quite wide (and larger than the 4FS or 5FS LO bands) -something that is to be expected from a LO prediction and a maximal variation of the two matching scales, eq. (4.12). Rather than being of concern, this observation gives us confidence that we do not underestimate the inherent uncertainties present and furthermore allows us to trust the size of the NLO+NLL band. We additionally observe that most of the LO+LL scale uncertainty is driven by the µ m variation, while at NLO+NLL the uncertainty band is dominated by the µ H variation. On the other hand, the 4FS displays a less significant reduction in its uncertainty band at NLO, and with its large correction shows relatively poor perturbative convergence due to the presence of unresummed logarithms. The 5FS band shrinks more visibly with increasing order, however the NNLO central value lies outside the NLO band. 14 The (N)LO+(N)LL results lie significantly higher than their respective 4FS (N)LO counterparts. This is due to the resummation contained in the former results. We also notice that the NLO+NLL results lie slightly above the NNLO 5FS predictions. The reason for this is that the former results contain the (positive) effects of light channels (gg, qg and qq) at O(α 3 s ) whilst the latter results contain (negative) two-loop corrections to the bb channel (see figure 8).
By construction the Santander matched result lies between the 4FS NLO and 5FS NNLO predictions (irrespective of the precise inputs for m b (m b ) used in the 4FS calculation). Our NLO+NLL result is therefore higher than the Santander matched results, and specifically we find an increase of 6% with respect to the magenta LHCHXSWG Santander point, and of 12% with respect to the brown Santander point (which is directly comparable to the NLO+NLL result). We also notice that the error bands of the NLO+NLL prediction are roughly of the same size as those obtained through Santander matching indicating that the size of the latter is likely to be realistic. There is a sizeable overlap in the uncertainty bands of the Santander and the NLO+NLL predictions, however the central Santander points lie towards the bottom edge (magenta) and outside (brown) of the NLO+NLL error band. We also point out that while our NLO+NLL result is stable upon hard scale variation, the Santander matched result would be significantly smaller for larger hard scales, as is clear by inspection of figure 9.
In figure 15 we present the same comparison as in figure 14 for larger Higgs masses, m H = 300 GeV and m H = 500 GeV. The overall pattern observed for a light Higgs remains qualitatively unchanged. The NLO corrections in the 4FS have grown noticeably, such that the 4FS NLO results are now outside the 4FS LO bands, which is likely due to the unresummed logarithms, which get larger at higher Higgs masses. While at lower m H the central values of the 4FS NLO and 5FS NNLO are relatively close to each other and have overlapping uncertainty bands, at higher m H they are further apart and have only barely overlapping uncertainty bands. Consequently, the Santander average becomes even less reliable. Encouragingly, the LO+LL and NLO+NLL results at large m H continue to display the good perturbative behaviour and convergence pattern present at lower m H values. They also remain systematically higher than the Santander matched predictions.
Finally, we briefly comment on the effect of including the known C (2) bb term 15 (i.e., the two-loop corrections to the bb-channel), which is formally of higher order in our approach. This is the only known contribution that we have not included in the NLO+NLL result. The NLO+NLL result for m H = 125 GeV shown in figure 14 is σ NLO+NLL = 0.224 ± 0.021 pb. With the addition of the higher-order C (2) bb -term this cross section becomes σ NLO+NLL+C (2) bb = 0.211 ± 0.010 pb, to be compared with the NNLO 5FS cross section of σ 5F,NNLO = 0.209 ± 0.010 pb. Clearly the addition of this higher-order term reduces the NLO+NLL cross section and additionally reduces its uncertainty band, bringing both central value as well as the size of the error bands closer to those of the 5FS NNLO this in order to directly compare to our method of estimating uncertainties. 15 Additionally, we include the lowest order bb-channel term C (2) bb , appearing at the same order, which however gives a negligible effect. result. This term is likely to be important at very high scales µ H 1 TeV, where the 5FS perturbative counting is expected to be more appropriate. However, we feel that including this term for a SM Higgs is rather ad-hoc given that there are a number of additional terms that would also contribute at the same order, but which are not known and have not been included here. Furthermore, the sizeable reduction of the error bars is slightly discomforting given that we have no control on the effects of these missing terms. We note that the error bars of the original NLO+NLL cross section nicely cover the effect of adding the C (2) bb term, which provides further support that the uncertainty bands presented are indeed reasonable.

Conclusions
We have presented a systematic EFT setup to derive heavy-quark initiated cross sections at hadron colliders. Our framework includes the resummation of potentially large logarithms ∼ ln(m b /Q). Furthermore, it consistently includes power corrections ∼ m 2 b /Q 2 , reproducing the full fixed-order (4FS) result. As such our final result gives predictions that are accurate in both of the limits m b Q and m b ∼ Q as well as in the transition region in between.
Our result is obtained via a two-step matching procedure, and variation of the scales at which these two matchings are performed allows us to obtain a robust estimate of the perturbative uncertainties. The construction of the coefficient functions of our result bears several similarities with existing VFNSs for DIS. A key difference in our approach is the different perturbative order counting. In particular, we argue that it is more appropriate to count the effective b-quark PDF, which is generated perturbatively at the scale µ m , as a perturbative object of O(α s ). This organization of the perturbative series leads to perturbatively stable results and allows for a smooth transition between fixed-order and resummation regions. The simplicity of the EFT approach for DIS makes its generalization to hadron-hadron collisions straightforward.
We have applied our framework to the case of the bbH cross section. We first studied the cross section as a function of m b , which served to demonstrate that our resummed and matched result satisfies all required properties. We then presented numerical results of phenomenological interest for the LHC. We have compared our results to the 4FS and 5FS result, as well as the Santander prescription, which combines the two results by taking a weighted average. Since our predictions are derived from a field-theory setup, consistently combining the 4FS and 5FS limits, it can be regarded as a definite improvement over this prescription. At NLO+NLL, we find a slightly reduced perturbative uncertainty compared to the Santander average. The Santander central value is lower than the NLO+NLL result by about 12% for m H = 125 GeV and lies outside our uncertainty band when consistent MS masses are used in both the 5FS and 4FS ingredients of the Santander result. The difference is reduced to 6% when using a larger MS mass in the 4FS result, as used by the LHCHXSWG, with the central value lying at the lower edge of our uncertainty band. We observe that the NLO+NLL result is stable upon variation of the hard scale, unlike the Santander matched results which would vary significantly under such variation.
The framework presented in this paper can be straightforwardly applied to other processes involving initial-state b-quarks, for example single-top or V + b-jet production. Additionally, extending the framework to study more differential observables of interest, such as the transverse momentum of the Higgs in bbH-production or jet-vetoed cross sections is possible. In making the calculation less inclusive, more scales appear in the perturbative expansion, which can be efficiently dealt with in an EFT setup. Finally, it would be very interesting to adopt our perturbative counting in DIS in the context of PDF fits, where the used variable flavor number scheme typically plays a central role in determining the accuracy and the goodness of the fit itself.
Note added: While this paper was being finalized ref. [74] appeared, which obtains the bbH cross section in the FONLL approach. As discussed in section 3, the FONLL approach adopts a different perturbative counting to the one we have presented here. In particular, the result of ref. [74] is computed at FONLL-A accuracy, that is, it combines the 4FS LO result with the 5FS NNLO result. In comparison to our NLO+NLL result, it does not include the 4FS NLO contributions fromC (3) gg ,C qq , andC (3) qg . However, it does include the C (2) bb and C (2) bb 5FS NNLO terms, which are higher-order terms in our approach (and whose impact on our result is discussed at the end of section 4.5). As a result, the FONLL-A result of ref. [74] is very close to the 5FS NNLO, as one might expect, since the difference to the latter is the inclusion of the numerically small O(α 2 s m 2 b /m 2 H ) power corrections from the 4FS LO result. Finally, ref. [74] does not include an estimate of the resummation and matching uncertainty (analogous to our µ m variation).

A Renormalization
In this appendix we briefly discuss aspects of renormalization important to our work. In the MS scheme ultraviolet (UV) divergences associated with quark lines are subtracted in the same way for all quarks, independently of their mass (mass doesn't play a role in the UV). However, a direct consequence of the MS scheme is that α s runs with n f = 6 flavors, irrespectively of the energy scale. A more physical renormalization scheme for heavy quarks is the so called CWZ or decoupling scheme [75], where all "light" quarks are renormalised with MS counter terms, while UV divergences associated with the heavy-quark loops are subtracted at zero momentum. This ensures decoupling of the heavy quarks at energies much smaller than their mass. Since the concept of light and heavy depends on the actual scale, in practice a variable flavor number renormalization scheme is used, where the number of active (light) quarks renormalised in MS depends on the hard scale and changes at the crossing of the heavy quark thresholds. We have used this in Sect. 2 when treating the fixed-order and resummation regions differently regarding the renormalization of b-quark loops.
At the threshold scale µ = µ m , matching conditions relate the value of α s above and below that scale. Denoting with a superscript the number of flavors used in the evolution of α s , we have α [5] s (µ 2 m ) α [4] s (µ 2 m ) = 1 + α [4] s (µ 2 m ) π where m is the heavy quark pole mass. These conditions are currently known through 4 loops [76], but for our applications we just need the 2 loop expression of eq. (A.1). A generic observable F (α s ) can be written as a perturbative expansion equivalently in both schemes, F (α s ) = k α [4] s k (µ)F [4](k) (µ) = k α [5] s k (µ)F [5](k) (µ), (A.2) and to all orders they are identical. The relation between the coefficients F [4](k) , whose heavy quark UV divergences are renormalised in the decoupling scheme, and F [5](k) , renormalised in MS, can be simply obtained by using eq. (A.1) to write α [4] s in terms of α [5] s in the first sum (or viceversa), re-expanding and matching order by order.

B Construction of PDFs with variable threshold
In this work we have considered PDFs in which the heavy quark thresholds are not fixed to the heavy quark mass, but can vary. Moreover, the effective PDFs defined in eq. (2.48) are required with mixed evolution and matching accuracies. If we consider the "universal" PDFs as those at a small scale µ Λ , then the dependence on the threshold and the details of the order at which each ingredient is retained are all in the perturbative evolution. Typically PDFs are used through the LHAPDF library, where evolved PDFs at any (available) scale are built from interpolation grids, previously created assuming a particular evolution with specific heavy quark thresholds.
For our purposes, performing the evolution each time picking the desired value of the heavy quark thresholds seems advisable. However, this approach faces speed problems, since performing the evolution is much more time consuming than interpolating a grid. Therefore, for the work in this paper we have created LHAPDF grids with different choices for the b-quark threshold, µ m , and with the required combinations of the perturbative ingredients at different orders. To do so, we have used the public code APFEL, which has the ability of creating grids after performing its own evolution. To accommodate the possibility of choosing a threshold different from the heavy quark mass, we have modified the code adding the µ m -dependent matching conditions through NNLO from Ref. [53]. Furthermore, we modified the code to produce PDF grids with the required combinations of orders of the matching coefficients, M ij . Practically, we proceeded as follows: • We start with a central member of a public PDF set, namely MSTW2008 with α s (m Z ) = 0.1171, both at NLO and at NNLO. The value of the bottom pole mass used in this work has been taken to be m b = 4.75 GeV for consistency with the chosen PDF set.
• We compute all the required PDFs as well as α s from this set at an initial scale µ 0 = m c = 1.4 GeV.
• We use APFEL to perform the forward evolution and create the corresponding grids for the different setups we are interested in.
Note that we choose to perform the evolution at (N)NLL, starting from a (N)NLO set, for producing the PDFs that we used in our (N)LO+(N)LL results. While this is somewhat in constrast with the discussion in section 2.3, using the evolution at one higher order has two advantages. The first is that the order of the evolution and the highest order in the matching functions are consistent, so that when we do not use the strict expansion we can use standard PDFs, as explained below eq. (2.50). The second is that our final NLO+NLL result in section 4 is more directly comparable to the 5FS NNLO result. As an illustration of what the modified code is capable of, in figure 16 we show the standard 5F bottom PDF f [5] b for a fixed value of x = 0.005. We plot this as a function of the scale µ close to the bottom threshold for three different values of the threshold itself, µ m = m b /2, m b , 2m b , as used in the phenomenology section 4.5. Changing the threshold of course shifts the value of the PDF at µ m , and in particular makes this value nonzero at NLO (at NNLO it is already nonzero even for µ m = m b ). The initial condition, shown as the gray dotted or dashed lines, is given by the product of 4F PDFs with the matching conditions. We observe that at NNLO the PDFs obtained with different values of the threshold are almost identical at a large scale, indicating that a NNLO cross section is only likely to have a mild dependence on µ m . We also note that the initial condition itself behaves in a nice perturbative way at large scales, where we expect higher order corrections to be small, while it deviates strongly at smaller scales, where α s is larger and higher-order corrections are not negligible.

C Massive kinematics in hadron-hadron collisions
In the case of an incoming massless parton its four-momentum is considered to be a fraction of the four-momentum of the proton. In case of massive incoming parton this formulation would lead to a mass that scales with the momentum fraction, which is clearly inconsistent. To overcome this, the proper approach [5] is to use light-cone coordinates, where momenta can be written as p = p + , p − , p t , p ± = (p 0 ± p 3 )/ √ 2, (C.1) where p i are the usual Minkowski components. We choose to orient the beam axis along the third spatial direction. A collinear particle with mass m has p t = 0, and the mass can be expressed as m 2 = 2p + p − . We can write the momenta of massless protons 16 in the 16 The proton mass can always be neglected compared to its momentum at the LHC.
When there is just one proton, as in DIS, the problematic configuration described above never takes place, the effect of the parton's mass being a further restriction to the accessible values of x with respect to the massless case, and standard collinear factorisation works even in presence of massive partons [5]. In the hadron-hadron collider case, only a systematic expansion in the heavy quark mass such as the one presented in this work allows the description of heavy quarks in the initial state.