Higgs production in bottom quark fusion: Matching the 4- and 5-flavour schemes to third order in the strong coupling

We present analytic results for the partonic cross sections for the production of a Higgs boson via the fusion of two bottom quarks at N$^3$LO in QCD perturbation theory in the five-flavour scheme. We combine this perturbative result with NLO accurate predictions in the four-flavour scheme that include the full bottom quark mass dependence by appropriately removing any double-counting stemming from contributions included in both predictions. We thereby obtain state-of-the-art predictions for the inclusive production probability of a Higgs boson via bottom quark fusion at hadron colliders.


Introduction
Measuring precisely the properties of the Higgs boson, and possibly establishing the Standard Model (SM) of particle physics as the correct mechanism to explain the electroweak symmetry breaking, is one of the primary goals of the third run of the Large Hadron Collider (LHC) and its future upgrades. Since the SM Higgs boson couples to other particle species with a coupling strength proportional to their mass, measurements of the couplings of the Higgs boson to massive electroweak bosons and third generation fermions -the τ lepton as well as top and bottom quarks -are promising candidates to probe its interactions. The Yukawa coupling of the bottom quark is of particular interest, as several models of New Physics -like for example minimal supersymmetric extensions of the Standard Model -predict enhanced bottom Yukawa couplings (see chapter IV.2.2 of ref. [1]).
The interactions of the Higgs boson and the bottom quark can be probed at the LHC either through processes in which the Higgs decays to a pair of bottom quarks, or through processes in which it is produced from bottom quarks. In principle it is possible to directly constrain the bottom quark Yukawa coupling by measuring the decay of a Higgs boson into a bottom quark pair. However, even though this decay benefits from a large branching fraction, it is challenging to measure it precisely at a hadron collider due to the purely hadronic final state signature [2,3]. Moreover, any measurement of a Higgs boson decay necessarily relies on a precise prediction for its inclusive production crosssection. It is thus beneficial to study Higgs production processes at the LHC that involve bottom quarks. To measure the Yukawa coupling in this fashion, one particularly relevant production mechanism is that of the annihilation of two bottom quarks extracted from the colliding hadrons. The goal of this paper is thus to perform a phenomenological study of the production of a Higgs boson through bottom quark fusion.
Due to the small, but non-negligible, value of the bottom quark mass, there are two different ways in which one can model theoretical predictions for LHC processes involving bottom quarks. In the five-flavour scheme, the bottom quark is considered a massless parton. Consequently, all finite-mass effects are neglected, except for collinear logarithms that are resummed into the parton density functions. The five-flavour scheme has the advantage that the computation of higher-order corrections in the strong coupling constant is greatly simplified because all relevant quark species are massless (we neglect all top quark effects in the computations performed in the five-flavour scheme). In this scheme the inclusive bottom quark fusion cross section was computed through next-to-next-toleading order (NNLO) already almost two decades ago [4][5][6]. Very recently, a subset of the authors have computed for the first time the next-to-next-to-next-to-leading order (N 3 LO) corrections [7] (for a combination of the N 3 LO cross section with resummation of threshold logarithms, see ref. [8]). In a first part of this paper, we give more details on the structure of the partonic coefficient functions of ref. [7]. In particular, we make all the partonic coefficient functions publicly available as ancillary material attached to the arXiv submission of this paper. We also perform a detailed phenomenological analysis of Higgs production in bottom quark fusion, and we investigate the main sources of uncertainty that affect the cross section at N 3 LO.
While effects due to the non-zero mass of the bottom quark are expected to be small, they can nevertheless lead to sizeable effects, especially when compared to the level of precision with which the QCD effects are incorporated at N 3 LO. In the four-flavour scheme the bottom quark is treated as massive and is produced in the hard process, leading to higher final-state multiplicities. Consequently, Higgs production in bottom quark fusion is only known through next-to-leading order (NLO) in the four-flavour scheme [9][10][11]. As massive quarks cannot appear as initial state partons, all bottom quarks are generated from gluon splittings. While the non-zero mass protects the gluon splittings from collinear divergences, the four-flavour scheme is plagued by large logarithms involving the bottom quark mass which may spoil the convergence of the perturbative series. It is therefore desirable to combine the two schemes into a single prediction. Several methods to perform this combination have been proposed in the literature, ranging from purely phenomenological prescriptions [12] to theoretically well-grounded matching procedures [13][14][15][16]. So far, however, all these prescriptions have suffered from the fact that the equivalent of the NNLO result in the five-flavour scheme is only the leading order cross section in the four-flavour scheme. No matched prediction including all ingredients consistently through third order in the strong coupling has been obtained.
One of the main results of our paper is the first consistent matching of the four and five-flavour schemes through third order in the strong coupling. This is made possible by combining the N 3 LO result for the cross section of ref. [7] with the matching procedure of refs. [15,16]. In this way we are able to obtain the most precise predictions for this process, where all QCD and mass effects are included through third power in the strong coupling, and all logarithms of the bottom quark mass are resummed at leading power through next-to-next-to-leading logarithmic (NNLL) accuracy. Our paper is organised as follows: In section 2 we review inclusive Higgs production in the four and five-flavour schemes, and we introduce our notations and conventions. In section 3 we discuss the analytic structure of the partonic coefficient functions at N 3 LO in the five-flavour scheme, and in section 4 we present a detailed analysis of the different sources of uncertainty that affect the N 3 LO cross section. In section 5 we review the FONLL matching scheme, and in section 6 we present our results for the combination of the two schemes. In section 7 we draw our conclusions.

