Bottom-induced contributions to Higgs plus jet at next-to-next-to-leading order

We present a next-to-next-to-leading order (NNLO) QCD calculation of the bottom-induced contributions to the production of a Higgs boson plus a jet, i.e. the process pp → H + j to Oyb2αs3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{O}\left({y}_b^2{\alpha}_s^3\right) $$\end{document}. We work in the five-flavor scheme (5FS) in which the bottom quark mass is retained only in the coupling to the Higgs boson. Our calculation uses N-jettiness slicing to regulate infrared divergences, allowing for fully-differential predictions for collider observables. After extensively validating the methodology, we present results for the 13 TeV LHC. Our NNLO predictions show a marked improvement in the overall renormalization and factorization scale dependence, the latter of which proves to be particularly troublesome for 5FS calculations at lower orders. In addition, using the same methodology we present a NNLO computation of bb¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ b\overline{b} $$\end{document} → H. Our results are implemented into MCFM.


Introduction
Since its discovery nearly a decade ago [1,2], the Higgs boson has become an established part of the particle physics landscape, and a significant amount of research effort has been devoted to the greater understanding of its properties. After the initial establishment of its intrinsic properties such as its CP and spin [3,4], the focus has shifted to obtaining precision measurements of the couplings of the Higgs boson to the other particles of the Standard Model (SM) as well as to itself. The Higgs boson self-coupling is a particularly pressing measurement to obtain, since it will allow for a more detailed study of the electroweak symmetry breaking potential in the Standard Model. Although the Higgs self-coupling is fully predicted from known parameters in the Standard Model, it is often sensitive to extensions of the SM (referred to as BSM) that change the nature of the electroweak potential. Similarly, it is also of vital importance to constrain the couplings of the Higgs boson to other particles in the Standard Model, namely the W and Z bosons and massive fermions. The LHC has made significant progress over the last decade [5,6], and will continue to improve upon existing results over the forthcoming Run III. Further in the future, measurements with sub-percentage uncertainties will require a collider with upgraded capabilities, for which serious planning is now underway [7,8]. In all of these endeavors, precision predictions for differential distributions in the Standard Model are critical in order to avoid the situation in which theoretical uncertainties become the dominant source of error.

JHEP05(2021)045
Of the Higgs interactions, one of the most fascinating to study is the coupling between the Higgs boson and third-generation fermions, the top and bottom quarks and tau leptons. Among these, the bottom quark is unique, as there are two ways to gleam insights into its coupling: through Higgs boson production or its direct decay to the quarks themselves. Given the large hierarchy in mass, the top quark dominates Higgs boson production at the LHC through the gluon-fusion mechanism and therefore probing the bottom Yukawa coupling through Higgs production is challenging. On the other hand, the SM Higgs boson copiously decays to bottom quarks with a branching fraction of around 50%, and therefore the bottom Yukawa coupling dominates the Higgs decay width, propagating to all other (on-shell) measurements of the Higgs boson at the LHC. Determining the bottom-Higgs coupling as precisely as possible is thus an essential requirement of the future experimental high-energy physics program. For example, in extensions of the SM, extended Higgs sectors typically modify the coupling of the 125-GeV Higgs boson to up-and down-type fermions, and can lead to enhanced production cross sections [9][10][11][12].
Given the immense interest, there have been many theoretical studies of processes involving the Higgs boson and the bottom quark at hadron colliders. When making predictions at the LHC, one must first decide how to handle the mass of the bottom quark. Since the bottom quark is heavier than the proton, a natural choice is to keep the mass of the bottom quark in the calculation and exclude it from initial-state contributions. This scheme is known as the four-flavor scheme (4FS), due to the number of active initial-state flavors in the proton. The leading-order production mechanism at the LHC in the 4FS is thus the process gg → bbH (plus a sub-dominant qq-initiated contribution). The 4FS has the advantage that no approximations are made in regards to the kinematics, which is particularly helpful in relation to final-state bottom-quark tagging, since single b-tagged jets can be isolated without theoretical issues, i.e. no jet cuts are required even though the final state contains two partons. However, a major drawback of the 4FS scheme arises from the occurrence of large logarithms of the form α s log(m 2 b /m 2 H ). When computing cross sections as perturbation series directly in α s , these logarithms induce large corrections and spoil the convergence of the series. One way to ameliorate this problem is to resum the logarithms when possible. Initial-state collinear logarithms can be resummed into the parton distribution functions (PDFs). This introduces the five-flavor scheme (5FS) in which the bottom quark contributes to the initial-state PDFs, and the mass of the quark is neglected in the remaining kinematics. At leading order, Higgs boson production in the 5FS thus proceeds directly through bottom-quark fusion bb → H. It is worth noting that, since final-state collinear splittings g → bb are not resummed, care must be taken in the 5FS when applying b-tagging requirements and comparisons are made to experimental data (i.e. one should try to remove jets arising from gluon splitting from the experimental analyses).

