Matched predictions for the $b\bar{b}H$ cross section at the 13 TeV LHC

We present up-to-date matched predictions for the $b\bar{b}H$ inclusive cross section at the LHC at $\sqrt{s}=13$ TeV. Using a previously developed method, our predictions consistently combine the complete NLO contributions that are present in the 4-flavor scheme calculation, including finite b-quark mass effects as well as top-loop induced $Y_b Y_t$ interference contributions, with the resummation of collinear logarithms of $m_b/m_H$ as present in the 5-flavor scheme calculation up to NNLO. We provide a detailed estimate of the perturbative uncertainties of the matched result by examining its dependence on the factorization and renormalization scales, the scale of the Yukawa coupling, and also the low b-quark matching scale in the PDFs. We motivate the use of a central renormalization scale of $m_H$/2, which is halfway between the values typically chosen in the 4-flavor and 5-flavor scheme calculations. We evaluate the parametric uncertainties due to the PDFs and the b-quark mass, and in particular discuss how to systematically disentangle the parametric $m_b$ dependence and the unphysical b-quark matching scale dependence. Our best prediction for the $b\bar{b} H$ production cross section in the Standard Model at 13 TeV and for $m_H$ = 125 GeV is $\sigma(b\bar b H) = 0.52 \,{\rm pb}\, [1 \pm 9.6\% {\rm (perturbative)}\,{}^{+2.9\%}_{-3.6\%} {\rm (parametric)}]$. We also provide predictions for a range of Higgs masses $m_H\in [50, 750]$ GeV. Our method to compute the matched prediction and to evaluate its uncertainty can be readily applied to other heavy-quark-initiated processes at the LHC.


Introduction
Heavy-quark initiated processes and accurate theoretical predictions for them are important for precision tests of the Standard Model (SM), searches for New Physics, and are also essential for the precise determination of parton distribution functions (PDFs).
Predictions for processes involving initial-state b-quarks at hadron colliders are typically made using factorization theorems that are valid for two different parametric hierarchies between the b-quark mass, m b , and the hard-interaction scale of the process, Q. In the 4-flavor scheme (4FS), one formally considers m b ∼ Q, while in the 5-flavor scheme (5FS) one formally considers m b Q. The predictions obtained in either scheme have different features. The 4FS predictions include power corrections in m b /Q as well as the exact massive quark final-state phase space, while logarithms of ln (m b /Q) are included at a given fixed order. The 5FS predictions do not include power corrections in m b /Q but resum the potentially large logarithms ln (m b /Q) to all orders via DGLAP into a b-quark PDF. The construction of matched predictions that include the merits of both schemes in DIS (sometimes called variable flavor number schemes) has been the subject of many years of work [1][2][3][4][5][6][7][8][9][10][11][12][13] and recently progress in this direction for hadron-hadron colliders has also been made [14][15][16].
The focus of this paper is on bbH production, for which predictions for the inclusive cross section are available up to NNLO in the 5FS [17][18][19][20][21] and NLO in the 4FS [22][23][24]. In the past, the results in the two schemes have been averaged using the pragmatic Santander-Matching prescription [25]. Given that this prescription amounts to a simple weighted average of the two results, it is clear that it does not constitute a satisfactory and theoretically-consistent matching. In particular, as already observed in ref. [16], the contributions present in one and not the other scheme should get added in their combination rather than averaged, which in this case leads to a noticeably higher cross section of the properly matched result compared to the Santander average.
In ref. [16], we used a simple effective field theory setup to systematically derive matched predictions for the bbH cross section, which combines the ingredients of both 4FS and 5FS schemes. The result contains the full 4FS result, including the exact dependence on the b-quark mass, and improves it with a resummation of collinear logarithms of m b /m H , as are present in the 5FS. An important aspect of our construction is the perturbative counting of the effective b-quark PDF, which is counted as an O(α s ) object for phenomenologically relevant hard (factorization) scales below ∼ 1 TeV. This counting rearranges the perturbative expansion of the cross section, such that in the massless limit the result corresponds to a reorganized 5FS result. An important advantage is that in this way the logarithms present in the resummed (5FS) result order-by-order exactly match those present in the fixed-order (4FS) contributions. This feature greatly facilitates the combination of both types of contributions.
In this paper, using the approach of ref. [16], we present state-of-the-art predictions for the bbH cross section for the LHC at 13 TeV, including a comprehensive study of its perturbative and parametric uncertainties. Our predictions include the effects of top-quark loop-induced interferences, proportional to Y b Y t , which are known to be important in the Standard Model [22][23][24]. (These were not yet included in ref. [16].) We also investigate in detail the effects of pure 2-loop terms that are present in the 5FS NNLO calculation but are formally of higher order in our perturbative counting.
The two-step matching used in ref. [16] makes it explicit that there are two relevant scales in the problem: the hard scale µ H ∼ m H and the b-quark scale µ b ∼ m b (also referred to as the b-quark threshold scale). Since a perturbative expansion is performed at both of these scales, a reliable theory uncertainty should take into account the perturbative uncertainties related to both. The uncertainty due to µ b has been neglected in the past in essentially all 5FS and matched cross section predictions, since all standard 5FS PDF sets make the fixed choice µ b = m b . However, the µ b uncertainty should be regarded as an additional resummation uncertainty, and in ref. [16] it was shown how to systematically estimate it and that it can indeed have a nonnegligible effect. We discuss and motivate an appropriate choice of the central scales and provide a full breakdown of theoretical uncertainties due to µ F , µ R , µ b variations.
We show that this more general setup also disentangles the dependence on µ b and m b and thus allows one to correctly evaluate the uncertainty in the predictions due to the parametric uncertainty in the value of m b . This discussion is relevant for any heavy-quark initiated process calculated in either the fixed 5FS or a matched approach.
The structure of the paper is as follows: in section 2, we recall the main features in the construction of the matched NLO+NLL result and discuss the extensions to include the Y b Y t interference terms and higher-order 2-loop terms. In section 3, we discuss in detail the perturbative and parametric uncertainties. In section 4, we present our numerical results for the inclusive cross section for a range of Higgs masses m H ∈ [50, 750] GeV, including the full breakdown of all uncertainties. We conclude and give a short outlook in section 5.