Higgs production in bottom quark fusion
In this section we review some basic facts about Higgs production in bottom quark fusion, and we introduce our notations and conventions. Using QCD factorisation, the cross section can be written as where µ F denotes the factorisation scale and the f i (x, µ 2 F ) denote the parton density functions (PDFs) to find a parton species i with momentum fraction x inside the proton. Thê σ ij denote the partonic cross sections to produce a Higgs boson from a collision of two partons i and j. Here we are interested in the production of a Higgs boson from the fusion of a pair of bottom quarks. More precisely, we focus on the part of the cross section proportional to y 2 b , where y b denotes the bottom quark Yukawa coupling. The sum runs over all active partons in the proton, i.e. gluons and all massless quark flavours.
Due to the small mass m b of the b quark compared to the mass m H of the Higgs boson, there are two ways in which eq. (2.1) can be interpreted. In the four-flavour scheme (4FS) the bottom quark is considered massive. Consequently, there is no PDF for the bottom quark and all finite mass effects are retained in the partonic cross sections. The non-zero mass also prevents the appearance of collinear singularities involving b quarks. Instead, the partonic cross sections develop collinear logarithms log Q 2 /m 2 b , where Q ∼ m H denotes the hard scale of the process. Given the hierarchy between the Higgs and the bottom quark masses, these logarithms may spoil the convergence of the perturbative series and need to be resummed to all orders in perturbation theory. This resummation is achieved by working in the five-flavour scheme (5FS), where the bottom quark is treated as massless and interpreted as a parton inside the proton. While the 5FS has the advantage that all collinear logarithms are resummed into the bottom quark PDF, it suffers from the fact that, unlike in the 4FS, the cross sections in the 5FS do not include any finite-m b non-logarithmic effects.
The 4FS and 5FS start to contribute at different orders in the perturbative expansion in the strong coupling constant α s . Indeed, in the 4FS (and under the assumption that there is no intrinsic bottom quark in the proton) the bottom quarks are generated perturbatively from gluon splittings, and therefore the perturbative expansion in the 4FS starts at order α 2 s . In the 5FS, instead, the bottom quark is considered a parton, and the leading-order cross section is proportional to α 0 s . Representative Feynman diagrams that contribute to each of the two schemes are shown in tab. 1.
The inclusive cross section in the 4FS can be written as ij τ, L f , L r , m b , α (4) s .

(2.2)
Here τ = s (µ 2 R ) denote the Yukawa coupling of the b-quark and the strong coupling constant for N f = 4 massless quark flavours. In the 4FS computation, the strong coupling is renormalised in the mixed scheme of ref. [17] in which the contribution from the four massless quark flavours is subtracted in the MS scheme, while the contribution from the massive bottom and top quarks running in the fermionic loop of the one-loop gluon selfenergy is subtracted on-shell. We define the normalisation factor Here, v is the vacuum expectation value of the Higgs field and n c refers to the number of colours. The renormalisation and factorisation scales are denoted by µ R and µ F respectively. Unless specified otherwise, all coupling constants are evaluated at a renormalisation scale µ R . The partonic luminosities are defined as the convolution of the corresponding four-flavour PDFs, L where the convolution is defined by (2.5) The sum in eq. (2.2) runs over all four massless quark flavours and the gluon. We find it convenient to use both integer numbers and explicit parton names as indices, e.g.: The partonic coefficient functions depend on the bottom quark pole mass and the logarithms of the factorisation and renormalisation scales: (2.7) They admit the perturbative expansion: s (µ 2 R )/π. The partonic coefficient functions in the 4FS are known (numerically) through NLO [9][10][11].
Similarly, the inclusive cross section for Higgs production in bottom quark fusion in the 5FS can be cast in the form Above, we again chose to normalise the partonic coefficient functions in the 5FS by the factorσ 0 defined in eq. (2.3). Throughout this paper we use the convention that X (n) denotes the quantity X computed in the n flavour scheme, and the notations introduced for the 4FS remains valid in the 5FS context. The main difference between the cross sections in the 4FS and 5FS in eqs. (2.2) and (2.9) is that in the 4FS the partonic coefficient functions have an explicit dependence on the bottom (pole) mass m b , and that the 4FS expression does not include the bottom quark into the sum over flavours. In particular, the coefficient functions η (5) ij admit the perturbative expansion The partonic coefficient functions in the 5FS are known at NLO [4,5] and NNLO [6]. Very recently also the N 3 LO corrections have become available [7]. We will review the results of ref. [7] in the next section.