JHEP05(2021)045
A further intricacy relating the 5FS to the 4FS comes from the dependence on the factorization scale through the bottom-quark PDF in the 5FS. It was noted in the earliest calculations of bottom-quark fusion at NLO in the 5FS [13,14] that the higher-order corrections were large and resulted in cross sections that could have differences of an order of magnitude from the 4FS LO result. Detailed analysis in ref. [16] (following the arguments of ref. [35]) illustrated that the choice of a central factorization scale of m H was too high for the process, and that a scale of around m H /4 was more appropriate to ensure the reliability of collinear factorization. Predictions made with a central factorization scale choice in this region showed much better perturbative convergence, and a broader compatibility with the 4FS result. Nevertheless, a strong dependence on the unphysical factorization scale should be seen as a negative feature of the 5FS when used at LO and NLO. Subsequent NNLO and N3LO predictions for bottom-quark fusion [21,22,33] show a significant reduction of this problem and motivate our computation of H + j in the 5FS at NNLO.
On the Higgs boson decay side the theoretical situation is also under good control. Again, there are two scheme possibilities which can be considered as the decay versions of those discussed above. A commonly-used approximation is to retain the mass of the bottom quark only in the coupling to the Higgs boson (the decay equivalent of the 5FS), and there have been many detailed theoretical studies of this process (see e.g. [36][37][38]), which is now known up to O(α 4 s ) inclusively [39] and at N3LO differentially [40]. Increasing the complexity of the calculation, one can include the mass of the bottom quarks fully (the 4FS equivalent), and in this setup recent calculations have pushed the accuracy to NNLO for fully-differential predictions [41,42].
The continuing maturation of the experimental analyses at the LHC has had a twofold impact on Higgs boson studies. Firstly, the increased statistics and precision have allowed for an extensive range of Higgs boson observables to be studied, including Higgs-plusmultiple-jet production (see e.g. [43,44]). Secondly, comparison of data to theory has highlighted the need for increased precision on the theoretical front, emphasizing the importance of NNLO predictions in QCD. Over the last few years significant progress has been made, resulting in several independent calculations of Higgs-plus-jet at NNLO in the effective field theory (EFT) in which the top quark is integrated out [45][46][47][48][49][50][51]. Impressively, in ref. [52] the accompanying jet was integrated out of the calculation, allowing for a computation of the pp → H + X differential cross section at N3LO accuracy.
Our aim in this paper is to provide a similar level of theoretical accuracy for the bottom quark-initiated contribution to H + j as it is present for the dominant EFT production mechanism. In order to do so, we will work in the 5FS and treat the bottom quark as massless everywhere except in the coupling between the Higgs boson and the bottom quarks. The different calculations of H + j at NNLO available in the literature used a variety of infrared-regulating techniques. For brevity we do not describe them all in detail here, but focus on the pieces pertinent to this paper. One of the initial calculations of H + j [47] used a non-local slicing procedure (based on the event shape N -jettiness). The results obtained in this paper seemed to be in conflict with those obtained in other studies based upon local subtraction techniques. A detailed study in ref. [51] showed that the two methodologies do indeed produce the same results, but that power corrections arising from JHEP05(2021)045 the approximate form of the factorization formula used in the slicing procedure must be handled carefully. In this paper we will use the same methodology and follow the same techniques as shown in ref. [51] to control and estimate the remaining power corrections. In ref. [51] the gg channel was shown to have the worst power corrections for EFT H + j production. Thankfully for our calculation, this channel does not appear at leading order and we therefore expect power corrections to be easier to control.
Our paper proceeds as follows. In section 2 we present a brief overview of the technical details of our calculation, while section 3 discusses its validation. We present results for the 13 TeV LHC in section 4, and finally we draw our conclusions in section 5.