Matched cross section
In this section we discuss the theoretical ingredients for our predictions of the bbH cross section. In sections 2.1 and 2.2, we summarize the main steps and results of ref. [16] in the construction of the fixed-order + resummation matched result, which is valid for any parametric hierarchy between m b and m H . For the detailed derivation we refer to ref. [16]. In section 2.3, we discuss the inclusion of the formally higher-order two-loop terms in the resummed results. We will denote the corresponding result as NLO+NNLL partial . In section 2.4, we discuss the inclusion of the Y b Y t interference terms.

Fixed order and resummed results
The fixed-order cross section that enters the matched result, and which is recovered when the resummation is turned off, is obtained in a factorization scheme that is formally valid for the parametric power counting m b ∼ m H . Hence, it is equivalent to a 4FS result, where bottom quarks do not appear in the initial-state and can only be produced via gluon splittings into b-quark pairs. It includes the exact dependence on the b-quark mass (i.e., it includes power corrections to all orders in m 2 b /m 2 H ) both in the partonic cross sections and in the phase space. Logarithmic terms ∼ ln(m 2 b /m 2 H ) arising from collinear gluon splittings are included at fixed order in α s . The fixed-order result is given in terms of coefficient functions D ij (m H , m b , µ F ), which depend explicitly on m b , and 4-flavor PDFs, ik (µ F , µ 0 )f [4] k (µ 0 ) U [4] jl (µ F , µ 0 )f [4] l (µ 0 ) . (2.1) For notational simplicity, we keep all Mellin convolutions between PDFs, evolution factors, and coefficient functions implicit and simply write products throughout the paper. The sum over partons only includes gluons and light (anti)quarks (q = d, u, s, c), µ F is the hard factorization scale of the process. In the second line the 4-flavor PDFs f [4] j (µ F ) are written in terms of the fitted PDFs at the low scale µ 0 1 and DGLAP evolution factors U [4] ij with n f = 4 active quark flavors. The fixed-order result of eq. (2.1) is equivalent to a 4FS result and the coefficients D ij for bbH are known to NLO.
The pure resummed result is based on the factorization theorem derived in the limit m b m H , i.e., in a power expansion in m b /m H , where the leading power term leads to the usual 5FS result, while subleading power terms correspond to subleading twist terms and are usually not considered. In this case, b-quarks are treated as massless at the hard matching scale µ F ∼ m H and appear in the initial state of the corresponding partonic matching calculation. The collinear logarithms of m b /m H are resummed by DGLAP evolution from the hard scale µ F ∼ m H to the b-quark matching scale µ b ∼ m b . At µ b , the b-quark is integrated out of the theory. These steps result in an effective perturbative bquark PDF. The resummed cross section is given by the convolution of coefficient functions C ij (m H , µ F ), which contain no m b dependence, with 5-flavor PDFs, In the second step, the 5-flavor PDFs at µ = µ F are written out explicitly in terms of the 5-flavor evolution from µ F to µ b , the matching at µ b yielding matching coefficients M kl (m b , µ b ) that explicitly depend on m b , followed by the 4-flavor evolution from µ b to µ 0 and the fitted 4-flavor PDFs at scale µ 0 . All mainstream 5FS PDF sets construct their b-quark PDFs in this way, however, with the fixed choice of µ b = m b . With this choice, the O(α s ) matching coefficients M ij in the m b pole-scheme are exactly zero, which somewhat simplifies the implementation. However, identifying µ b = m b confuses the parametric physical dependence on m b and the unphysical dependence on the matching scale µ b , which controls the resummation of logarithms and should cancel to the order one is working. We will discuss how we rectify this situation in section 3.