Partonic coefficient functions in the 5FS
One of the main results of this paper are expressions for the N 3 LO corrections to the partonic coefficient functions η (5) ij for the production of a Higgs boson in bottom quark fusion. In this section we first discuss the general structure and computation of the partonic cross sections. We then explain the function space needed to represent the partonic coefficient functions. Finally, we give an alternative representation of our partonic coefficient functions in terms of expansions around different expansion points.

Structure of the partonic coefficient functions
At LO the only non-vanishing partonic coefficient functions have a bottom and anti-bottom quark in the initial state: The variable z is defined by where the x i are defined in eq. (2.1). Up to third order in the strong coupling constant there are eight distinct functions necessary in order to describe all partonic coefficient functions for different initial states. These eight functions are given by Above, g,b and b refer to a gluon, anti-bottom quark and bottom quark respectively, and q andq refer to a single quark and anti-quark that is not a bottom (anti-) quark. Results for the partonic coefficient functions at NLO and NNLO were computed in refs. [4][5][6].
The above functions were obtained by a subset of the authors at N 3 LO for the purposes of ref. [7]. Here, we present explicit results for these functions and make them publicly available in computer-readable form as ancillary material of this article. The computation of the the N 3 LO partonic coefficient functions follows the same strategy as that of the computation of the inclusive cross section for Higgs boson production through gluon fusion [18,19] and the inclusive Drell-Yan cross section [20]. In particular, the results were obtained by using the framework of reverse unitarity [21][22][23][24][25] in order to compute all required interferences of real and virtual amplitudes contributing to the N 3 LO cross section. The required phase-space and loop integrals were carried out implicitly by using integration-by-part (IBP) identities [26][27][28] together with the method of differential equations [29][30][31][32][33]. This method allows one to represent the required integrated and interfered amplitudes in terms of linear combinations of master integrals. Purely virtual amplitudes were first computed in ref. [34] using the master integrals from refs. [35][36][37][38][39][40][41], and recomputed and confirmed in ref. [7]. Contributions with one real parton in the final state were considered in refs. [42][43][44][45][46][47] and the master integrals we used for our calculation were documented in refs. [42,46]. Master integrals with two and three real partons were obtained for the purpose of ref. [19] and are based on results from refs. [18,[48][49][50][51][52].
We work in the MS-scheme in conventional dimensional regularisation. The ultraviolet (UV) counterterm for the strong coupling constant has been determined through five loops in refs. [53][54][55][56][57]. The renormalisation constant for the Yukawa coupling is identical to the quark mass renormalisation constant of QCD in the MS-scheme [6,55,[58][59][60]. Infrared (IR) divergences are absorbed into the definition of the PDFs using mass factorisation at N 3 LO [61][62][63]. The mass factorisation involves convoluting lower-order partonic cross sections with the three-loop splitting functions of refs. [64][65][66]. We have computed all the convolutions analytically in z space using the PolyLogTools package [67]. After combining our interfered matrix elements with the UV and PDF-IR counterterms we send the dimensional regulator to zero and obtain our final results.
The partonic coefficient functions for a bottom and anti-bottom quark in the initial state contain distributions in the variable z that were already obtained in ref. [68]. We checked that our computation agrees with this result. We refer to these contributions as soft-virtual (SV) contributions and to the non-distribution-valued part of the partonic coefficient functions as the regular part. Consequently, we split our partonic coefficient functions into regular and SV parts.
The coefficients of the leading two powers of logarithms log 5 (1 − z) and log 4 (1 − z) of the regular part can be derived using the method of physical evolution kernels of refs. [69][70][71] and agree with our results. Furthermore, we investigated the structure of the partonic cross section in the high energy limit. The leading logarithmic behaviour of the partonic coefficient function could be computed along the lines of ref. [72] for the Drell-Yan cross section. To the best of our knowledge, for the bbH cross section this computation currently does not exist. However, the structure we observe agrees with our expectation as we observe only a single logarithm at N 3 LO and the coefficient of this logarithm appears to be universal. Explicitly, we find at NNLO η (5,2) bb , η (5,2) bg , η (5,2) bq , η (5,2) bq , η (5,2) bb , η (5,2) gg , η

Analytic results for the partonic coefficient functions
Our partonic coefficient functions can be expressed in terms of the same set of functions used to represent the results of ref. [19]. For convenience, we repeat here the most essential definitions. We define an iterated integral as Our partonic coefficient functions can be expressed in terms of linear combinations of the above iterated integrals with algebraic functions in z as prefactors. The required integration kernels ω i (z) are drawn from the set The functions t ij are the solutions to the differential equation These functions can be represented in terms of elliptic integrals. If an iterated integral only contains integration kernels corresponding to the first three elements of eq. (3.8) then it belongs to the class of well known harmonic polylogarithms [73] (HPLs). More generally, if no integration kernel involving the functions t ij (z) appears, then the iterated integral can be expressed in terms of multiple polylogarithms (MPLs) [74] evaluated at algebraic arguments. If also integration kernels involving some t ij (z) appear, the iterated integral cannot be expressed in terms of MPLs alone, but it belongs to a more general class of functions related to elliptic curves. Currently it is unknown if these iterated integrals can be expressed in terms of elliptic multiple polylogarithms [75] or iterated integrals of modular forms [76,77], which have recently appeared in the context of multiloop calculations. For the purposes we choose to represent our partonic coefficient functions in terms of HPLs and iterated integrals as in eq. (3.7). In order to evaluate the partonic coefficient functions numerically, we find it useful to express them in terms of generalised power series expansions. In ref. [19] it was discussed how such iterated integrals relate to one another and how they can be expanded around different numerical points. The physical domain for our partonic coefficient functions is given by z ∈ [0, 1]. By studying the singularities of the functions expressing the partonic coefficients, we can deduce that a generalised power series expansion of the coefficient functions around the point z = 1 is convergent within the entire physical domain z ∈ [0, 1]. However, in order to reduce the number of terms required to evaluate the partonic coefficient functions to a given numerical accuracy, we choose to expand them around two additional points.
1. z ∈ [ 3 4 , 1]: In this interval we expand around the point z = 1 and define the variablē z = 1 − z for convenience. The power series inz is convergent within the entire unit interval but further sub-divisions are desirable in order to avoid loss of numerical accuracy when including only few orders in the expansion. We provide 50 terms in the series expansion aroundz = 0.
2. z ∈ [ 1 13 , 3 4 ]: Within this interval we expand around the point z = 1 2 and define the variable w = 1 2 − z for convenience. We provide 200 terms in the expansion around w = 0. Formally, this expansion around w = 0 is convergent in the entire interval z ∈ [0, 1].
3. z ∈ [0, 1 13 ]: In this interval we expand our partonic coefficient functions around the point z = 0 and we provide 100 terms in this expansion. Contrary to the previous two expansions, this one is only convergent within the interval z ∈ [0, − 1 2 11 − 5 With the provided number of terms in the different series expansions the partonic coefficient functions can be evaluated with a relative numerical precision of at least 10 −10 . While the formal radius of convergence of the different expansions listed above refers to the validity of the expansions, we advise to stick to the suggested intervals in order to achieve a numerical accuracy of the partonic coefficient function of at least ten significant digits. We provide digital files containing the partonic coefficient functions through N 3 LO as ancillary material of this article. Figure 1 shows the individual regular partonic coefficient functions for the eight different partonic initial states. We use the PDF4LHC15 nnlo mc set [78] parton distribution functions if not stated otherwise explicitly. Throughout this article we only consider contributions proportional to O(y 2 b ). We however remind the reader that bottom quark fusion contributions proportional to O(y b y t ) and O(y 2 t ) are relevant as already discussed in refs. [11,79].

Dependence on the perturbative scales
Through N 3 LO our cross section is independent of the factorisation and renormalisation scales. However, the numerical values for cross section predictions will vary depending on the choice of the values for the perturbative scales since the evolution of the PDFs, the strong coupling and the Yukawa coupling are performed in a resummed fashion. At NLO it was argued in refs. [80][81][82][83] that the t-channel singularity in the gluon-initiated process gb → bH leads to a collinear logarithm of the form log(4µ F /m H ) in the inclusive cross section and that consequently a low value for the factorisation scale should be preferred. In refs. [1,6,84] it was observed that choosing low factorisation scales leads to faster stabilisation of the perturbative series. We consequently follow this approach and choose the as the central values for our perturbative scales: (4.1) Figure 2 shows the dependence of the hadronic cross section on the factorisation (left) and renormalisation (right) scales. The bands in the two figures are obtained by varying one particular scale up and down by a factor of two around the central value. We observe in fig. 2 that including higher-order perturbative corrections reduces the dependence of the hadronic cross section on both perturbative scales since the span of the bands is reduced by the inclusion of higher-order corrections. We also notice that the perturbative series is relatively well behaved for low values of the factorisation scale. This strengthens the case for our choice of central value for the factorisation scale. Figure 3 shows the cross section for the production of a Higgs boson in bottom quark fusion for various hadron collider energies. Different colours refer to different orders of the perturbative expansion, and the bands correspond to varying the perturbative scales by a factor of two around their central value while satisfying the inequality (7-point variation)   = m H . We find that the nominal value of the cross section at N 3 LO is comparable for these two choices. However, the perturbative corrections are much larger in the latter case, thus further supporting our choice of a low factorisation scale for this process.

PDF and α s uncertainties
We take the PDFs and the strong coupling constant as external input. These quantities are naturally associated with an uncertainty that we asses following the guidelines of the providers of these quantities. In particular, we use the PDF4LHC15 nnlo mc set [78] as our default PDF set and follow the Monto-Carlo prescription outlined in ref. [78] in order to determine the PDF uncertainty of our cross section. In particular, following this prescription the hadronic cross section is computed with 100 different PDF sets and the resulting values are then ordered by nominal size. The PDF uncertainty is then determined by δ(PDF) = ± σ where, the σ (5) i corresponds to the i th member of the ordered set. As a central value for cross section predictions is recommended to bē 16 .
(4.4) Figure 4 shows the resulting PDF uncertainty as a function of the collider energy. Furthermore, we compare different PDF sets with prediction based on the PDF4LHC15 set. In particular we study the sets • NNPDF30 nnlo as 0118 [88] , • NNPDF31 nnlo as 0118 [89] .
We observe a sizable PDF uncertainty from 7−9%. Comparing the predictions based on the PDF4LHC15 set with the other PDF set we see significant differences. The PDF4LHC15 set itself is a statistical combination of the CT14, MMHT and NNPDF3.0 sets, and we observe in fig. 4 that indeed the resulting prediction is in between the three input sets. NNPDF3.1 is an updated version of NNPDF3.0 and technically supersedes the latter. Consequently, it is possible that a combination of CT14, MMHT and NNPDF3.1 into an updated version of a PDF4LHC combination would lead to a significantly lower central prediction of the bbH cross section. However, such a study is beyond the scope of this article.
In order asses the uncertainty due to the imprecise knowledge of the strong coupling constant, the authors of ref. [78] provide two PDF sets within the PDF4LHC15 nnlo mc pdfas set that allow to vary the strong coupling constant by ±0.0015 in a correlated fashion. The associated uncertainty is computed as (4.5) Following the recommendation of ref. [78] this uncertainty can then be combined in quadrature with the PDF uncertainty: The definition of the value for the prediction of the inclusive cross section in eq. (4.4) can be compared with the prediction that is obtained with the central member of the PDF4LHC15 nnlo mc set. Their ratio is shown in fig. 4 in green on the right. While there is a non-negligible difference the two predictions are compatible within the PDF uncertainties.

PDF theory uncertainty
PDFs are currently determined using NNLO cross sections as input for their extraction from a wide set of measurements. Consequently, we refer to these PDFs as NNLO PDFs. Since our cross section is computed at N 3 LO this leads to a mismatch that can ultimately be remedied by using N 3 LO cross sections for the PDF extraction. In the meantime we estimate the potential impact of this mismatch on our cross section predictions. In ref. [90] a prescription was introduced that studies the variation of the NNLO cross section as NNLO or NLO PDFs are used. This defines the PDF theory uncertainty Here, the factor 1 2 is introduced as it is expected that this effect becomes smaller at N 3 LO compared to NNLO. Figure 5a displays δ(PDF-TH) as a function of the collider energy. Throughout this uncertainty is smaller than the PDF uncertainty. We interpret the numerical crossing  point at about 60 TeV as a coincidence and a simple consequence of the method we use to estimate this uncertainty. Consequently, this does not mean that there is no PDF theory uncertainty for a 60 TeV collider and we assign always at least a 1% uncertainty whenever the prescription of eq. (4.7) falls below.

Bottom quark mass uncertainty
According to the PDG [91] the bottom quark mass in the MS-scheme is determined to be Since the cross section in the 5FS is proportional to the square of the bottom quark mass the hadronic bbH cross section is affected by the corresponding uncertainty The bottom quark mass evaluated at the renormalisation scale is completely factorised from the partonic coefficient functions as can be seen in eq. (2.9). We perform the scale evolution via a numerical solution to the evolution equation using anomalous dimensions at (n + 1) perturbative order in order to compute the N n LO cross section: (4.10) The constants γ (i) are taken from ref. [92]. Overall, we find that truncating the anomalous dimension at the (n+1) th order slightly improves the rate of convergence of the perturbative expansion. However, we find that the value of m b (µ cent.

R
= m H ) changes at the sub-permille level if we are using three-loop or four-loop anomalous dimensions, cf. tab. 2. Consequently, we do not assign an additional uncertainty for the exact implementation of the bottom quark mass.  Alternatively to the MS-scheme, we derive predictions for the bbH cross section using the on-shell bottom quark mass. Using the three-loop conversion relation of refs. [93,94] we find that the on-shell bottom quark mass is given by (4.11) Figure 5b shows the ratio of the bbH cross section with computed with on-shell bottom quark mass at different perturbative orders to the same computed with MS mass at N 3 LO. We observe that as the perturbative order is increased the predictions based on different mass schemes approach each other. However, the perturbative convergence of the cross section predictions using on the on-shell mass is quite slow. In part this can be attributed to the fact that we are not resumming the mass evolution as in the MS-scheme. At LO the bottom quark mass in eq. (2.3) is now evaluated with its on-shell value and the ratio of the normalisation factorsσ 0 of the two different schemes is ∼ 2.67. Furthermore, it is well known that the conversion from MS to on-shell scheme is affected by large perturbative corrections (see for example refs. [93,94]). Based on the above observations we recommend the treatment of the bottom quark mass as in our default set-up.

The FONLL matching procedure
In order to have precise theoretical predictions it is desirable to combine the 4FS and 5FS into a single prediction which retains finite mass effects through a certain order in perturbation theory while at the same time resumming the collinear logarithms to all orders in the strong coupling. Various methods have been proposed in the literature to combine the two schemes [12][13][14][15][16]. Here we focus on the so-called FONLL scheme, first introduced in refs. [95,96] for hadron production in hadronic collisions and deep inelastic scattering and recently applied to Higgs [15,16] and Z-boson [97] production in bottom quark fusion in proton collisions. The original versions of refs. [15,16], however, contained some misprints, and we therefore reproduce all formulas here for completeness.
At all perturbative orders, the cross sections in the 4FS and 5FS in eqs. (2.2) and (2.9) are identical up to power suppressed terms (and possibly up to non-perturbative effects encoded in the different PDFs), A similar relation, however, does not hold at the level of the partonic coefficient functions calculated in the two schemes. Indeed, the coefficient functions in the 4FS develop logarithmic divergencies in the limit of a vanishing bottom quark mass, which are not captured by the coefficient functions in the 5FS. Instead, these m b -dependent logarithms are encoded (and resummed) into the PDFs and the strong coupling constant in the 5FS. The starting point of the FONLL method is to express both computations in terms of a common set of PDFs and α s , namely the ones in the 5FS. The relation between the strong coupling constant and the PDFs in the two schemes takes the form, . The explicit form of the kernels K ij relevant here can be obtained from ref. [98]. In particular, they have the property that K ij = δ ij δ(1 − x) + O(α s ) for |i| = 5 and K ij = O(α s ) for i = ±5. This allows us to invert eq. (5.2) order by order in the coupling, and to express the cross section in the 4FS in eq. (2.2) in terms of the coupling and the PDFs in the 5FS, where the partonic coefficient functions admit the perturbative expansion: Through third order in the strong coupling, the relation between the partonic coefficient functions in eqs. (2.2) and (5.3) reads, with T f = 1 2 . By inserting the expression of the PDFs in the 4FS in terms of those in the 5FS back into eq. (5.2), we can re-express the b-PDF entirely in terms of the PDFs for the other parton flavours in the 5FS. Through the order we need it, this relation reads Note that the bottom and anti-bottom distributions are only identical through the first two orders in the strong coupling constant, and they will start to differ starting from O(a 3 s ) (cf., e.g., ref. [99]). The kernels A (k) bg and A (2) bΣ can be found in ref. [98]. Inserting this relation into eq. (2.9), we can write the cross section in the 5FS as σ (4−5) in a way that does not involve the b-PDF and which is formally equal to σ (5) up to third order in α (5) s , The partonic coefficient functions A ij can be expressed in terms of the partonic coefficient functions in the 5FS in eq. (2.9) and the kernels in eq. (5.6). In the following we only show this relation for µ R = µ F , and we suppress the dependence of all functions on their arguments for readability. If we denote the coefficient of = 0 for all other channels. We have performed all these convolutions analytically using the PolyLogTools package [67]. The analytic expressions for the convolutions in terms of multiple polylogarithms are provided as ancillary material with the arXiv submission.
Using these definitions, we can write the cross section in the FONLL scheme as The fact that σ (4−5) removes the overlap between the cross sections computed in the 4FS and 5FS is guaranteed by noting that Using a straightforward rearrangement of terms, we can cast eq. (5.9) into the alternative form, σ matched = σ (4) +σ (5) −σ (4−5) , andσ (5) collects only those channels in the 5FS that have a b-quark in the initial state (we suppress again the dependence on all arguments for readability) With the completion of the N 3 LO corrections in 5FS, we have now for the first time the possibility to compute all ingredients in eq. (5.9) consistently through third order in the strong coupling. The phenomenological impact of these corrections will be explored in the remainder of this paper.

Phenomenological results
In this section we present our results for the inclusive cross section matched according to the FONLL procedure through third order in the strong coupling. We work with a Higgs mass of m H = 125 GeV and the pole mass of the bottom quark is m b = 4.58 GeV. The strong coupling and the Yukawa coupling are evaluated at the renormalisation scale µ 2 R using three-loop running in the MS-scheme [6,[53][54][55][56][57][58][59][60], and we start the evolution from α s (m 2 Z ) = 0.118 and m b (m b ) = 4.18 GeV. We choose to work with the PDF set of ref. [14,100], which is based on the combined PDF4LHC15 nnlo mc set [78], but starting from a low scale where there is no bottom quark, and then performing the evolution to higher scales using a consistent value of the bottom pole mass throughout.
The 4FS results are generated using MadGraph5 aMC@NLO [101]. The computation of the one-loop amplitudes is carried out with the module MadLoop [101,102], which generates the loop integrand using an in-house implementation of the OpenLoops optimisation [103]. The loop integrals are then evaluated by switching dynamically between two one-loop reduction techniques: OPP [104] or Laurent-series expansions [105] that are performed at the integrand level, and methods applied at the tensor integral level [106][107][108]. These reduction techniques have been automated in tools that MadLoop interfaces to: CutTools [109], Ninja [110,111] and COLLIER [112]. The renormalisation of the bottom quark Yukawa coupling is performed by default in the on-shell scheme in Mad-Graph5 aMC@NLO [101]. In order to renormalise this quantity in the MS-scheme instead (and correctly account for the running of y b (µ R ) in this case), we must perform adjustments 1 of the process output identical to those considered in ref. [11]. Finally, we note that the top mass contributions of order O(y 2 b ) (i.e. but not the ones involving y t ) are included in the NLO 4FS computation (whereas they are not in the N 3 LO 5FS computation). These top-quark contributions come in through corrections of the triple-gluon vertex as well as the gluon propagator (and therefore its wavefunction counterterm). We stress that considering the top-quark contribution only in the 4FS part of the computation does not spoil the consistency of the matching procedure presented in section 5. In addition, we have verified that its numerical impact is at the permille level only.
Before we present our results, let us briefly comment on different ways to implement the FONLL matching procedure. More specifically, in refs. [15,16] three different scenarios were considered: • FONLL-A: All ingredients in eq. (5.9) are included through O(α 2 s ). This corresponds to matching the 5FS at NNLO to the 4FS at LO, and all collinear logarithms are resummed at next-to-next-to-leading logarithmic accuracy (NNLL   [14,100] with the bottom mass set to infinity. In contrast, the 4FS cross-sections entering the FONLL matching procedure in eq. (5.5) were computed using the same PDF set as for the 5FS computation (i.e. PDF set evolved using a bottom mass set to 4.58 GeV).
results for Higgs production in bottom-quark fusion using the FONLL-A prescription have first been obtained in ref. [15].
• FONLL-C: All ingredients in eq. (5.9) are included through O(α 3 s ). This corresponds to matching the 5FS at N 3 LO to the 4FS at NLO, so that all collinear logarithms are resummed at NNLL. Phenomenological results using the FONLL-C are presented for the first time in this paper.
In figs. 6a and 6b we show the variation of the 4FS, 5FS, and matched results with the renormalisation or factorisation scale, with the other scale held fixed. We observe that FONLL-C prediction increases the value of the N 3 LO 5FS result by roughly 2% over the whole range of scales considered, while maintaining the very reduced sensitivity to the residual scale dependence of the N 3 LO result. This is at variance with the matching at the previous order (FONLL-A), where the matched prediction only resulted in a tiny increase of the 5FS cross section at NNLO [15]. Finally, we observe that the FONLL-B prescription leads to a substantial increase of the cross section compared to the 5FS NNLO result. The FONLL-B prescription misses the contributions the b-initiated channels at N 3 LO, which give large and negative contributions to the cross section. More precisely, the FONLL-B prescription does no satisfy eq. (5.10) since it considers all 4FS contribution B (3) ij while ignoring (i.e. effectively setting to zero) the 5FS counterpart pieces η (5,3) bi , η (5,3) bi , η ij . As a consequence, it seems that for this particular process the FONLL-B prescription does not give a reliable estimate of the value of the cross section at O(α 3 s ). This underlines the need to include the N 3 LO 5FS prediction.  Table 3: FONLL-C (N 3 LO 5FS matched to NLO 4FS) predictions for the bbH cross section at different collider energies and associated uncertainties.
In tab. 3 we present results for the matched cross section for various representative collider energies. We estimate the uncertainty due to the truncation of the perturbative series by varying the factorisation and renormalisation scales independently up and down by a factor around the central values (µ F , µ R ) = ((m H +2m b )/4, m H ) within the constraint of eq. (4.2). This choice for the central scales was discussed in section 4. Furthermore, we quote the PDF and strong coupling uncertainty δ(α S + PDF), the PDF theory uncertainty δ(PDF-TH) and the bottom quark mass uncertainty δ(m b ) that we asses based on the five-flavour cross section as outlined in section 4.

Conclusion
In this paper we have performed a detailed phenomenological study of Higgs production in bottom quark fusion. In a first part of the paper we have focused on the N 3 LO cross section in the 5FS. We described the structure of the analytic partonic coefficient function for this cross section as well as for the matching contributionσ (4−5) and include it in electronically readable form together with the arXiv submission of this article. Next, we elaborated on the phenomenological analysis of ref. [7]. We have studied the dependence of the cross section of the renormalisation and factorisation scales. We observe a convergent behaviour of the perturbative series, provided that the factorisation scale is set to a relatively low value. This corroborates similar conclusions drawn based on the behaviour of the cross section at lower orders, and gives further support for this unconventionally low choice of the factorisation scale. We have also studied other sources of uncertainty that may affect our prediction for the cross section, including the effects due to PDFs and the strong coupling constant, as well as the value of the bottom quark mass that is used in the computation.
In a second part of the paper we have combined our N 3 LO computation in the 5FS with the NLO cross section in the 4FS computed with MadGraph5 aMC@NLO. The overlap between the two schemes is removed using the FONLL matching procedure, first applied to Higgs production in bottom quark fusion in refs. [15,16]. The novelty of our computation lies in the fact that for the first time we can compute all quantities that enter the combination consistently through third order in the strong coupling. We find that the effect of the matching is non-negligible, increasing the value of the 5FS N 3 LO cross section by roughly 2%. We note that this increase is of the same order as the scale dependence at N 3 LO. We also find that previous attempts to match the two schemes through third order in the strong coupling without including the complete N 3 LO calculation had led to a substantially different answer. The reason is that the b-initiated channels at N 3 LO give a large and negative contribution to the cross section, an effect which was not captured by previous calculations.
To conclude, we have presented the most precise prediction for the inclusive bottom quark fusion cross section by combining the most precise calculations in both the 4FS and 5FS. The non-negligible effect of the N 3 LO corrections underlines once more the need for calculations at this order for the precision physics program at the LHC, and we expect that our results will play a role in the study of the interactions of the bottom quark and the Higgs bosons, both at the LHC and at future hadron colliders.