General overview
The primary focus of this paper is the calculation of the NNLO QCD corrections to the bottom-induced contributions to Higgs plus one jet at the LHC. Given its phenomenological relevance and role as a check of our calculation, we will also present results for the bottom quark fusion process at NNLO (i.e. bottom-induced contributions to Higgs plus zero jets).
As discussed in the introduction, the most important theoretical choice when considering bottom-quark processes at hadron colliders is how to treat m b , i.e. whether to work in the 5-flavor (5FS) or 4-flavor (4FS) scheme. In this paper we work in the five-flavor scheme, which will allow us to extend the computation to NNLO accuracy. Representative Feynman diagrams relevant for our calculation at this order are shown in figure 1. We recall that in the 5FS the bottom quark mass is taken to zero and bottom quarks have a non-zero contribution to the PDFs. At first glance, the 5FS scheme may appear not to be useful for computing H + b related processes, since by setting the bottom quark mass to zero the bottom Yukawa coupling should also be taken to zero. In order to circumvent this problem, we work in the mixed-renormalization scheme, in which the bottom quark Yukawa is taken in the MS scheme, and the bottom quark mass, used in propagators and in the relativistic kinematics, is taken in the on-shell scheme. This scheme allows one to JHEP05(2021)045 take the limit m OS b → 0 while keeping the Yukawa coupling non-zero. This scheme has two advantages in QCD calculations. Firstly, it allows for a robust definition of the 5FS for H + b amplitudes. Secondly, by evolving the scale in the running Yukawa coupling to µ R (i.e. m H ), one avoids large logarithms which arise in the OS scheme at higher orders, and as a result the perturbative corrections are under better control. Downsides of the mixed scheme include breaking the relationship between the OS mass and the MS mass [53,54] and an inability to consistently renormalize higher-order corrections in the electroweak coupling [31]. Nevertheless, the reduction in sensitivity to collinear initial-state logarithms (at the cost of a strong dependence on the factorization scale at LO), and the ability to pursue higher-order corrections, renders the 5FS along with the mixed renormalization scheme a very useful theoretical construct for LHC computations.

Technical details
For the bottom-induced H + j process at NNLO, three phase-space topologies contribute (see figure 1), corresponding to the double-virtual, real-virtual, and double-real corrections to the underlying LO topology. UV and IR singularities are present at this order and must be appropriately renormalized and regulated. We describe the calculation of the various UV-renormalized matrix elements for each phase-space configuration in ref. [55] for the decay H → bbj at NNLO. This leaves the discussion of the IR regulation, which is different from that described in ref. [55] due to the LHC kinematics.
In order to regulate the IR divergences present at this order we use the N -jettiness slicing approach [56,57]. This method has become an established technique for evaluating NNLO cross sections involving final-state jets at the LHC [51,[57][58][59], and we provide a brief overview in this section. The central idea is to separate the (differential) cross section of a process into two pieces, where the variable τ N is the N -jettiness variable [60]. For our 1-jet example, this variable is defined as where {p m } is the set of all partonic momenta in an event, while {k i } are the momenta of the two incoming beams and the hardest jet present in the event (after clustering). The quantity P i is a somewhat arbitrary choice of hard scale, and in our calculation we take P i = 2E i (known as the geometric measure [61,62]). The above-cut term σ(τ N > τ cut n j ) has sufficiently large value of the N -jettiness variable to have at most one unresolved parton, and therefore corresponds to a NLO computation of the cross section with an additional parton present. The below-cut term σ(τ N ≤ τ cut n j ) contains all of the double-unresolved limits at NNLO. However, in the limit τ cut 1 → 0 the cross section can be approximated using the following factorization theorem from Soft-Collinear Effective Field Theory (SCET): 3)