Matching fixed order and resummation: NLO+NLL
In all practical applications we are aware of, the evolved PDFs are always treated as external O(1) objects and the perturbative expansion of the cross section is based solely on the perturbative expansion of the coefficient functions D ij and C ij in eqs. (2.1) and (2.2). For the fixed-order and resummed results this leads to the usual 4FS and 5FS predictions. The corresponding contributions are schematically depicted in figure 1 by the first column for the 4FS (green dotted boxes) and by the three diagonals (blue shading) for the 5FS. As noted above, the b-quark PDF is itself a perturbative object with the expansion In principle, each of its terms should be included in the perturbative counting and be regarded as part of the perturbative expansion of the cross section. Counting it as an O(1) object would be justified in the limit where the off-diagonal evolution factor U [5] bg (µ F , µ b ) ∼ 1. However, U [5] bg (µ F , µ b ) is suppressed by an overall α s ln(µ F /µ b ) relative to the diagonal evolution factors and vanishes in the limit µ b → µ F , and therefore this only holds for scales µ F ≫ µ b . Numerically, for µ b ∼ m b this is only attained for scales µ F 1 TeV. Hence, for the scales of interest here it is more appropriate to count U [5] bg (µ F , µ b ) as O(α s ), as indicated in eq. (2.3), in which case the whole f becomes an O(α s ) object. The perturbative expansion in α s of the resummed cross section in eq. (2.2) is then performed by expanding the product of coefficient functions C ij together with the terms making up the b-quark PDF including the b-quark matching coefficients M ij and U [5] bg ∼ α s . For a more detailed discussion we refer to ref. [16].
The 4FS and 5FS results can significantly differ from each other and in particular display different patterns in their factorization-scale dependence due to the different logarithmic terms present at each order in the two schemes. Hence, a consistent matching appears to be nontrivial. The above treatment of the resummed result has the important added advantage that it reorganizes the resummed series into a form that is consistent with the logarithms present in the fixed-order result. The key feature is that order by order in α s the limit µ b → µ F in the resummed cross section now exactly reproduces all the logarithmic terms (and nothing more) that are present in the m b → 0 limit of the fixed-order cross section. In other words, the reexpansion of the resummed result to fixed order is simply given by setting µ b = µ F . This in turn means that for µ b < µ F the evolution factors U [5] ij in this expansion precisely resum the singular logarithms present in the fixed-order result. Hence, all that is missing in the resummed result compared to the fixed-order result are purely nonsingular terms proportional to m 2 b /m 2 H , i.e. terms that vanish in the limit m b → 0 given by The complete matched cross section is then simply given by adding the nonsingular fixedorder terms to the resummed result, By construction, it satisfies σ FO+resum → σ FO in the limit µ b → µ F where the resummation is turned off, as required for a consistently matched prediction. On the other hand, it reduces to σ resum in the limit m H m b . We emphasize again that this crucially relies on the fact that the nonsingular terms vanish for m b → 0, which in turn relies on adopting the perturbative counting for the resummed result described above. 2 The corresponding terms included in the matched result at each order are depicted by the rows (red boxes) in figure 1.
As discussed in ref. [16], for the practical implementation, the nonsingular contributions can be conveniently absorbed into modified gluon and light-quark coefficient functions, C ij (m H , m b , µ F ), which now carry an explicit dependence on m b , convolved with effective 5F PDFs. 3 The final matched result is then written as where f [5] i,b are perturbative objects, and an expansion of C ij andC ij against f [5] i,b as discussed above is implicit.
The strict expansion of the coefficient functions against the individual terms making up the b-quark PDF is quite inconvenient for the practical implementation, as it requires performing the entire n f = 5 DGLAP evolution above µ b by hand. However, as long as we are only interested in the phenomenologically relevant region µ b ∼ m b µ F ∼ m H , we can also keep formally higher-order cross terms in order to simplify the practical implementation. Doing so allows us to use common preevolved 5FS PDFs, under the condition that we count f [5] b ∼ O(α s ) while the light quark and gluon PDFs are counted as ∼ O(1). In addition, we have to use PDFs of sufficiently high order such that they include all matching corrections required by our perturbative counting. Specifically, at (N)LO+(N)LL this requires the use of at least (N)NLO PDFs. It was explicitly checked in ref. [16] that for values of m b /m H 0.1 this implementation gives practically the same numerical results as the strict expansion. On the other hand, the strict expansion is required if one wishes to explicitly study the limit µ b → µ F and obtain a smooth transition of the matched result into the fixed-order result.
In summary, with this simplification, expanding the matched cross section in powers of α s = α s (µ F ), the following perturbative expansion is obtained The factors of two and four account for the exchange of partons among the two protons and (to a first approximation) the equality f . A sum over light quarks and gluons is implicitly assumed for repeated indices i, j, k. The superscripts on the coefficient functions indicate the order in α s to which these are computed. The first two orders in eq. (2.7) are illustrated by the red boxes in figure 1. As seen there, our perturbative counting implies that we include bb, bg and gg initiated contribution consistently at the same order.
Finally, we note that the construction of the coefficient functionsC ij is formally the same as the corresponding construction in the FONLL approach [15] (and in a hypothetical S-ACOT construction). There are, however, two main differences between these approaches. First, as explained above, we use the fact that the effective b-quark PDF counts as an O(α s ) perturbative object to construct the perturbative expansion of our matched result. As a result, it contains the complete fixed-order result at each perturbative order, and (with the strict expansion) smoothly merges into it. Secondly, as discussed further in section 3, in our approach we explicitly distinguish µ b and m b allowing us to include an explicit estimate of the resummation uncertainty associated with the 5F resummation by varying the (in principle arbitrary) matching scale µ b .

