NNLO+PS predictions for Higgs production through bottom-quark annihilation with MINNLO$_{\text{PS}}$

We consider Higgs production through bottom-quark annihilation at hadron colliders and we calculate next-to-next-to-leading-order (NNLO) corrections in QCD perturbation theory matched to parton showers (NNLO+PS). To this end, we have adapted the MINNLO$_{\text{PS}}$ method to account for the extra scale dependence induced by an overall Yukawa coupling that is $\overline{\rm MS}$ renormalized. We compare our results against state-of-the-art fixed-order predictions at NNLO as well as resummed predictions at next-to-next-to-leading-logarithmic (NNLL) accuracy.


Introduction
The Higgs boson provides one of the cornerstones of the Standard Model (SM) of particle physics.With its discovery about a decade ago [1,2] the measurement of the properties of the Higgs boson has become a major quest in the rich physics programme at the Large Hadron Collider (LHC).The exploration of the Higgs sector is of utmost importance not only in the context of the SM, but also in the search for new-physics phenomena.So far the characterization of the Higgs coupling to top (t) and bottom (b) quarks, W and Z bosons, and tau leptons is fully consistent with the SM picture [3,4].However, the coupling measurements are becoming continuously more precise with the rapid data taking by the LHC experiments, becoming increasingly more sensitive to small deviations from the SM expectations.At the same time, other couplings that are currently restricted due to large statistical uncertainties are anticipated to become more accessible in the future.An example of this is the self-interaction of the Higgs boson.
The accurate simulation of all relevant Higgs-production and decay modes at the LHC is henceforth a crucial requirement for successfully finding deviations from the SM predictions.In this context, the associated production of a Higgs boson with bottom quarks (b bH) plays a special role.Although its total rate (of less than 1 pb) is only about two percent of the dominant Higgs production mode through gluon-fusion, it is still large enough so that its cross section has to be accounted for in precision measurements of the Higgs boson at the LHC.On the other hand, the direct detection of a b bH signal (by tagging the bottom quarks) is extremely challenging at the LHC, due to large backgrounds and because its rate is substantially reduced by requiring one or two bottom quark tags.Even if the measurement of b bH production would be possible, its main purpose of extracting a bottom-quark Yukawa coupling is essentially hopeless due to a large contamination from other production mechanisms that do not contain the bottom-Yukawa coupling, see e.g.Ref. [5].Nevertheless, b bH production (and its precise simulation) is particularly relevant in two further respects.Firstly, it is the main Higgs production mechanism in beyond-the-Standard-Model (BSM) theories with enhanced bottom-Yukawa coupling, for instance for the production of a heavy Higgs boson in a Two-Higgs-Doublet-Model (2HDM) like the Minimal-Supersymmetric-SM (MSSM) with a large value of tan β.More importantly, b bH production is the dominant irreducible background in searches for Higgs-boson pair (HH) production in the SM in the most sensitive search channels where at least one Higgs boson decays to bottom quarks (for a review, see Ref. [6]).A reliable modeling of the b bH background to HH measurements becomes indispensable at the High-Luminosity phase of the LHC (HL-LHC), when the Higgs-pair production cross section in the SM is expected to be measured with a significance of 3.4σ (4.9σ) [7] in the combination of all search channels. 1  The dominant contributions to the b bH process are those proportional to the bottom Yukawa coupling (y b ) where the Higgs couples to a bottom-quark line, see Figure 1 (a) and (b), as well as those proportional to the top Yukawa coupling (y t ) where the Higgs boson couples to a closed top-quark loop, see Figure 1 (c).In fact, the latter, which corresponds to the gluon-fusion process with a b b pair originating from a QCD splitting, has a slightly larger cross-section yield, and its relative size further increases when tagging the two bottom quarks.Nevertheless, both production mechanisms are relevant in the SM and they receive particularly large QCD perturbative corrections, so that higher-order calculations for these processes are crucial.Other b bH production modes through V H associated production with V → b b and b-associated vector-boson fusion have a subleading impact on the cross section from a few percent up to a few tenths of percents depending on the category, and simulations for these production channels exist [5].
Different schemes can be employed for the calculation of the b bH process, since the bottom quark can be considered both a massless or a massive quark at the typical scale of the b bH production at the LHC.Therefore various predictions have been obtained in both a five-flavour scheme (5FS) with massless bottom quarks, see Figure 1  progress in higher-order calculations has been made for the y 2 b contribution due to the much more involved structure of the LO process [5,22,23,[28][29][30][31][32][33][34].Next-to-LO (NLO) corrections in QCD (matched to parton showers) combined with NLO electroweak (EW) corrections are still the state-of-the-art.On the other hand, a combined study of the y 2 b and y 2 t contributions to b bH production modes has been performed only in the 4FS [35], which included NLO QCD corrections in all relevant coupling structures (y 2 b , y 2 t and y b y t interference contributions).The matching to parton showers of this full NLO QCD calculation has been considered in Ref. [36] and studied in the context of backgrounds to HH searches.Differences between 4FS and 5FS results have been studied in various works, see e.g.Refs.[37,38], and consistent combinations of the two schemes have been obtained in Refs.[39][40][41][42][43].
In this paper, we focus on the 5FS calculation of the b bH process proportional to y 2 b and perform the first fully-differential calculation of next-to-NLO (NNLO) QCD corrections matched to parton showers (NNLO+PS).To this end, we exploit the MiNNLO PS method for colour-singlet production presented in Refs.[44,45] and we adapt the method such that it can account for an overall scale-dependent Yukawa coupling renormalized in the MS scheme.We compare our predictions against reference predictions computed at NNLO in QCD [10,46] and against the resummed transverse-momentum (p T ) distribution at NNLL+NNLO accuracy [17].

Outline of the calculation
We consider the production of a Higgs boson through bottom quark annihilation in the 5FS, where the bottom quarks are treated as being massless, as given in Figure 1 (a) for the LO process, which is proportional to the bottom-quark Yukawa coupling y b .From relative order α 2 s on, i.e.NNLO, a loop-induced gluon-fusion contribution proportional to the top-quark Yukawa coupling y t enters, see Figure 1 (c).We refrain from including this contribution throughout this paper, as it can be considered as a completely independent process, which can be obtained via higher-order simulation of Higgs-boson production in gluon fusion.Moreover, in the 5FS, any y b y t interference contribution between the bottom-quark annihilation and gluon-fusion Higgs production processes vanish to all orders in perturbation theory, since it is proportional to the bottom-quark mass.
We implement a fully differential computation of Higgs production in bottom-quark annihilation in the 5FS up to NNLO in QCD perturbation theory and consistently match it to a partonshower simulation.To this end, we have adapted the MiNNLO PS method for colour-singlet production [44,45] to account for an overall scale-dependent Yukawa coupling renormalized in the MS scheme.And we have implemented an alternative scale setting in all contributions up to NNLO QCD that are regular in the transverse momentum of the Higgs boson (p T ).Both adaptations of the MiNNLO PS method are described in detail in Section 3.
The MiNNLO PS method has various positive features.In particular, it provides physically sound results without relying on an unphysical slicing scale to separate events with different jet multiplicities.It is numerically very efficient, since the NNLO corrections are included directly in the event generation without any post-processing or reweighting of the events.And, not least, it preserves the leading-logarithmic accuracy of the parton shower by keeping the appropriate ordering of the emissions and scaling intact.The MiNNLO PS was not only extended beyond 2 → 1 production processes in Ref. [47] and applied to several colour-singlet production processes in Refs.[47][48][49][50][51][52][53][54][55], it was even reformulated for the case of heavy quark-pair production [56][57][58], being the first method applicable to processes with colour charges in both initial and final state thus far.
Our MiNNLO PS b b → H generator has been implemented within the Powheg-Box-Res framework [59].First, we have implemented a NLO+PS generator for Higgs plus jet production in bottom-quark annihilation using the Powheg method [60][61][62], see Figure 2 for examples of respective LO diagrams.For the evaluation of the tree-level amplitudes of the Higgs plus jet (HJ) and Higgs plus two-jet (HJJ) processes we employ OpenLoops [63][64][65], using its interface Powheg-Box-Res developed in Ref. [66].For the one-loop virtual corrections to the HJ processes we have used the analytic results from Ref. [12], which substantially improve the numerical performance of the code.We have cross-checked this implementation numerically against OpenLoops for several phase space points, finding full agreement at the level of the machine precision.In a second step, we have extended the HJ NLO+PS implementation to NNLO accuracy for b b → H production through the newly extended MiNNLO PS method described in the next section.
3 Revising the MINNLO PS method

Original method
In the following, we briefly summarize the MiNNLO PS formalism for colour-singlet production, which has been introduced in Ref. [44] and optimized in Ref. [45].The MiNNLO PS cross section for b b → H production can be expressed through the standard Powheg formula for HJ production with a modified content of the Powheg B function [60][61][62]: where Φ HJ is the HJ phase space.In the above expression, ∆ pwg is the Powheg Sudakov form factor with a cutoff Λ pwg = 0.89 GeV, while dΦ rad and p T,rad are the phase space measure and the transverse momentum of the real radiation with respect to HJ production.In this context, Powheg takes care of the matching of fixed-order HJ calculation with a parton shower by producing the first additional radiation through the ratio of the tree-level matrix elements for HJJ (R HJ ) and HJ productions (B HJ ).All the subsequent radiations with smaller transverse momenta are generated by a Shower Monte Carlo.
The key ingredient of the MiNNLO PS method is the modified Powheg B function which is denoted as BMiNNLO PS in Eq. ( 2).The derivation of this function stems from the following formula for the NNLO cross section of Higgs production differential in the Born phase space (Φ H ) and in the Higgs transverse momentum (p T ): In the above equation, the cross section is divided into the singular component and the regular part R f (p T ), which is finite in the p T → 0 limit.The luminosity L(p T ) includes the squared virtual matrix elements for b b → H production and the convolution of the parton densities with the collinear coefficient functions ( Cij ), which can be written as while the exponent of the Sudakov form factor is defined through where M is the invariant mass of the colour-singlet final state, i.e.M = m H for Higgs production in bottom-quark annihilation.Notice that the hard function ( H) and the parton densities are computed at the scale p T , which is a feature of the MiNNLO PS approach.In bottom-quark annihilation the Higgs is produced either via b b or bb initial state partons, which is why the index c denotes either a bottom quark or antiquark.Since the Born and virtual matrix elements for this process are invariant under charge conjugation of the initial-state bottom quarks, we will drop the index of |M (0) cc |2 and Hcc in the following.The expansion of the relevant resummation coefficients and coefficient functions in the MiNNLO PS scheme is given by with 2 where P (0) ij (z) is the leading-order regularised splitting function and β 0 = (33 − 2n f )/(12π), with n f being the number of light quark flavours.The latter replacements of the resummation coefficients in the MiNNLO PS formalism (marked by the tilde symbol) have been derived in detail in Ref. [44].They originate from the translation from b-space to direct space, while the H (1) term in B( 2) is due to the evaluation of the hard function at the scale p T .The coefficients without tilde are the standard ones for transverse-momentum resummation production for quark-induced processes, which have been summarized in Ref. [44].In particular, for quark-initiated processes C (2) has been obtained in Ref. [67], and the resummation coefficients A (1,2,3) and B (1,2) can be found in Ref. [68].
The hard-virtual coefficients H (1) and H (2) are determined from the one-loop and two-loop amplitudes for b b → H production, which have been computed for the first time in Ref. [69] and Refs.[70,71], respectively.In the MiNNLO PS resummation scheme they are expressed as where the Casimir factors for SU(3) are C F = 4 3 and C A = 3.We now return to our starting formula in Eq. ( 3) to derive the MiNNLO PS master formula.The regular part R f (p T ) can be written as where the first term on the right-hand side of the above equation is the NLO differential cross section for the production of Higgs boson in association with one jet (HJ) in bottom-quark annihilation, given by (2) The notation X(µ) indicates that the quantity X is evaluated at the scales µ R = µ F = µ.Moreover, we can express the coefficients of the expansion of the singular terms as dσ We can now rewrite Eq. ( 3) by factoring out the Sudakov exponential and then substitute R f (p T ) by its expression in Eq. ( 15) (2) Factoring out the Sudakov exponential entails two positive features of the MiNNLO PS method.Firstly, it improves the numerical stability when integrating over the low-transverse momentum region.Secondly, it removes the need for a slicing cutoff at small transverse momenta.
Given that we have modified our NNLO accurate starting formula in Eq. ( 3) only by terms beyond accuracy, our MiNNLO PS master formula in Eq. ( 19) includes NNLO accuracy (upon integration over p T ) by construction.Notice that up to α 2 s Eq. ( 19) corresponds exactly to the MiNLO ′ formula [72], which is NLO accurate in both the H and the HJ phase spaces, while adding the relevant singular terms α 3 s (and beyond), which are generated by the total derivative in our starting equation ( 3), adds the relevant contributions required to reach NNLO accuracy in the H phase space upon integration over p T .
Finally, we apply the same concept to render the Powheg B function in Eq. ( 2) NNLO accurate, by deriving it as (2) where the function F corr encodes spreading of the MiNNLO PS corrections in the full Φ HJ phase space, given that the function D depends only on the kinematical variables of the colour singlet.
At last, we note that the logarithmic terms contained in S(p T ) and D(p T ) are switched off at p T ∼ M with the following replacement, and taking into account the required jacobian factors.In our calculation, we will use these modified logarithms with p = 6 in order to match smoothly with the fixed order prediction at high transverse momenta.

Accounting for an overall Yukawa coupling
We now determine the scale dependence of the MiNNLO PS formulae and coefficients in the presence of an overall MS-renormalized Yukawa coupling.In Appendix D of Ref. [44], the scale dependence of the original MiNNLO PS formulation was presented for processes with an overall power of the strong coupling at Born level.In the case of Higgs production in bottom-quark annihilation in the 5FS, the cross section instead involves two powers of the bottom-Yukawa coupling at Born level.
For the sake of generality, we present all relevant formulae for a process with the LO coupling structure i.e. involving both an arbitrary overall power n B of the strong coupling and an arbitrary overall power m B of the bottom-Yukawa coupling at Born level.This is the case for instance for the 4FS process of Higgs production in association with bottom quarks (where n B = m B = 2).The bottom-quark Yukawa can be expressed through with m b being the mass of the bottom quark and v being the vacuum expectation value of the Higgs field.Below, we present all formulae needed to implement (independent) scale variations in the strong and the Yukawa couplings.
Different renormalization schemes can be employed for the bottom-quark mass in the Yukawa coupling.Given that the natural scale of the Yukawa coupling is a hard scale (for instance of the order of the mass of the Higgs boson), the most appropriate choice is the MS renormalization scheme, which introduces a renormalization scale for Yukawa coupling (or more precisely its mass) that can be set to a suitable scale.In the MS scheme, the scale dependence of Yukawa coupling can be deduced by solving the renormalization group equation (RGE) for the bottom-quark mass, which reads Here, the anomalous dimension can be expressed as an expansion in α s (µ) where we recall that n f denotes the number of light quark flavours.For a generic power m B of the Yukawa coupling at Born level, we derive the relevant scale-compensating terms up to NNLO when performing a change of the Yukawa renormalization scale from M to µ (0),y R by solving Eq. ( 24) In the above equation, we introduced an arbitrary scale µ (0),y R for the choice of Yukawa coupling, where the index (0) indicates that this scale is related to the Born cross section.It represents a generic scale that can be varied in order to study the theoretical scale uncertainty related to the Yukawa coupling.
Accordingly, changing the scale for an arbitrary power n B of the strong couplings at Born level from M to µ (0),α R takes the following form: with Here, we introduced an arbitrary scale µ (0),α R for the strong coupling constant in the Born cross section.Also in this case, µ (0),α R can be varied in order to probe the uncertainties related to missing higher-order terms.
Within MiNNLO PS , the scale-compensating terms originating from the variation of the overall Born couplings are implemented at the level of the hard-virtual coefficient function in Eq. ( 6).Introducing explicitly the scales µ (0),α R and µ (0),y R , the squared virtual matrix element can be written as Here, M (0) (µ ) is the tree-level amplitude with the strong and Yukawa couplings evaluated at µ (0),α R and µ (0),y R , respectively.In addition, we have introduced a generic symbol µ R for the renormalization scale of the extra powers of the strong coupling in the expansion of H, which is set to µ R = K R p T following the MiNNLO PS prescription.Using the identity while inserting the relations in Eqs. ( 26) and ( 27), we can absorb the logarithmic scale-compensating terms into the hard-virtual coefficient function For Higgs production in bottom-quark annihilation the coefficients β 0,1 and γ 1,2 are calculated considering five massless quarks (n f = 5), and there is no strong coupling (n B = 0), but two powers of Yukawa couplings (m B = 2) at the Born level.Moreover, the invariant mass M corresponds to m H for this process.As a result of the modification of H (1) , also the B (2) coefficient in the Sudakov receives a µ (0),α R and µ (0),y R dependence through its dependence on H (1) within MiNNLO PS , see Eq. (12).For completeness we provide here also the standard µ R dependence of the Sudakov coefficients, whose complete scale dependence is implemented through To assess the theoretical uncertainty of our MiNNLO PS predictions, we can now vary µ R , µ (0),y R around their central scales, either simultaneously by a common factor or independently.Our default choice and the impact on Higgs production in bottom-quark annihilation will be discussed in detail in Section 4.
In Ref. [57], a scale Q = K Q M in the modified logarithm was introduced, dubbed resummation scale, which determines the region at large transverse momentum where resummation effects are smoothly turned off within MiNNLO PS .Since there is an interplay with scale of the Yukawa coupling µ (0),y R , we report the full scale dependence for of the hard-virtual coefficient function below.The resummation-scale dependence is derived by splitting the integral of the Sudakov form factor in Eq. ( 5) into one contribution from p T to Q and a second one from Q to M , where the second integral is then expanded in powers of α s (K R /K Q p T ) up to second order.While the logarithmic contributions are absorbed into a redefinition of the B(2) coefficient, see Eq. (4.8) of Ref. [57], the non-logarithmic terms are expanded outside the exponential factor and are absorbed into the hard-virtual coefficient function H.For consistency the scale of the strong coupling in the expansion of H is changed as follows: while the complete scale dependence of the expansion coefficients of H reads −n B πβ 0 A (1) log (µ We note that we refrained from reporting the factorization scale (µ F = K F p T ) dependence of the MiNNLO PS method here, which is absorbed into the collinear coefficient functions, see Eq. ( 7), as there is no direct interplay with µ (0),y

R
. Instead, we refer to Ref. [57] for the relevant formulae.

Alternative scale setting for the non-singular contribution (FOatQ)
In the following, we explore a different scale choice for the non-singular part R f in Eq. ( 15), which can improve the comparison with the fixed-order computations by reducing differences in the treatment of terms beyond accuracy.In particular, we rederive the MiNNLO PS formulae for the case where R f is evaluated at a hard scale Q h instead of the transverse momentum p T as used in the original formulation of the MiNNLO PS method. 3The singular part encoded through the D function remains to be evaluated at the scale p T .More precisely, we will set the scale of the PDFs to µ F = K F Q h and the scale of the extra powers of the strong coupling (those beyond the ones appearing at LO) to µ R = K R Q h , in order to facilitate the usual scale variations from the central scale Q h through K F and K R .
We start by writing the regular part as Here, the argument of R f being "(Q h )" simply denotes that µ R /K R = µ F /K F = Q h as opposed to "(p T )", which implies µ R /K R = µ F /K F = p T .As a result and to keep the singular terms evaluated at a scale p T the B function has to be modified as (2) We note that, as opposed to the implementation used in Ref. [52], we perform a factorization of the full Sudakov form factor, also in front of the regular part.By contrast, Ref. [52] employed a less accurate form of the Sudakov form factor, denoted as S in Eq. (2.30) of that paper, to be factored out in front of the regular part.The two approaches are equivalent up to the desired accuracy of the MiNNLO PS method.However, to us it seems more natural and practically simpler to use the same exponential factor for both the singular and the regular part.
This implementation of the scale choice can be turned on through setting the flag FOatQ 1 in the code.The two new ingredients that are required to implement Eq. ( 39) are D (1) (Q h ) and D (2) (Q h ), whose formulae we will derive and provide in the following.They are defined by writing the expansion of D with a scale setting of whereas the corresponding expression for µ To derive D(Q h ) we simply start from D(p T ) and consistently change the scale setting in the parton distribution functions (PDFs) and the strong coupling.To this end, we employ the first-order expansion of the DGLAP evolution of the PDFs and of the RGE evolution of the strong coupling Note that the function D is defined in a scale invariant way, so that D(Q h ) and D(p T ) is simply the same function, but with a different scale choice, and they differ only by terms beyond accuracy.When imposing the expansion of the DGLAP and the RGE evolution in Eq. ( 42) and ( 43), respectively, on the first-order coefficient D (1) to change its scale from p T to Q h , the logarithmic scale dependence is absorbed into the second-order coefficient D (2) .Thus, the expansion coefficients of D(Q h ) can be expressed as D (1) where the notation D (i) (p T ) µF=KFQ h simply denotes to take the exact functional form of D (i) (p T ), but setting the scale of the PDFs to K F Q h .Note that D(p T )| µF=KFQ h is not scale independent up to higher orders, while the function D(Q h ) is scale invariant due to the inclusion of the logarithmic scale-compensating terms.Moreover, we have introduced a symbolic shortcut in Eq. ( 45) for the effect of the convolution in the DGLAP evolution of the PDFs in D (1) whose explicit form reads P (0) ⊗ D (1)   symb Lastly, we only need to specify the relevant formulae to determine D (1) (p T ) and D (2) (p T ), for which we refer the reader to Eq. ( 27) and ( 28) of Ref. [45] and to Ref. [57] for the relevant K R and K F dependence.We stress that, apart from the formulae specified above, all scale dependent terms in K R and K F remain exactly as in the original formulation of MiNNLO PS .
Our notation here, in particular the simplicity of Eqs. ( 44) and ( 45), reflects exactly how the modifications for FOatQ have been implemented in the code.The implementation of the FOatQ option has been performed within the Powheg-Box-Res code, which will be included in a future release, so that this feature can be used for all publicly available MiNNLO PS codes.We have performed additional checks of our FOatQ implementation by comparing both prescriptions for Drell-Yan production and for Higgs-boson production in gluon fusion.

Results
We present numerical predictions for Higgs production in bottom-quark annihilation (b b → H) for the LHC at 13 TeV centre-of-mass energy.We choose the following input parameters: For the (stable) Higgs boson the mass is set to m H = 125 GeV and its width to Γ H = 0 GeV.The 5FS is used with massless bottom quarks, but with a non-vanishing Yukawa coupling y b , renormalized in the MS scheme.The Yukawa coupling is evaluated from an input value m b (m b ) = 4.18 GeV and evolved to its respective central scale µ (0),y R = m H using four-loop running, while scale variations are obtained from that central value with a three-loop evolution, consistent with the order of our calculation.This procedure follows the recommendation by the LHC Higgs cross section working group [73].We use the NNLO set of NNPDF40 [74] parton densities with 5 active flavours with α s (m Z ) = 0.118 via the LHAPDF interface [75] as our default setting unless specified otherwise.The central factorization and renormalization scales are set following to the MiNNLO PS method [44], and we also present results with the FOatQ 1 scale setting discussed in Section 3.3.The associated scale uncertainties are determined through the customary 7-point envelope obtained by varying independently the corresponding factors K R and K F by a factor 2, where the scale of y b is varied simultaneously with K R , i.e. µ (0),y Moreover, we choose K Q = 0.25, while we have checked that setting K Q = 0.5 leads to rather small differences in the results.These scale settings are in line with ones used for the analytic resummation for Higgs production in bottom-quark annihilation presented in Ref. [17], which we will compare to at the end of this section.
For the predictions matched to a parton shower we employ Pythia8 with the A14 tune (py8tune 21 in the input card).Since we are interested in inclusive Higgs production and comparisons to other theory calculations, the Higgs boson is kept stable, and the effects from hadronization, multi-parton interactions (MPI) and QED radiation are kept off.
We start the discussion of our phenomenological predictions by comparing predictions for the total inclusive cross section from MiNLO ′ and MiNNLO PS with fixed-order results at NLO and NNLO in Table 1.The fixed-order rates for b b → H production are obtained with the program SusHi [10,46] and they have been obtained with the input parameters specified above, while setting renormalization and factorization scales to µ R /K R = µ F /K F = m H .For MiNNLO PS we quote the results for both the standard and the FOatQ 1 scale setting, and find their total bb → H b b → H 0.646(0) +10.4% −10.9% pb 0.518(2) +7.2% −7.5% pb 0.571(1) +17.4% −22.7% pb 0.509(8) +2.9% −5.3% pb 0.508(4) +3.6% −4.3% pb Table 1: Total cross sections of Higgs-boson production in bottom-quark annihilation.The number in brackets denotes the numerical uncertainty on the last digit, while the uncertainty quoted is the theory uncertainty estimated via scale variations as described in the text.
cross sections to be very close.We further observe that the effect of NNLO QCD corrections, both when comparing the fixed-order numbers at NNLO QCD with NLO QCD from SusHi and when comparing MiNNLO PS to MiNLO ′ , is to reduce the cross section by more than 10%.Moreover, we observe a significant reduction of scale uncertainties once NNLO QCD corrections are included.Finally, our MiNNLO PS predictions are in agreement with the NNLO QCD cross section within the quoted uncertainties from missing higher-order corrections.We stress that due to the different scale setting and a different treatment of terms beyond accuracy MiNNLO PS predictions are not expected to coincide with NNLO QCD results, but rather to agree with them within the given uncertainties.
Although all the results in this paper are obtained with the correlated scale variation (K y R = K R ) for the Yukawa and strong couplings, we implemented in the code the possibility to change the two renormalisation factors in an independent way in order to study the Yukawa effects.If we vary K y R and K F keeping K R = 1, we observe a more symmetric scale variation for the total cross-section, precisely 0.509(8) +4.7% −5.1% pb with the standard settings without FOatQ 1.We continue by analyzing differential distributions.In Figure 3, we compare our MiNNLO PS predictions for the distribution in the rapidity of the Higgs boson (y H ) against the NNLO results [21]. 5To facilitate this comparison, we adapted the input settings from our default ones, which were outlined before, to the ones employed in the NNLO computation of Ref. [21].Specifically, we use the NNLO set of the CT14 PDFs [76] and we set the scales to µ R = K R m H and µ F = K F m H /4. We find a good agreement, both in terms of normalization and in terms of shape, between the two central predictions, with any difference being fully covered by the scale uncertainties.When looking at the scale-uncertainty bands, we notice some differences: the bands of MiNNLO PS result are more symmetric and they are slightly smaller than the NNLO ones, where the central prediction lies at the upper edge of the prediction for central rapidities.Apart from that, we have checked that impact of the parton shower on the y H distribution is very moderate.We note that for the results presented in the remainder of this paper we will employ again our default input settings outlined at the beginning of this section.
In Figure 4, we consider distributions that require the presence of at least one jet and compare MiNNLO PS (blue, solid curve) with MiNLO ′ (black, dotted curve) predictions.The left plots show the rapidity distribution of the leading jet (y j 1 ) and the right plots show the rapidity difference between the Higgs boson and the leading jet (∆y H,j 1 ).In the upper (lower) plots a transversemomentum cut of p T,j 1 > 30 GeV (p T,j 1 > 120 GeV) is imposed on the leading jet.For these distributions MiNNLO PS and MiNLO ′ predictions are both formally NLO QCD accurate.The purpose of this comparison is to validate that the NNLO corrections added by MiNNLO PS do not alter significantly the MiNLO ′ result.Indeed, we observe that the MiNLO ′ and MiNNLO PS results have very similar shapes and that they are fully consistent within the quoted scale uncertainties.Moreover, the harder the required jet is, i.e. by increasing p T,j 1 , the more similar MiNLO ′ and MiNNLO PS predictions become.The distribution in the rapidity difference between the Higgs and the leading jet is, as expected, peaked around zero for p T,j 1 > 30 GeV.This is due to the fact that the jet and the Higgs boson are approximately balanced in transverse momentum, which typically leads to the rapidity difference being centered around zero.By contrast, when the leading jet is boosted (i.e. for p T,j 1 > 120 GeV), the Higgs boson and jet tend to be slightly farther apart in rapidity and we observe a dip in the ∆y H,j 1 distributions at central rapitidies.
Next, we consider the transverse-momentum spectrum of the Higgs boson (p T,H ), focusing on the large-p T,H region.Here we validate our MiNNLO PS generator, which formally is again only NLO accurate in QCD for p T,H ≳ m H , against appropriate fixed-order calculations.In Figure 5 (left) we compare our Powheg implementation for H+jet production (dark-blue, double-dash-dotted curve) with a fixed-order calculation for H+jet production obtained from Ref. [12] (brown, long-dashed curve) requiring p T,j > 10 GeV, while Figure 5 (right) shows the analytic p T,H spectrum up to α 2 s from Ref. [17] (red, dashed) and our MiNNLO PS prediction with (magenta, dash-dotted curve) and without the FOatQ 1 scale setting (blue, solid curve).As expected, we find full agreement between our predictions and these fixed-order predictions in the regime of large p T,H .By contrast, the fixed-order calculations yield a divergent cross section for p T,H → 0, while the MiNNLO PS prediction is finite in this region.We will study the low p T,H region in more detail at the end of this section.
We now turn to the comparison of MiNLO ′ with MiNNLO PS predictions at the fully differential  Figure 4: The left plots show the rapidity distribution of the leading jet (y j 1 ) and the right plots show the rapidity difference between the Higgs boson and the leading jet (∆y H,j 1 ).In the upper (lower) plots a transverse-momentum cut of p T,j 1 > 30 GeV (p T,j 1 > 120 GeV) is imposed on the leading jet.level in Figure 6, which allows us to assess this size of the MiNNLO PS corrections on top of MiNLO ′ .Specifically, the transverse momentum distributions of the leading jet (p T,j 1 ), subleading jet (p T,j 2 ) and the Higgs boson (p T,H ) as well as the rapidity distribution of the Higgs boson (y H ) are shown.We observe that the MiNNLO PS corrections have a significant impact at small transverse momenta, dampening the spectrum significantly more in this region.Moreover, scale uncertainties are substantially reduced in that region.On the contrary, at large transverse momenta the two predictions essentially coincide with each other, given that they are of the same formal Figure 5: The left plot shows our Powheg implementation for H+jet production (dark-blue, double-dash-dotted curve) and a fixed-order calculation for H+jet production obtained from Ref. [12] (brown, long-dashed curve) requiring p T,j > 10 GeV.The right plot shows the analytic p T,H spectrum up to α 2 s from Ref. [17] (red, dashed) and our MiNNLO PS prediction with (magenta, dash-dotted curve) and without the FOatQ 1 scale setting (blue, solid curve).accuracy in that region.We stress that these results are obtained using the correlated scale variation (K y R = K R ).However, we also isolated the effects of the Yukawa scale variation fixing K R = 1 and varying just K y R .In the low p T,H region, we obtained more symmetric uncertainties reflecting the previous observation on the total cross-section, while we notice that in the tail of the transverse momentum spectrum the scale variation is symmetric for both the settings, with a bigger uncertainty band for the correlated case K y R = K R .Finally, looking at the y H distribution, we observe a rather constant and flat negative correction of about −12% due to MiNNLO PS and a significant reduction of the scale uncertainties.
We end this section by comparing the predictions for the p T,H spectrum from our MiNNLO PS generator against the analytic resummation at NNLO+NNLL of Ref. [17] (green, double-dashdotted curve) in Figure 7.The left plots show this comparison zoomed into the low p T,H region, while the right ones include the p T,H spectrum up to 270 GeV.The upper plots include MiNNLO PS results at the LHE level and the lower ones after showering with Pythia8.We find that again MiNNLO PS results with and without the FOatQ 1 setting are very close.By contrast, the NNLO+NNLL prediction is softer than the MiNNLO PS ones in the low p T,H region, where peak of the spectrum is shifted towards higher values for MiNNLO PS .At high p T,H , where the predictions have the same accuracy, they coincide again, when considering the LHE-level MiNNLO PS results.Including the parton showering effects from Pythia8 tends to slightly worsen the agreement at small p T,H , moving MiNNLO PS away from the NNLO+NNLL scale band.Considering the fact that the MiNNLO PS predictions feature smaller scale uncertainties at small p T,H than the more accurate NNLO+NNLL prediction, it is clear that the MiNNLO PS scale band does not reflect the actual size of uncertainties at small p T,H , which would require additional variations within the shower settings.Since the assessment of the Pythia8 uncertainties is an issue intrinsic to the shower itself and not related to the NNLO+PS matching, we refrain from further studying these aspects here.We note, however, that in the small p T,H region a massless approximation misses potentially relevant mass effects, such that a combination with a massive 4FS calculation might be needed, which is beyond the scope of this paper and left for future work.The shower also induces some effect at large p T,H moving the MiNNLO PS prediction up by about 10%, which is however well within the given scale-uncertainty bands.Figure 7: Predictions for the p T,H spectrum from our MiNNLO PS generator compared to the analytic resummation at NNLO+NNLL of Ref. [17] (green, double-dash-dotted curve).

Summary
We have presented a first NNLO+PS calculation of the b b → H production process at the LHC.
To this end, we have revised the MiNNLO PS method to account for an overall Yukawa coupling that is renormalized in the MS scheme.Moreover, we have implemented an alternative approach to set the renormalization and factorization scales of the NLO colour-singlet plus jet calculation within MiNNLO PS .
We have performed an extensive validation of our MiNNLO PS predictions, by comparing against inclusive NNLO QCD fixed-order results for the total cross section as well as the rapidity distribution, against lower-order MiNLO ′ results for jet-related quantities, and by comparing to NLO QCD accurate calculations at large values of the Higgs transverse momentum.In all cases, we found consistency of our MiNNLO PS results in the relevant kinematical regions.At last, we considered the Higgs transverse-momentum spectrum and showed a comparison of MiNNLO PS against an analytically resummed calculation at NNLO+NNLL.Uncertainties of the MiNNLO PS predictions seem to be underestimated at small-p T,H , given that they are smaller than the ones of the more accurate NNLO+NNLL prediction and that the two results barely agree within the quoted uncertainties.The quoted uncertainties for MiNNLO PS predictions however do not include uncertainties related to the parton shower.The latter requires a dedicated study which we leave to the interested reader.
We reckon that this MiNNLO PS generator could be useful not only in the context of direct searches for b bH final states, in particular in some BSM context, but also for HH measurements, where the b bH process constitutes a major irreducible background.Moreover, the present 5FS calculation of this process shall be considered as a first step to a full NNLO+PS description of b bH production, with the prospect of a 4FS b bH calculation with massive bottom quarks and eventually a full 4FS-5FS combination at NNLO+PS accuracy, which we leave for future work.Such combination could be particularly relevant in the small p T,H region where mass effects of the bottom quarks are potentially relevant.

Figure 1 :
Figure 1: Sample Feynman diagrams for Higgs production in association with bottom quarks.

Figure 2 :
Figure 2: LO diagrams for HJ production in 5FS.

Figure 3 :
Figure 3: Comparison of MiNNLO PS predictions (blue, solid) with the NNLO results of Ref. [21] (red, dashed) for the rapidity distribution of the Higgs boson.

Figure 6 :
Figure 6: Comparison of MiNLO ′ (black, dotted) with MiNNLO PS predictions (blue, solid) for the transverse momentum of the leading jet (p T,j 1 ), the subleading jet (p T,j 2 ), and the Higgs boson (p T,H ) as well as the rapidity distribution of the Higgs boson (y H ).