JHEP05(2021)045
where in our case n j = 1. The above equation is valid up to power corrections (denoted by the F (τ cut n j ) term), which vanish in the limit τ cut n j → 0. At NLO the leading power corrections are well described by the form τ cut n j log(τ cut n j /Q), and at NNLO the leading power corrections have the form τ cut n j log 3 (τ cut n j /Q) (where in both cases Q is the hard scale associated with the process). The general terms that enter the SCET factorization theorem are the soft (S), jet (J ), and beam (B) functions, for which calculations accurate to O(α 2 s ) needed for our calculation can be found in refs. [63][64][65][66][67][68].
There are several alternative choices [51,59] one can make when applying the jettinessslicing method. Firstly, one can choose whether to work with a fixed version of τ cut 1 , in which all events are compared to a given energy scale, or with a dynamical definition, in which the final-state kinematics (of the clustered system) generates different τ cut 1 values for each phase-space point. Typically, for 1-jet NNLO processes it is more prudent to use the latter option. Since power corrections are sensitive to the overall hardness of the system through the expansion parameter τ cut 1 /Q, very energetic jets have suppressed power corrections. By using a fixed τ cut 1 , the calculation for these terms includes points that are very soft and collinear (relative to the hard scale), resulting in large Monte Carlo uncertainties and code instabilities. On the other hand, using a dynamic τ cut 1 ensures a more relaxed τ cut 1 for more energetic jets, reducing this problem and producing more stable results, without increasing the impact of unwanted power corrections.
In order to obtain the remaining process-specific hard function (H) appearing in eq. (2.3), we use our double-virtual calculation for the decay amplitude H → bbg presented in ref. [55]. 1 The result for LHC kinematics is obtained by performing the relevant crossing, moving the desired final-state partons to the initial state. In practical terms, this involves taking the appropriate analytic continuation of the various harmonic polylogarithms that appear in the virtual amplitudes as described in section 4 of ref. [70]. After crossing the relevant final-state partons to the initial state, we have checked that our results have the correct factorization properties in the relevant soft and collinear limits [71,72], finding excellent agreement.

Matching to the EFT
The Standard Model does not allow for the consideration of the impact of a single fermion generation in isolation. For the purposes of this calculation, in order to completely specify our theoretical framework we must also address the role of the top quark in the computation. This is because at O(α 3 s ) the cross section becomes sensitive to the presence of the top induced production. One must therefore specify whether one works in the effective field theory or full Standard Model. Precision calculations in the full Standard Model are made considerably more difficult by the presence of the additional mass scale and are currently known to NLO accuracy for H +j [73]. On the other hand, NNLO predictions are available in the EFT [45][46][47][48][49][50][51]. Typically, in EFT calculations the top mass effects are included via a rescaling of the cross section by those computed in the full theory at lower orders.

JHEP05(2021)045
For H + j in the 5FS at O(α 3 s ) accuracy there are two contributions, the pure bottominduced and the LO top-induced piece, and since the top-induced contribution is leading order one could easily consistently work in the full SM or the EFT. However, it is known that the higher corrections to the EFT pieces are large, and therefore including only the LO piece makes little phenomenological sense. Our strategy in this paper is to ignore the topinduced pieces altogether and focus only on the technical aspects of the NNLO calculation of the bottom-induced contribution. In order to obtain reliable phenomenological predictions at "NNLO", one would therefore wish to combine the O(α 3 s y 2 b ) pieces with the O(α 5 s ) EFT results. To avoid having the bottom-induced component be entirely overwhelmed by the EFT piece, one would also wish to apply b-tagging requirements (and consider other sources of Higgs plus heavy flavor [31] arising from VBF and V H processes). We postpone such detailed phenomenology study to a future publication. Additionally, we note that when working in the EFT the bottom Yukawa coefficient is matched to that of the full SM as follows, This means that, when working in the context of the EFT, we should include a term proportional to ∆ F multiplying our LO predictions. In this paper we remain agnostic to the exact implementation of the top quark and therefore choose to present results in terms of the unmatched y SM b . The impact of adjusting the coefficient to y EFT b is a rather small (sub-percentage) effect and does not affect the conclusions presented in this paper.
Finally, we note that there are interference terms between the top quark (or EFT) initiated contributions to H + j and the bottom-induced contributions. This interference requires a helicity flip in order to be non-zero, inducing an overall scaling of the form y MS b m OS b (since the helicity flip is a kinematic mass). As a result, the interference vanishes in the 5FS. However, there is an ambiguity in the mixed-renormalization scheme which renders this argument not quite complete, since one can relate the OS mass to the MS mass changing the scaling to y , such a procedure is well-defined since the interchange of the mass schemes is trivial. However, at higher orders this procedure is much more delicate due to the presence of IR logarithms in m OS b , and rich UV structure. Very recently, this limit was studied in the context of extracting a sensible result at NLO in the interference for Higgs-plus-charm production [74] (where the even larger hierarchy between y c and y t makes these terms more important). In addition, these pieces were studied in ref. [75] for the decay of H → bb and H → cc at O(α 3 s ) in the "4FS" in which the mass was fully retained. In keeping our focus on the technical aspects of the NNLO computation, and being agnostic regarding the top quark implementation, in this paper we take the first limit, in which the interference is set to zero. However, when pursuing a full LHC phenomenology study and given the size of y t , we advocate including the term using JHEP05(2021)045 the limit extraction in which the helicity flip mass is coverted into the MS mass prior to taking the limit. We leave this to a more detailed future study.