Higher-order two-loop terms: NLO+NNLL partial
At present, all coefficient functions in eq. (2.7) required up to NLO+NLL are known [19,20,22,23]. Going to NNLO+NNLL is not yet possible and would require the full NNLO 4FS result, corresponding to the unknownC (4) ij and C (3) bk coefficients in eq. (2.7), together with the three-loop matching coefficients M ij . The two-loop coefficients C (2) bb and C (2) bb are known from the NNLO 5FS result [19] and are part of the NNLL resummation in our counting. Adding them to the NLO+NLL result provides a partial NNLL result, see figure 1, which we denote as NLO+NNLL partial . (In ref. [16], this result was called NLO+NLL+C (2) bb . It corresponds to what in ref. [15] would be called FONLL-B.) Including the partial NNLL terms violates the exact correspondence between the resummed and fixed-order results. That is, the fixed-order limit of these resummed terms are only reproduced by the m b → 0 limit of the fixed-order result at NNLO. Hence, in the limit µ b → µ F these terms spoil the smooth matching of the matched result into the fixed-order result. In this regard, these terms are analogous to the higher-order cross terms that are kept when the strict expansion is not implemented. Therefore, whenever the strict expansion is used these terms should also be consistently dropped. This is the case when going toward larger values of m b /m H , where the logarithms become small and fixed-order contributions become dominant, and the matched result is transitioning into the pure fixed-order result.
For intermediate values of m H (small values of m b /m H ) where our perturbative counting is most appropriate, including these terms does not necessarily improve the overall accuracy of the prediction, since other terms of the same perturbative order are not included and including only a partial set of terms might bias the result in the wrong direction. At the same time, their inclusion can lead to a reduced scale dependence, which would then potentially underestimate the perturbative uncertainties. For these reasons, we do not take the NLO+NNLL partial prediction to be our default result. Nevertheless, we also provide it in section 4, as it can provide an indication of the numerical size of the next-order correction. This can for example be an additional useful cross check on the estimate of the perturbative uncertainties of the NLO+NLL result or alternatively guide the choice of the central scales.
Finally, in the limit of very large m H (very small values of m b /m H ), including these terms becomes beneficial once the size of the resummed logarithms grows, α s ln(µ 2 F /µ 2 b ) ∼ 1, and the original strict 5FS counting applies.

Fixed
Starting at NLO, the cross section receives contributions proportional to Y b Y t due to the interference of the Born-level gg → bbH diagrams with diagrams where the Higgs is radiated from a closed top-quark loop; some examples of the latter diagrams are shown in figure 2. The fixed-order (4FS) cross section can be written schematically as where the interference terms are included in σ correspond to an interference between the bbH and gluon-fusion processes. Since they involve b-quarks in the final state, they are usually regarded as part of the bbH cross section.
Curiously, these terms can be treated as purely nonsingular terms, even though they may, at first sight, appear to contain large logarithms of m b /m H in the m b → 0 limit. On closer inspection, these interference terms turn out to vanish to all orders in the m b → 0. The reason is that in the interference between bbH-like and gluon-fusion-like diagrams the Higgs boson must be attached to two different closed fermion lines. This requires a helicity flip on the b-quark line, which is not allowed for m b = 0. Equivalently, for m b = 0, they will always contains a trace over an odd number of Dirac matrices and thus vanish [19]. For the same reason, such terms are absent in the 5FS.
For us, this means that these terms are purely nonsingular and can be straightforwardly added by including them in σ FO in eq. (2.8), which then enters into theC (3) ij in eq. (2.7). In practice, we extract the numerical result for σ and σ (1) Y b Yt in the pole scheme from Madgraph5 aMC@NLO [26] by generating the process pp → bbH at NLO with Y t turned on and off. These are then used to constructC (3) ij in the MS scheme for the Yukawa couplings. The Y b Y t interference terms have a noticeable numerical effect (∼ 5%) in the SM, while in beyond-the-Standard-Model (BSM) scenarios such as SUSY with large tan β their relative effect compared to the dominant Y 2 b contribution tends to be much milder. In section 4 we therefore provide the results both with and without the Y b Y t terms included, which we denote as NLO[Y 2 b ] and NLO[Y 2 b +Y b Yt] in the following. For our choice of central scales, the Y b Y t terms reduce the cross section for m H 300 GeV and increase it for m H 300 GeV, see figure 4 below.