Validation
The calculation described in the previous section has been implemented into the Monte Carlo code MCFM [76][77][78][79]. We make extensive use of the code's ability to handle processes involving a final-state jet at NNLO, and particularly important for this paper are the MCFM developments outlined in refs. [51,59]. This section details the various checks we have performed on our computation (in addition to the analytic soft and collinear checks previously mentioned).

bb → H at NNLO
We validate our calculation of bb → H at NNLO by comparing our results to those known in the literature. This process has been well studied and public codes are available for the computation of cross sections at NNLO accuracy. We use the SusHi framework [22,80], which can compute a variety of Higgs production cross sections at NNLO accuracy in the SM and its supersymmetric extensions. In this comparison we run MCFM with parameters set to match the default implementation in SusHi. We use the MMHT14 [81] PDF sets, (taking the NNLO set for all predictions) and the following setup for our comparison: As discussed in previous sections, the choice of a factorization scale around the Higgs boson mass is somewhat of a problem for phenomenology, since the perturbation theory is subject to large corrections. However, in this case the large corrections act in our favor when attempting to validate our implementation, since the larger NNLO coefficient allows us to separate scales associated with the pure coefficient, power-suppressed corrections, and numerical Monte Carlo uncertainties. With the parameter choices listed above, the values of the bottom quark mass used in the Yukawa couplings are m MS b (µ R ) = {2.868, 2.666, 2.645} GeV at LO, NLO, and NNLO respectively. By matching the MCFM parameters to these values we observe excellent agreement at NLO: the result from Sushi is 701.97 fb, while MCFM gives 701.95 fb. We then proceed to compare our result for the NNLO coefficient in figure 2.
Here we use the τ 0 parameter in the laboratory (unboosted) frame and present results over the range 0.001 ≤ τ cut 0 ≤ 0.05 GeV. As is by now well known in the literature, the leading power corrections at NNLO can be described parametrically as follows, where the ellipses indicate sub-leading contributions of the form τ log 2 τ etc., and Q is a hard scale associated with the process (e.g. m H ). The residual power corrections in our results are clearly well described by this parametric form, and by fitting our results accordingly we are able to simultaneously extract the coefficient in the limit τ cut 0 → 0 and JHEP05(2021)045 which is in excellent agreement with the coefficient obtained from Sushi, δ NNLO Sushi = −100.14 fb. Our results clearly show that for τ cut 0 in the region 10 −3 GeV the residual power corrections are significantly less than 1% of the NNLO coefficient, and subsequently per-mille level relative to the total physical prediction. In addition to the detailed comparison described above, we have performed a similar fit to our calculation with the canonical scale choice of µ R = m H = 4 µ F , obtaining δ NNLO 0 = 18.39 ± 0.18 fb, which is again in excellent agreement with the result obtained from SusHi, 18.52 fb.

bb → H + j at NLO
Before studying the slicing dependence of the main result of our paper, the bottom Yukawa contributions to H +j at NNLO, we study the process at NLO. The primary area of interest is to study the different options for the definition of τ cut 1 and their associated asymptotic regions of validity. In order to test the various ingredients of our calculation, we begin by computing cross sections for H + j at NLO using the different IR-regulating prescriptions described in the previous sections. For these comparisons we use the following setup: with jets clustered using the anti-k T algorithm with R = 0.4. Additionally, we will briefly study the power corrections with the more central jet requirement of |η j | < 2. 5 Next, we turn our attention to validating the NLO cross section, which we have computed using dipole subtraction [83] and the jettiness-slicing approach. As part of the validation of the dipole calculation we have checked the (in)dependence of our result on the unphysical α parameters [84], which control the amount of non-singular phase space utilized in the dipole subtractions. Using the dipole method and the parameter choices above, the corresponding NLO cross section is 144.98 fb. We now consider the various implementations of the jettiness-slicing method. Our results are presented in figure 3, where the panels on the left side show the ratio of the cross section obtained using a fixed value for τ cut 1 to the dipole result, for both the boosted and traditional definitions. The data points on the figure show the results obtained with the full phase-space cuts described above as well as a fit to the data of the form In order to quantify the impact of forward radiation on the power corrections we additionally show a fit to similar results obtained with a tighter jet requirement |η j | < 2.5, although for readability we suppress the Monte Carlo output. The difference between the solid and dashed lines on the figure is therefore indicative of the sensitivity of the power-suppressed terms to the presence of forward jets. The panels on the right side of figure 3 show the same cross section ratios, computed using a dynamic version of τ cut 1 , which we define as As before, results are evaluated in both the laboratory and boosted frames. The corresponding fit for this setup is as follows, We observe the same pattern for the impact of the power corrections as reported in ref. [51]. By evaluating in the boosted frame, the size of the power corrections is significantly reduced [85], especially when the phase space includes contributions from regions in which the jet has larger pseudo-rapidity (|η j | > 2.5). Using the dynamic version of τ cut 1 also results in smaller power corrections, and in particular the boosted-dynamic definition is the least sensitive to power corrections. We therefore employ the boosted-dynamic version of the slicing in our subsequent studies at NNLO.

bb → H + j at NNLO
In this section we discuss the validation of our primary result, the NNLO predictions for the bottom-induced contributions to H + j. We begin by presenting a check of the H + 2j NLO result, which forms the above-cut piece of our NNLO prediction. As before, we check the α (in)dependence of the calculation, for which results are presented in figure 4. These predictions were obtained using the NNLO CT14 PDF set [86], µ R = µ F = m H , and the two-loop running of the bottom Yukawa coefficient. By comparing results at sub per-mille level accuracy we are able to rigorously test the cancellation of IR singularities at one loop and the subsequent cancellation of dipole-related terms from the real-virtual and double-real contributions at NNLO. Following the notation of ref. [51], we define ab = σ(α ab = 1) − σ(α ab = 10 −2 ) σ(α ab = 1) , (3.6) where the indices a and b correspond to either initial-(I) or final-(F ) state dipoles. Figure 4 illustrates that our results are insensitive to the choice of the α parameter at the level of ∼ 0.0005 for the qg, qg, and qq channels. We have studied these channels in greater detail since they receive contributions from both the ggbb amplitudes and the four-quark amplitudes, and therefore have the most intricate IR structure. Results are also shown for gg and qq fluxes, which we have constrained to the (still stringent) level of ∼ 0.001, ∼ 0.002, or better. We are therefore confident that the cancellation between the unintegrated and integrated dipoles has been correctly implemented in our NLO H +2j calculation.
Finally, we arrive at the main result of this section, namely the validation of the τ cut  Our results are presented in figure 5. We study the dependence of the NNLO coefficient on the dynamic version of τ cut 1 , which we recall is defined as

JHEP05(2021)045
The results of ref. [51] for H + j in gluon fusion, and our preceding study for this process at NLO, clearly demonstrate that this choice results in the smallest power corrections, particularly when the associated jet is not required to be central (as in our case). The results in figure 5 span the range 2×10 −5 ≤ ≤ 1×10 −3 , which is approximately equivalent to a fixed τ cut 1 in the range 0.0025 − 0.12 GeV (setting p H T = 30 GeV, which corresponds to the minimum jet transverse momentum). As expected, our results are well described by the following approximation for the power corrections, where δ NNLO 0 represents the physical correction obtained in the limit → 0. We find  Results are presented for the dynamic version of τ cut 1 , evaluated in the rest frame of the H + j system. Also shown is a fit to the results parametrizing the residual (τ cut 1 ) power corrections. The dashed line corresponds to the limit → 0 of the fit, and the shaded band represents fitting uncertainties on the asymptotic limit. which means that we are able to control the remaining unknown power corrections to the level of 0.2% on the NNLO cross section. For the remainder of this paper we will use the boosted dynamic τ cut 1 with set to 1 × 10 −4 . From our preceding study we can estimate that for this value the remaining power corrections should be at the level of a few percent on the NNLO coefficient, and hence around the per-mille level on the full physical NNLO prediction. Such a level of accuracy should be more than adequate for the phenomenology presented in the next section.