Estimate of perturbative and parametric uncertainties
We now turn to the discussion of the theoretical and parametric uncertainties. The estimation of the perturbative uncertainties by variations of the hard scales µ F and µ R , and low matching scale µ b are discussed in section 3.1. The parametric uncertainty from the value of the b-quark mass is discussed in section 3.2. In section 3.3, we discuss the PDF uncertainty and the construction of modified PDF sets that are required to properly separate the m b and µ b uncertainties.

Scale choices and perturbative uncertainties
As discussed in detail in ref. [16], we can distinguish two different sources of perturbative uncertainties. One is an overall "fixed-order" uncertainty (within the resummed or matched results), which can be estimated by exploiting the dependence on the hard matching scale. The second is a "resummation" uncertainty related to the uncertainty in the resummed logarithmic series, which can be estimated by exploiting the dependence on the low µ b matching scale.
In ref. [16] we considered a common hard scale. Here, we additionally study the dependence of the cross section on the renormalization scale µ R at which α s is evaluated and on the renormalization scale µ Y at which the Yukawa coupling is evaluated. In this case, the role of the hard matching scale in the resummation is played by the factorization scale µ F at which the PDFs are evaluated.

Central scale choices
For the factorization scale µ F we use the central choice µ F = (m H +2m b )/4, as is commonly used in both the 4FS and 5FS calculations. This choice is motivated by the well-known observation that in bbH such a small factorization scale leads to an improved perturbative convergence, see e.g. [16,[27][28][29][30]. We point out that the matched NLO+NLL result turns out to be significantly less sensitive to the central value of For the renormalization scale we use a somewhat larger central value µ R = m H /2. This is motivated by the fact that kinematical arguments for a small scale ∼ m H /4 are related to the collinear factorization (µ F ) and not the renormalization (µ R ). On the other hand, choosing µ R = m H , which would be the canonical renormalization scale, produces somewhat artificial leftover ln(µ R /µ F ) ln 4 terms in the cross section, which become even larger under scale variations. The value µ R = m H /2 is a reasonable compromise that lies halfway between the standard 4FS and 5FS choices of µ R = (m H + 2m b )/4 and µ R = m H . Also, at the Higgs masses of interest, the matched result is dominated by the resummation contributions and as shown in ref. [30] the 5FS tends to favor µ R values between m H /2 and m H . Finally, this choice has the convenient side effect that the NLO+NNLL partial result turns out to be a very small correction over the NLO+NLL result.
The Yukawa couplings Y b (µ Y ) and Y t (µ Y ) are defined in the MS scheme are obtained by evolving from m b (m b ) = 4.18 GeV [31] and m t (m t ) = 162.7 GeV [32] to the central Yukawa scale µ Y with 4-loop evolution, while µ Y variations are computed using 2-loop evolution. While both µ R and µ Y are renormalization scales, they do not necessarily need to be the same. It is always possible to evolve α s and the Yukawa coupling to different scales using their own renormalization group evolution, compensated by including the appropriate fixed-order logarithms in the partonic coefficients. In figure 3 we study the dependence of the cross section on µ R and µ Y at different orders, always keeping µ F fixed at its central value. We find that varying µ R and µ Y together gives the largest scale variation, so for our numerical results and uncertainty estimation we identify µ Y ≡ µ R as is usually done.  We also observe that at LO+LL and NLO[Y 2 b ]+NLL the renormalization scale dependence comes entirely from the µ Y dependence, which reduces significantly at each higher order. (This is also another motivation to choose a higher central scale for µ R since m H is the (only) relevant scale seen by the bbH vertex.) Note that at NLO[Y 2 b ]+NNLL partial the µ Y dependence reduces, which is precisely to be expected from the included two-loop virtual corrections. On the other hand, the µ R dependence for fixed µ Y at NLO[Y 2 b ]+NNLL partial actually increases, which could well be related to only including a partial set of higher-order terms.
We note that this observed pattern for the µ dependence is somewhat changed by the inclusion of the Y b Y t interference terms. This is not unexpected, since these terms introduce new LO dependence on both µ R and µ Y . However, since their absolute correction to the cross section is small, we base our discussion of the scale choices on the pattern observed without interference terms.
Finally, for µ b we take the canonical central value µ b = 4.58 GeV, which corresponds to the central pole mass value we use, but importantly is kept fixed under m b variations. The canonical scale choice µ b = m b is appropriate in the resummation region where µ b µ F , which is the case for all Higgs masses we consider in this paper. As explained in ref. [16], for larger values of m b /µ F 0.3 (smaller m H ) one enters the transition region, where the strict expansion should be used and µ b should be chosen via a more general profile scale to allow for a smooth turning off the resummation and transition into the fixed-order limit m b ∼ µ F .

Estimate of perturbative uncertainties
In ref. [16], µ R = µ Y = µ F were all varied together, which already yields an excellent perturbative convergence. Here we follow an even more conservative approach and also explore independent variations of µ R and µ F . We consider the usual 7-point variation where the two scales are varied independently up and down by factors of two, excluding the two cases where they are both varied in opposite directions. Whenever varying µ F up or down, the low matching scale µ b is varied up or down by the same factor, such that the ratio µ F /µ b and therefore the resummed logarithms remain fixed. As discussed in ref. [16], this allows us to interpret the hard scale variations as an estimate of the overall fixed-order uncertainty. The final fixed-order uncertainty is then obtained by the maximal envelope, i.e., we use the absolute value of the largest deviation from the central value as the symmetric uncertainty. Doing so explicitly avoids attributing any physical meaning to accidentally small one-sided scale dependence (that would yield asymmetric uncertainties), which just results from a nonlinear scale dependence as frequently encountered at higher orders or near the points of minimal scale dependence. Next, the resummation uncertainty is computed by varying the matching scale µ b by a factor of 2 up and down about its central value. For this variation we keep all the other scales fixed (and also the value of m b ). Doing so explicitly changes the ratio µ F /µ b and thus directly probes the size of the resummed logarithms.
The fixed-order and resummation uncertainties are considered as independent uncertainty sources, and the total perturbative uncertainty is obtained by adding them in quadrature. A simple alternative approach, which however lacks the physical interpretation of the source of uncertainty, would be to consider all possible independent variations of µ F , µ R , µ b by factors of two, eliminating all cases where the ratio of any two variation factors exceeds 2, and taking the total envelope. This turns out to be more aggressive and produces a smaller total uncertainty, because some of the variations (the µ b variations in particular) that in our approach are considered independent and added in quadrature, are simply contained within the overall envelope and thus have no effect on the final uncertainty.

Parametric uncertainties due to m b
We now discuss the settings we use for the b-quark mass. The b-quark mass is most precisely measured when defined in a renormalon-free short-distance scheme, such as the MS or 1S schemes. As our starting point and central input value we thus take the measured value of the MS mass m b (m b ) = 4.18 ± 0.03 GeV [31].
In principle, the best option would be to always use a renormalon-free mass renormalization scheme along with the corresponding measured value in all the places where the mass appears. These are the Yukawa coupling Y b , the m b dependence in the nonsingular parts of the coefficient functions, and the m b dependence of the PDF matching coefficients M ij (m b , µ b ) [see eq. (2.2)]. As mentioned in section 3.1 above, the Yukawa coupling is renormalized in the MS scheme and is directly obtained by evolving from µ = m b to µ R .
Unfortunately, most current PDF fits are performed with pole-scheme masses, including all PDF sets that currently enter the PDF4LHC15 combination, which is what we will utilize, see section 3.3 below. Furthermore, the fixed-order computations that we use as input and hence our nonsingular coefficient functions are currently obtained in the pole scheme. To be consistent, we thus also require m b in the pole scheme.
Since the pole mass has a leading renormalon ambiguity, the choice of its value is somewhat delicate. To reproduce as closely as possible the renormalon cancellation that would happen when properly translating the perturbative expressions from the pole scheme to the MS scheme, we have to convert the m b input value to a pole-mass value at the same loop order at which the perturbative series where m b appears (and that contains the cancelling renormalon) is used. In our case, this means we should translate it at 1-loop order because the proper NLO+NLL MS result would require an NNLO MS PDF set and would be affected by the 1-loop pole to MS conversion, To see this, note that the m b dependence of the fixed-order contributions first appears at LO and we work to NLO, so the scheme translation requires the 1-loop conversion. Equivalently, in the PDFs, the b-quark mass enters in the b-quark matching coefficient , which first appears in the NLO PDFs, while the NNLO PDFs contain its 1loop correction. Note that here it is again important that in our perturbative counting the fixed-order contributions and their singular limit contained in the resummed contributions are consistently included at the same order.
As a cross check, we have explicitly verified that evolving (with APFEL To evaluate the parametric uncertainty due to the uncertainty in the measured value of m b , we first note that the current world average for m b has an uncertainty of ±30 MeV [31]. Given the current tensions in different extractions of m b , this uncertainty might be considered too optimistic and one might want to consider the 2σ variation. For our purposes however the resulting uncertainty in Y b is small compared to other uncertainties, and so we will use the 1σ variation which is directly translated into the variation of the MS Yukawa coupling. Since the conversion to the pole scheme is performed at 1-loop, we also have to take into account the intrinsic uncertainty in the conversion, which is much larger than the uncertainty on m b itself. are always used in conjunction with the lower (upper) variation on m b in eq. (3.2). As we will see, the uncertainty due to m b will be small.
Finally, as mentioned earlier, when varying m b , the b-quark matching scale is kept fixed at its central value µ b = 4.58 GeV and µ F is also kept fixed at its central value µ F = (m H + 2m b )/4 with m b = 4.58 GeV.

Input PDFs and PDF uncertainties
As discussed in section 2, the resummation of collinear logarithms is entirely contained in the evolved 5FS PDFs, and these carry a physical dependence on m b and an unphysical dependence on the b-quark matching scale µ b . When computing the cross section, it is important that the input parameters used are consistent with the used PDF set. In particular, in our matched predictions we have to ensure that the value of m b in the computation of the fixed-order contributions is equal to the one present in the resummed contributions, i.e., in the PDF set, since otherwise the nonsingular terms would receive residual singular logarithmic terms arising from miscancellations.
To be able to use in a fully consistent way all the values we want for m b and µ b , for the central value predictions as well as the uncertainty estimation, we require dedicated PDFs that are not available by default. In principle, when changing the internal parameters of the PDFs, they should be refitted. In practice, the chosen value of m b when fitting the PDFs has a very small effect for all PDFs except for the b-quark PDF itself [34]. The reason is that the presently fitted data provides only a very weak direct constraint on the b-quark PDF or m b , which means that for all practical purposes the b-quark PDF is essentially being calculated from the fitted gluon and light-quark PDFs at the low scale µ 0 < µ b . Hence, we can also safely assume that refitting the PDFs for µ b = m b will not have much effect on the light-quark and gluon PDFs at µ 0 .
Given the above, we can take any PDF set with a given value of µ b = m b , and re-evolve it starting from a low scale µ 0 < µ b but using different values for m b , µ b , and the mass renormalization scheme. In other words, we use the light-quark and gluon PDFs at µ 0 as the input quantity and compute the b-quark PDF ourselves. (Note that this procedure is also typically used by PDF fitting groups when constructing fixed-flavor PDF sets from the fitted variable-flavor sets.) This approach is very useful because it opens the possibility of using any desired values for µ b and m b with any particular input PDF set.
Regarding the input PDFs at µ 0 , we use the combined PDF4LHC15 nnlo mc set [35][36][37][38][39]. All 101 PDF members are re-evolved from an initial scale µ 0 = 2 GeV for all the values of m b and µ b that we need. This is done using the latest version (≥ 2.8) of APFEL [33]. 4 We emphasize that re-evolving the PDF4LHC15 PDFs in this way is in fact more consistent   for the case of bbH than directly using the b-quark PDF from the combined set, because the prior sets (MMHT2014 [38], CT14 [39], NNPDF30 [37]) from which the PDF4LHC15 set is constructed use different values of m b . The re-evolved PDF sets for various different values of µ b and m b are publicly available at http://www.ge.infn.it/∼bonvini/bbh.
To compute the PDF uncertainties we follow the PDF4LHC prescription [35], using our re-evolved set with our central values for m b and µ b . That is, we compute the cross section for all 100 replicas, order the results in ascending order, and obtain a 68% confidence level interval by using the 17th result as the lower variation and the 84th result as the upper variation. We note that the distribution of cross sections we obtain is Gaussian to a very good approximation.

Results for the 13 TeV LHC
In this section, we present our numerical results for the inclusive bbH cross section for values of the Higgs boson mass in the range m H ∈ [50, 750] GeV. We also consider values of the Higgs mass different from m H = 125 GeV, as these are relevant for BSM scenarios in which the Higgs coupling to the bottom quark is enhanced. For simplicity we always use SM couplings -the cross section in many BSM scenarios can be obtained by rescaling Y b by an appropriate model-dependent factor [40][41][42].
We always use the simplified implementation without the strict expansion, which should still be valid at the lowest considered value of m H . The highest value is chosen only semirandomly [43,44]. The results for the cross section at NLO Yt]+NNLL partial are given in tables 1, 2, 3, and 4, respectively. The precise definitions of the different orders are dis- cussed in section 2. In the cross section tables we give the central value for the cross section together with the full breakdown of the perturbative uncertainties due to {µ F , µ R } and µ b as well as the parametric uncertainties due to m b and PDFs.
For each of the fixed-order (µ F and µ R ), resummation (µ b ), and parametric m b uncertainties we use the absolute value of the maximum deviation from the central result as the symmetric uncertainty. The PDF uncertainty is computed according to the PDF4LHC15 prescription as described in section 3.3 and is kept asymmetric. As discussed in section 3.1, the full perturbative uncertainty is obtained by adding the fixed-order and resummation uncertainties in quadrature. If one is only interested in the total bbH cross section, the  perturbative and parametric uncertainties can be added in quadrature. In a more complicated setup, e.g., a global fit, the total perturbative, m b , and PDF uncertainties should be treated as independent uncertainty sources (e.g. they correspond to independent nuisance parameters). This allows one to properly take into account their correlations with other predictions affected by the same physical uncertainty sources. For example, the parametric uncertainty due to m b or Y b should be correlated with the corresponding parametric uncertainty when bottom loop effects are included in ggH.  Our results are also illustrated in three figures. In figure 4, we show the LO+LL, 5 Yt]+NNLL partial cross sections as a function of the Higgs mass. Here, the uncertainty bands show the total perturbative uncertainty adding in quadrature the {µ F , µ R } and µ b uncertainties. The important message that follows from this figure is that we can see an excellent perturbative convergence between the orders. Additionally, the band for the NLO[ [pb] [pb]    300 GeV, and a positive one for larger masses. (The small fluctuations are due to the finite integration statistics in MC@NLO when including the interference terms.) Note that for large Higgs masses the uncertainties at NLO+NNLL partial get significantly smaller than at NLO+NLL. As discussed in section 2.3, this is expected since for large m H the 5FS perturbative counting is more appropriate and in this limit the NLO+NNLL partial formally improves the accuracy of the result, becoming as accurate as the full NNLO 5FS cross section. For smaller values around the physical Higgs mass of m H = 125 GeV, our counting is more appropriate, and only the full NNLO+NNLL should be regarded as a complete next-order result. Therefore, in this region the larger NLO+NLL uncertainty should be regarded as a safer estimate of the residual theory uncertainty. This is also nicely confirmed by the fact that in this range the NLO+NNLL partial uncertainty is essentially as large as at NLO+NLL.
In figures 5, 6, and 7 we give a visual breakdown of the results for m H = 125 GeV, m H = 475 GeV, and m H = 750 GeV. The matched results and uncertainty contributions are equivalent to the numbers provided in the tables. In addition, we also give a comparison with the corresponding pure fixed-order results at NLO with and without Y b Y t terms (green points). 6 By comparing the green NLO with the red NLO+NLL points we can see that the effects of the resummation are significant, resulting in a ∼ 30% increase for m H = 125 GeV, and even more at the higher Higgs masses, and moreoever this effect is not covered by the fixed-order scale variation band. This clearly shows that the resummation of b-quark collinear logarithms cannot be neglected for these mass values. Comparing the red points for NLO[Y 2 b +Y b Yt]+NLL and NLO[Y 2 b +Y b Yt]+NNLL partial we see the same features as we saw in fig. 4.
The four points on the right of figures 5, 6, and 7 show the breakdown of the uncertainties for our default result at NLO[Y 2 b +Y b Yt]+NLL. Clearly, the fixed-order uncertainty estimated by the µ F and µ R variations is the largest source of uncertainty in the predictions, both for moderate as well as large values of m H . The orange points show the resummation uncertainty estimated from the µ b variation, which, although smaller than the fixed-order uncertainty by roughly a factor of two, is nevertheless not negligible. We also recall that at the lower LO+LL order the resummation uncertainty plays a significant role, as was shown in ref. [16].
The parametric uncertainties due to m b (purple points) and PDFs (brown points) are subdominant. The m b uncertainties are small around 2% and decrease for higher Higgs masses to around 1%. The (asymmetric) PDF uncertainties, computed as described in section 3.3, are around 4% and smaller than both the {µ F , µ R } and µ b uncertainties.
One point to emphasize is that the uncertainty which is solely due to the parametric uncertainty in m b is much smaller than the perturbative µ b uncertainty, especially for larger Higgs masses. This shows the importance of distinguishing these two effects. If we were to identify µ b ≡ m b , as is commonly done, we would be left to choose between two undesirable options. Either we could vary their common value in the range eq. (3.3), which would essentially set the µ b uncertainty to zero. Or, we could vary their common value in a much larger range to account for the µ b uncertainty, which however would be unjustified for m b and blow up the parametric m b uncertainty. In contrast, by identifying and separating these two uncertainty sources, we are able to properly estimate each of them.

Conclusions
We have presented state-of-the-art predictions for the bbH cross section at 13 TeV obtained from a matched calculation [16] that consistently combines the fixed-order (4FS) contributions (which include the full b-quark mass dependence) with the all-order resummation of collinear logarithms of m b /m H . We provide results with and without including the effect of the interference of top-loop induced Higgs production process with the pure bottom-induced production proportional to Y b Y t . We also study the effect of two-loop contributions that formally contribute at NNLL order, finding that they are small and their effect fully captured by the uncertainty of our default NLO+NLL result.
We perform a detailed study of several sources of uncertainty in our results, both theoretical and parametric. The perturbative uncertainty from missing higher orders is estimated by varying the hard scales µ F and µ R , as well as the resummation scale µ b , which represents the threshold scale at which 4FS evolution is matched to 5FS evolution in the PDFs. We consider our resulting theory uncertainty as a reliable estimate, which is neither aggressive nor overly conservative. Furthermore, the parametric uncertainties due to the b-quark mass value and PDFs are evaluated. In particular, we discuss how to disentangle the unphysical dependence on the b-quark matching scale µ b from the purely parametric dependence on the b-quark mass m b , which requires the construction of dedicated 5F PDFs.
Our methodology to compute the matched prediction and to evaluate its uncertainties can be readily applied to other heavy-quark-initiated processes at the LHC. The code for our matched predictions will be available at http://www.ge.infn.it/∼bonvini/bbh. Our results represent the currently most complete predictions for the bbH cross section in the Standard Model and we are looking forward to a first measurement of this process during the coming LHC Run 2.