Results
In this section we present our results for the NNLO predictions for H + j at the 13 TeV LHC. We use the CT14 PDF sets [86], matched to the appropriate order in perturbation theory (the running of α S and m MS b (µ R ) therefore occurs at the next perturbative order). We use the same fiducial cuts as in section 3, namely we cluster jets using the anti-k T algorithm with R = 0.4 and require them to have p j T > 30 GeV and |η j | < 4.5.

Factorization and renormalization scale dependence
We begin by investigating the factorization and renormalization scale dependence of the total cross section for H + 1j at NNLO accuracy. Additionally, we also investigate the cross section for H + 0j production at the same order. Although higher-order predictions are now available [33], it is nevertheless interesting to compare the two predictions with and without the additional jet requirement at NNLO. Our results for H+0j and H+1j are shown in figure 6. We set a central renormalization scale of µ R = m H and vary the factorization scale in the range m H /16 ≤ µ F ≤ 4m H . As is well known, the factorization scale dependence for the H + 0j cross section is dramatic at LO and NLO, while at NNLO the behavior is somewhat improved (and even more so at N3LO [33]). The maximum value around µ F = m H /3 adds weight to the historical argument of using µ F = m H /4 as the central scale choice in NLO predictions [16].

JHEP05(2021)045
The presence of initial-state gluons and a final-state jet conspire to decrease the dependence of the H + 1j cross section on the factorization scale when compared to the equivalent H + 0j result. However, it is clear that the H + j cross section still bears a striking dependence on the factorization scale. At LO, across the range studied the cross section increases by a factor of 7, while at NLO the increase is a factor of two. The NNLO results from our calculation lead to a substantial improvement. By including second-order corrections, the overall increase in the cross section over the range of µ F drops to a factor of 1.36. Indeed, the vast majority of this increase occurs at lower scale choices, while at larger scales m H ≤ µ F ≤ 4m H the NNLO cross section changes only by around 6% (compared to a change of 36% at NLO over the same range of µ F ). It is therefore clear that NNLO accuracy is mandated for a robust estimate of rates, free from large uncertainties induced by the unphysical dependence on the factorization scale. As the dependence on µ F significantly drops for µ F ≥ m H /4, we will choose a central scale choice that reflects this in our subsequent predictions.
In figure 7 we turn our attention to the renormalization scale dependence of the NNLO H + 0 j and H + 1 j cross sections. For the bb → H process we make the customary choice µ F = m H /4 and, motivated by the results discussed in the preceding paragraphs, we choose µ F = m H /2 for the H + 1j predictions. We present the dependence of the cross JHEP05(2021)045 sections over the range m H /8 ≤ µ R ≤ 8m H . For the bb → H process, at leading order the only dependence of the cross section on µ R arises from the evolution of m MS b . Higherorder corrections induce a rather mild dependence through α s , and the NNLO prediction is already rather insensitive to the renormalization scale. For the H + j cross section the situation is rather different, since here the LO result depends on both α s and m MS b . On the smaller end of the considered range µ R < m H /2, the perturbation theory becomes rather unreliable, with large corrections at each subsequent order. However, in the region m H ≤ µ R ≤ 8m H the perturbation theory becomes well-behaved. There is also a significant reduction in the residual scale uncertainty at NNLO: σ pp→H+j (µ R = m H )/σ pp→H+j (µ R = 8m H ) is ∼ 1.37 at NLO, but reduces to ∼ 1.11 at NNLO.
The results of this subsection indicate that predictions obtained using central scale choices comparable to (µ R , µ F ) = (m H , m H /2) should demonstrate convergent behavior in the perturbative expansion, with a reasonably small residual dependence for excursions from this central choice. We therefore choose these values for the differential predictions presented in the next section.  increased and decreased in opposite directions). For each bin in the distribution the largest deviations from the central value are taken as upper and lower estimates of the scale variation. Figure 8 shows the results for the rapidity distribution of the Higgs boson (y H ) and leading jet (y j ) at LO, NLO, and NNLO accuracy. By comparing the two distributions it is clear that the Higgs boson is produced with a more central distribution than the jet, which has a broader distribution (and hence more forward jets). This can be traced back to the underlying kinematics, since the Higgs boson is a scalar particle and therefore is produced more isotropically, while the leading jet favours the collinear (forward) region in which the quark-gluon splitting is enhanced. The pattern of higher-order corrections is broadly similar for both distributions, with a significant shape change from LO to NLO and a much smaller change from NLO to NNLO. We observe that the scale variation decreases from around ±10% at NLO to around ±4% at NNLO. In the central region |y X | < 2.5 (X = H, j), the corrections are relatively flat, whereas in the larger rapidity regions they become more sizable. We note that some care should be taken in this region, since it could be prone to larger power corrections in the N -jettiness slicing method. However, we estimate that remaining power corrections enter at around the percent level in the tails JHEP05(2021)045  of the two distributions, which should not substantially change the interpretation of the plots. In the lower panel we additionally include an estimate of the uncertainties due to the PDF extractions, obtained at 68% C.L. using LHAPDF [87]. We see that at NNLO the uncertainties from the PDFs and the six-point scale variation are of the same order (∼ 5%). For very forward Higgs bosons the PDF uncertainties become very large, but there is very little cross section in this region. Figure 9 presents the transverse momentum distribution of the Higgs boson and the leading jet. For the transverse momentum of the Higgs boson the softest bin p H T < 30 GeV corresponds to an observable one perturbative order lower than the rest of the calculation (since there exists no 2 → 2 underlying topology) and this is reflected in the larger NNLO/NLO ratio and overall scale variation of this bin. Focusing on the change from NLO to NNLO, we see that the ratios are rather dynamic (while remaining within the uncertainty band of the NLO calculations), especially for the transverse momentum of the Higgs boson. In the region p H T < m H /2 there is a decrease in the prediction of around 5 − 10%, while in the region around p H T ∼ m H the prediction increases at NNLO by around 10 − 15% before asymptotically approaching a smaller ratio (around 5%) at larger trans-JHEP05(2021)045 verse momentum. The leading-jet transverse momentum distribution is different: the main impact of the NNLO corrections is to soften the spectrum, especially at high p T . Finally, for both distributions the PDF uncertainties are again comparable in size to those obtained using a six-point scale variation.

Conclusions
We have presented a NNLO calculation of the bottom-induced contributions to pp → H +j in the 5FS. Our results have been implemented into MCFM and, as a useful by-product and cross-check of our setup, we have also implemented Higgs boson production through bottom-quark fusion bb → H at the same order. Analytic results for the various Higgsplus-parton amplitudes needed in this paper have been obtained from our previous study of the Higgs boson decay to three partons, crossed to the appropriate LHC kinematic configurations. We have performed many cross-checks of our calculation. At the analytic level we have checked the collinear and soft factorization properties of our two-loop amplitudes, and have carried out numerical checks of our H + 4, 5 parton amplitudes. For the H + 2j NLO computation we have done extensive testing on the exact cancellation of the integrated and unintegrated dipole subtractions. We have compared our results for bb → H to the public code Sushi, finding excellent agreement.
Higgs plus jet is fast becoming a standard candle process to study the Higgs boson at the LHC. While the bottom-initiated contributions remain a small piece of the total cross section, their study is motivated by the strong desire to constrain the bottom Yukawa coupling wherever possible. Unfortunately, bottom-induced processes in the 5FS have a strong dependence on the factorization scale, making a precision computation of these cross sections challenging. The results of this paper show a dramatic stabilization of the cross section at NNLO, particularly with regard to the overall dependence on the factorization scale.
While our results focused on the process pp → H + j, there are several interesting extensions of this study to pursue. Firstly, in order to target the bottom Yukawa interaction, experimental analyses would need to impose b-tagging requirements. Therefore, adjusting our computation to target pp → H + b at NNLO is a logical next step. In order to do a meaningful phenomenological study for the LHC, the bottom-induced contributions should be compared to the dominant production mechanism through gluon fusion (with a b-tagging requirement applied). Dealing with the interference term in a sensible limit in the 5FS at this order is also interesting for a full phenomenology study. The inclusion of the decay of the Higgs boson is also a must in order to adequately describe fiducial volumes, for which a particularly interesting example is when the Higgs boson decays to bottom quarks (which must also be matched to NNLO accuracy). With these alterations in hand, it would also be interesting to study Higgs-plus-charm production. Finally, in addition to the SM Higgs boson, the calculation described here could also be used in BSM extensions. For instance, the Higgs-bottom quark coupling could be modified as in the MSSM or SMEFT, or the mass of the scalar particle itself could be increased, e.g. in dark matter searches where the scalar particle acts as a mediator of putative dark forces (with potential missing energy decays). We leave these exciting and detailed studies to future work.