The fully differential hadronic production of a Higgs boson via bottom quark fusion at NNLO

The fully differential computation of the hadronic production cross section of a Higgs boson via bottom quarks is presented at NNLO in QCD. Several differential distributions with their corresponding scale uncertainties are presented for the 8 TeV LHC. This is the first application of the method of non-linear mappings for NNLO differential calculations at hadron colliders.


Introduction
The Large Hadron Collider is now at its third year of successful operation and both ATLAS and CMS report tantalizing hints of a Higgs boson at about 125 GeV. By the end of the 2012 run the experiments are likely to be able to either confirm those hints as a firm discovery or else exclude any Standard Model (SM) Higgs boson. In the event of a firm discovery further detailed examination of various production and decay channels will be necessary to determine the nature of the Higgs sector.
The dominant production channel in the SM, but also in all non-fermiophobic models of new physics, is single Higgs hadroproduction. Within the SM the production mechanism is dominated by gluon fusion, since the alternative mechanism of quark annihilation is severely suppressed by the small Yukawa coupling of bottom and light quarks to the Higgs boson. However, if the Higgs sector is non-minimal, as is the case in any two-Higgs-doublet model (among which the MSSM is the most studied example) the Yukawa coupling to down-type quarks is enhanced by a factor of tan β (the ratio of the vacuum expectation values of the two doublets) and the contribution of the bb → H process becomes significant. Furthermore the production cross section through gluon fusion decreases due to the enhanced, negative top-bottom interference diagrams. In such a scenario, the production of a Higgs boson via bb pairs contributes much more than in the SM, and a detailed description of this process is desirable. In other BSM models, for example in models with dynamically generated Yukawa couplings [1,2], both the bottom and the charm quarks have enhanced couplings to the Higgs boson and charm annihilation becomes important as well.
The experimental searches are currently focused on measuring an enhanced production rate via bottom annihilation in the τ + τ − decay channel with the MSSM as the default BSM model [3,4]. There are, moreover, several studies on measuring single Higgs decaying to bottom quarks in more generic models in which bottom annihilation is the dominant production channel, using, for example, three b-tagged jets [5,6], or measuring the ratio of three heavy (c-or b-) jet events to three b-jet events to discriminate between classes of models with two Higgs doublets [7].
Bottom quark annihilation has been the subject of much theoretical discussion in the last decade due to the freedom in treating the initial state bottom quarks. Bottom quarks lie in an intermediate mass range between the non-perturbative regime of the proton mass and the typical scale of a hard scattering event at the LHC. One can retain their small mass in the calculation, and exclude them from the proton constituents (four flavor scheme -4FS) or treat them as massless partons with their own parton distribution functions (five flavor scheme -5FS). In the 4FS the inclusive cross section develops large logarithms ∼ log( m b m H ) due to the collinear production of b-quarks which is regulated by the bottom mass. In the 5FS these logarithms are re-summed to all orders by the DGLAP evolution inside the bottom PDFs, for all scales up to the factorization scale adopted in the calculation. Improved convergence of the perturbative expansion is an advantage of the 5FS approach, but at the same time it makes the 5FS prediction very sensitive to the choice of factorization scale. It has been realized that if the factorization scale is set to low values ∼ m H /4, both the 5FS and the 4FS predictions for the inclusive cross sections agree with each other within their respective uncertainties [8][9][10], and there is an open discussion as to how one would combine information from both approaches [11,12].
In the 4FS, the lowest order process would be gg → bbH which begins at order α 2 s in the QCD perturbative expansion and is known at next-to-leading-order (NLO) in QCD [13][14][15][16]. The process bg → bH, which starts at order α s , has also been studied at NLO in QCD [17] and with electroweak (EW) corrections [18]. In the 5FS the lowest process is bb → H. Hence the LO 4FS process where a non-collinear bottom pair is observable, is only reached for the first time at NNLO in the 5FS. The inclusive cross section, in the 5FS, of bb → H is known at NNLO in QCD [19] as well as at NLO in EW [20]. NNNLO threshold re-summed soft and collinear terms are also known [21] and the transverse momentum distribution of the Higgs boson has been studied with re-summation techniques [22,23] . Also known at NNLO are the zero-, one-and two-jet rates and related distributions [24], quantities which can be obtained already from the differential H + jet computation at NLO [25] in combination with the fully inclusive NNLO cross section.
In this paper we present the fully differential NNLO QCD cross section for bb → H in the 5FS within the SM. NLO computations are currently performed with very well automated methods. Obtaining fully differential cross sections and decay rates at one order higher in the perturbative expansion requires the solution of new challenging problems. Regarding the treatment of the real emissions, pioneered for NLO computations in [26,27], rapid progress has been made in the last decade , mainly focusing on the treatment of the double real emission 1 , which resulted in the fully differential calculations of Higgs production via gluon fusion [54][55][56][57], Drell-Yan [58][59][60][61][62], associated Higgs production with a vector boson [63], three jet production from e + e − [64][65][66][67] and diphoton production [68].
However, further development of methods and new ideas are necessary for efficient cancellations of infrared singularities and evaluations of novel two-loop amplitudes in more complicated LHC processes.
With this paper we also take the opportunity to complete the second NNLO application, after the fully differential decay H → bb [69], using the method of nonlinear mappings to factorize singularities in the double real corrections [70]. The double real contributions have often been regarded as the bottleneck of NNLO, and this paper therefore also demonstrates the validity of the approach as a method for NNLO corrections in hadronic collisions.
The paper is organized as follows: in section 2 we set up the notation and describe the main components of the calculation. In section 3 we provide some detail about the treatment of the separation of soft and hard contributions. In section 4 we describe the treatment of the double real and the real-virtual components. In section 5 we present the way we perform the (non-trivial at NNLO) convolutions for the collinear subtraction terms in mass factorization. In section 6 we provide various numerical results both on jet rates, p T and rapidity distributions; demonstrate the completely differential nature of our calculation and provide typical results for the case in which the Higgs boson decays to two photons, including standard experimental cuts on photon momenta and isolation.

Fully differential calculations
One of the merits of fully differential calculations is the possibility to arrive at theoretical predictions for observables in the presence of final state phase-space cuts, like those used in experimental analyses, under the precondition that the observable defined is infra-red safe. Throughout this article the dependence on such arbitrary phase-space constraints will be contained in the jet-function J ({p} f ), where {p} f denotes the set of final state momenta in the lab frame. We will refer to the fully differential cross section as σ[J ], which we schematically define as where the sum is over all final states f . The usual role of the jet function J is to apply arbitrary final state phase-space cuts while ensuring infra-red safety. Here we promote it to a further task, which is to keep track of the bin-integrated cross section for any given differential observable with or without applying phase-space cuts. This can be achieved simply at the level of Monte Carlo integration by passing to J not only the set of final state momenta but also the weight of the given event. The role of the jet function becomes crucial in all amplitudes that have soft and collinear singularities which are regulated by counter terms. In such cases the jet function is keeping track of the kinematics of every subtraction term.

Hadronic cross section
We consider the following hadronic process where P 1 , P 2 are the incoming hadrons, H denotes the Higgs boson and X generically denotes surplus QCD radiation in the final state. The Higgs boson is assumed to couple only to bottom quarks via the SM Yukawa interaction. Assuming the usual factorization, the fully exclusive hadronic cross section can be written as where the f i (x) denote the bare (unrenormalized) parton distribution functions (PDFs) in the 5FS, x 1 and x 2 are the usual Bjorken-x momentum fractions of the partons i 1 and i 2 respectively, and τ = m 2 H S , where m 2 H is the (on-shell) mass of the Higgs-boson and S is the square of the total center of mass (CoM) energy of the colliding hadrons. By σ i 1 i 2 →HX we denote the partonic cross section for the processes The PDFs we have inserted in eq.(2.3) are bare and we still have to rewrite them in terms of the renormalized PDFs. This step will introduce collinear counter terms that cancel the initial state collinear singularities of the partonic cross section, which remain after all real and virtual corrections are added together. This cancellation is achieved fully numerically in our calculation. We outline the way collinear counter terms can be computed processindependently in section 5.

Partonic cross sections
Expanding the partonic cross section to NNLO in QCD we obtain where a sum is implied over all final state flavors k and l leading to a possible subprocess.
Here y b = y b (µ) and α s = α s (µ) are the MS renormalized bottom Yukawa and strong couplings, with 5 active flavors. We set the dimensional regularization scale µ to be equal to both the renormalization and factorization scales, µ R and µ F . Separation of the two scales can be easily achieved via the relations given in appendix B. The leading order contribution, fig. 1, is denoted by σ B ij→H . At NLO there are two separate contributions (see fig. 2): Figure 2. Some of the diagrams contributing to bb → H at NLO. These contributions are denoted as real (left) and virtual (right). Figure 3.
Some of the diagrams contributing to bb → H at NNLO. These contributions are denoted as double real (left), real-virtual (center) and double virtual (right).
corrections are absorbed into the renormalized couplings, see e.g. [19,69] or any textbook on QCD. In contrast, infra-red divergences cancel only after summing all real and virtual corrections contributing to a given infra-red safe observable.
Factorized singularities on the unit hypercube may be dealt with using the plusdistribution expansion where D m (x) = ln m (x) x + and the plus-distribution is defined through Beyond NLO a factorization of singularities is highly non-trivial. In this work it is achieved systematically using the method of nonlinear mappings [70]. Care must be taken when dealing with infrared singularities of real emission amplitudes: the plus-distribution, eq.(2.6), also acts on the jet function, such that cancellations happen also at the differential level and are therefore fully local.
We now give a brief overview of the matrix elements and phase-space measures required for the computation of the partonic cross section. Here we assume the amplitudes to be color and spin summed. Averaging and phase-space symmetry factors will be explicitly factored out.

i) Purely virtual corrections:
The purely virtual corrections include the Born, virtual and double virtual contribu-tions and are of the form where s 12 = (p 1 +p 2 ) 2 is the partonic CoM energy and N bb = 1/36. The corresponding phase-space volume element is trivially given by dΦ 2→1 = 2πδ(s 12 −m 2 H ), constraining Regarding the computation of amplitudes we refer the reader to [69] where the full virtual matrix-elements can be found. Explicit expressions for these contributions are given in section 4.4.
ii) Single real emissions: The single real emissions include real and real-virtual corrections and are of the form Further details are given in section 4.1.
iii) Double real emissions: These are of the form Further details are given in section 4.3.

Separation of soft and hard
Since in the soft limit of a 2 → 1 process the produced particle is at rest in the partonic center of mass frame, there is no difference between the soft piece of a fully differential partonic cross section and that of a fully inclusive partonic cross section. It is therefore very convenient to isolate the soft contribution (σ S ) to the partonic cross section (σ) from the hard one (σ H ), i.e.
This allows for a fully analytic treatment of σ S , while σ H must, as far as external kinematics are concerned, be treated numerically. Let us introduce the variable Then the soft limit of all real emission amplitudes corresponds to z → 1, which identifies the production threshold. Given that infrared singularities are of logarithmic nature, the divergence at z = 1 can be exposed as follows where σ V denotes the purely virtual correction, while σ (n) R (z, ) denotes real corrections collectively (at NNLO this includes both real-virtual as well as double real corrections).
Separation into soft and hard parts can now be achieved by adding and subtracting the soft limit from the second term in the above, yielding Of course this decomposition of the partonic cross section into its soft and hard components is not unique: one could use any other subtraction term with the correct limit, thereby including, for example, the luminosity function. Our choice, however, has the nice property that the soft part σ S can be expanded purely in terms of δ-and plus-distributions via eq.(2.6), Thereby all threshold divergences between σ V and σ (n) R are canceled analytically, leaving only a finite threshold contribution. Furthermore this framework provides a natural way to incorporate threshold re-summation in fully differential calculations.

The single real
The single real partonic cross section may be expressed as We define d Φ 2 to be the conventional phase-space volume dΦ 2 up to some renormalisation constants. Here we have to consider 6 separate channels: The corresponding amplitudes may all be found in [69]. A convenient phase space parametrization is given by where λ ∈ [0, 1], with the Lorentz invariants taking the simple form Note that the singularities of s 13 and s 23 are factorized in λ, (1 − λ) and (1 − z) which allows for a simple subtraction of the poles using eq.(2.6). This also allows us to identify The calculation of the hard part then trivially follows from eq.(3.4).

The real-virtual
The real-virtual partonic cross section may be expressed as where we have taken the liberty to define d Φ 2 to equal eq.(4.2) up to some renormalization constants. Then The real-virtual amplitude can be obtained from the corresponding one from the decay process H → bb published in [69] by crossing particles to the initial state. The box integrals we encounter in this amplitude are entirely expressible in terms of Gauss' hypergeometric function 2 F 1 (1, − , 1 − , z) where z can be in any of the three sets S f ine , S inv and S nl : When attempting a direct subtraction of the singularities created by the real emission, the points of subtraction overlap with singular points of the hypergeometric functions in the box integrals. It was found in [69] that one can apply transformations on the argument of the functions to circumvent this difficulty. Since here we are no longer in the euclidean regime of this amplitude, the required transformations are different than in [69]. Analyzing integral representations, we find that we have to apply the following identities: • If z ∈ S f ine the soft-collinear limits are well defined.
• If z ∈ S nl we apply • If z ∈ S inv we employ the argument inversion, After these transformations are applied, the singularities corresponding to the real emission are factorized in λ, (1 − λ) and (1 − z). The soft singularity structure of the real-virtual may then be extracted as In the soft limit only the m = 2, 4 coefficients survive and the integration over λ can be done analytically. The explicit expressions for the soft limit can be found in appendix A. The computation of the hard part then follows from eq.(3.4). While the structure is more complicated than in the case of the single real, a direct subtraction via eq.(2.6) can still be achieved in a straightforward manner. In order to obtain the final Laurent expansion in we employ the library HypExp [71] to expand the hypergeometric functions in terms of polylogarithms.

The double real
The double real partonic cross section can be written as where dΦ 3 is equal to the conventional three-particle phase space element dΦ 3 up to renormalization constants. Using the discrete symmetries of the squared amplitudes we are able to considerably reduce the number of independent channels, which one has to implement separately. These symmetries are due to the charge invariance of all the bb → H double real amplitudes (exchanging q ↔q or b ↔b leaves the amplitudes invariant). This leaves us with the following list of channels (4.8) By crossing partons from the initial to the final state, we can obtain all of the above from the three amplitudes |M bbggH→0 | 2 , |M bbbbH→0 | 2 and |M bbqqH→0 | 2 published in [69]. In order to deal with the intricate singularities, their factorization and subtraction, we refer the reader to the methods developed in [70], which we have implemented faithfully. As in the single real emissions, the double soft singularity occurs at the threshold. Its structure may be identified as (4.9)

The soft
Let us expand σ S in the strong coupling where The NLO correction ∆ S,NLO may be expressed as while the NNLO correction ∆ S,NNLO can be expanded as follows   The explicit soft limits of the real-virtual and double real pieces that are included in σ S can be found separately, and with their full color-dependence, in appendix A.

Collinear factorization
Parton distribution functions are renormalized to absorb initial state collinear singularities viaf where µ is the factorization scale and f i are the bare parton densities. In the following discussion summation over indices will always be assumed unless explicitly stated. We will also need the convolution integral, which is defined as The kernel Γ ij is defined in the MS scheme by where the coefficients of the expansion in the strong coupling involve the Altarelli-Parisi splitting functions P n ij . Specifically, with β 0 = 11 4 − 1 6 N F . Let us define the inverse of the kernel Γ ij as such that it satisfies the condition (Γ ik ⊗ ∆ kj ) (z) = δ ij δ(1 − z). Solving for the coefficients yields The strong coupling expansion of the bare PDFs then reads In evaluating the collinear counter terms we encounter convolutions of the type (f ⊗ ∆)(x), where the function f is regular and ∆(x) can in general be written as Care has to be taken with convolutions over D n . Since the integration does not start at zero, a boundary term must be included Because of the downward sloping shape of all parton distribution functions, a quadratic remapping of the integration variable y was found to optimize the convergence behavior, i.e. we parametrized the integral like with z uniformly distributed between 0 and 1. In our code, this integration is carried out numerically. The integration is onedimensional, which makes a simple deterministic trapezium integration with about 50.000 points the simplest option. The result of the integration is accurate to at least 5 digits, which is usually below the precision of the Monte Carlo integration. The precision of the integration can be arbitrarily increased by increasing the number of points used. For every bare PDF used, we construct a one-dimensional grid in the Bjorken-x variable and interpolate from it during runtime.
An alternative to constructing a grid is to perform the integration numerically along with the phase space ones, thereby increasing the dimensionality of the Monte Carlo integration by one (or by two in the case of double NLO kernels convoluted with the Born). We have implemented this as well and found that it yields the same results as the grid approach.
This procedure allows us to expand the (singular) bare PDFs via eq.(5.10) order by order in the dimensional regulator and substitute them directly in eq.(2.3). The singularities in the resulting convolutions, appearing as poles in the -expansion, cancel the initial state collinear singularities of the partonic cross section. This cancellation is achieved numerically in our calculation and can be observed bin by bin in e.g. the rapidity distribution of the Higgs boson. One can achieve this cancellation in each initial state channel separately, at the cost of separating the convolution integrals depending on the initial state parton in the convolution, i.e. by not performing the implicit j-summation in eq.(5.11).
It is worth pointing out that the procedure described here is entirely generic, i.e. it provides the collinear counter terms for any NNLO process numerically. Moreover, we thereby circumvent the usual insertion of eq.(5.1) in the equivalent of eq.(2.3) for renormalized quantities and the resulting cumbersome and process specific analytic treatment of the convolutions.

Numerical results
We have performed a number of tests to ensure that our results are consistent with each other and with results available in the literature: • We have implemented the entire calculation in two different computer codes, one in Fortran and one in C++, and all results agree within their respective Monte Carlo errors, both inclusively and differentially.
• The coefficients of all poles in the -expansion of all cross sections cancel both inclusively and differentially for the entire process and also for all individual initial state channels.
• The inclusive cross section agrees with the one available in [19] and from ihixs [73] and so does the inclusive cross section per initial state channel. This is the first independent check of the inclusive cross section published in [19] and adopted in [73].
• The soft limit of both real-virtual and double real contributions were computed both numerically (as a limiting case of the generic matrix elements) and analytically. Moreover the integrated double real contributions were found to agree with an analytic computation provided by [72].
• The subtraction process for every double real integral was implemented in two different ways and were found in complete agreement. We present results for the LHC with a center of mass energy of 8 TeV. We fix the mass of the Higgs boson at 125 GeV. We have used the MSTW2008 (68%CL) PDFs for all results presented here. The value of α s at m Z that we use is the best-fit value of the PDF set at the corresponding order. We use µ R = m H as the central renormalization scale. The value of α s used is run from m Z to µ R through NNLO in QCD. The mass of the bottom quarks is set to zero in all matrix elements, consistently with the 5FS choice. The bottom Yukawa coupling, however, depends on the mass of the bottom. The Yukawa coupling at µ R is obtained from the Yukawa coupling at µ * = 10 GeV, using m b (µ * ) = 3.63 GeV. We do not vary µ R in what follows, since the µ R scale dependence of the total cross section has been found to be very mild. We have also checked that the µ R -dependence of differential distributions is very small.
Previous studies have shown that the inclusive cross section is very sensitive to the choice of factorization scale. Arguments related to the validity of the 5FS approximation with respect to the collinearity of final state b-quarks, as well as to the matching to the 4FS calculation or to the need for a smoother perturbative expansion, point to factorization scales that are much lower than the Higgs boson mass. We adopt the choice µ F = m H 4 as a central scale and vary it in the range [ m H 8 , m H 2 ] to estimate the related uncertainty. All Monte Carlo integrations was performed with the Cuba [74] implementation of the Vegas algorithm.
The rapidity distribution of the Higgs boson is shown at fig. 4. As expected, the perturbative expansion is converging smoothly for this choice of central µ F and the NNLO uncertainty band is entirely engulfed by the NLO one.
The transverse momentum distribution for the Higgs boson is shown in fig. 5. This observable starts at NLO in QCD in the 5FS, and the fixed order prediction fails, as usual, to describe the very low p T spectrum due to the related large logarithms. At the large p T range we see that the NNLO calculation leads to a harder spectrum than the NLO one and the NLO scale uncertainty fails to capture this feature. This implies that great care should be taken when relying on NLO predictions for observables that are highly exclusive in the transverse momentum of the Higgs boson. The differential distribution in both the rapidity and the p T of the Higgs is shown in fig. 6, both in a three-dimensional lego plot and in a density plot. We see that the bulk of the events are produced centrally (with |y| < 2.5) and at relatively low p T ( 35 − 50GeV). In fig. 7 we show the cumulative distribution of the Higgs transverse momentum. This observable is equivalent to the cross section in the presence of a jet veto at NLO, but only related indirectly at NNLO. In fig. 8 we present the cross section in the presence of a jet veto. We see again that the perturbative description for high p T cut-offs is satisfactory (despite the discrepancy in high p T between NLO and NNLO, which is, in absolute terms, unimportant), while for cut-offs lower than 20 GeV the NLO description does not coincide with the NNLO one. The vanishing of the uncertainty around 15 GeV (which in the case of the jet veto is taking place at a slightly lower p T -veto value) is a feature reminiscent of a similar situation in Higgs production via gluon fusion [75]. The fixed order prediction in this region is very stable under varying the factorization scale, and any residual uncertainty in quantities like the acceptance in the presence of a veto is driven by the uncertainty in the total cross section. Various approaches to assign a larger uncertainty to similar observables involving re-summation exist, see for example [76]. An important observable in bb → H is the cross section for zero, one and two jets. We use the anti-k T algorithm [77] for jet clustering 2 with a cone in the y − φ plane of radius R = 0.4. We show in fig. 10 the jet rates as a function of the jet p max T used to define them. Here we do not distinguish between b-jets and light jets. We find the jet rates for p max T = 20GeV to be in agreement with those published in [24]. A wealth of information can be derived from examining the contribution of the different initial state channels to differential distributions. The six initial state channels that contribute to our NNLO calculation have singularities in various collinear regions that are canceled against the collinear counter terms from mass factorization. In order to make the cross section per channel finite one has to use collinear counter terms that include Γ the collinear counter terms numerically this modification was relatively easy to achieve. Initial state channel contributions to differential distributions have a strong dependence on the factorization scale, as do initial state channel contributions to the inclusive cross section. In fig. 11 we see the contributions to the Higgs boson p T distribution from each channel, for various factorization scales ranging from m H /16 to 2m H .
Within the 5FS, the factorization scale regularizes the collinear singularities which in the 4FS are regularized by the bottom mass. At NNLO, three initial state channels, bb, bg and gg share common collinear configurations whose leading logarithms cancel each other in different bins of the Higgs p T distribution. In the zero p T bin, in particular, squared log-  fig. 9. These features are also seen in the rapidity distribution of the Higgs boson per initial state channel, shown in fig. 12 for various values of µ F . There it is clearly seen that a scale like µ F = m H /4 eliminates the cross-channel cancelations but a lower scale µ F = m H /16 leads to a reduced, bg-dominated prediction.
We turn now to more exclusive observables. In large tan β models where the Higgs boson production gets significant contribution from the bottom quark annihilation process, one would like to examine differential distributions involving decay products of the Higgs boson, with cuts necessary in the experimental analyses. We focus here, for demonstration purposes, on the case where the Higgs boson decays to two photons. In such an analysis the minimal cuts used by CMS and ATLAS include: • A cut on the p T of the leading photon: p T ;1 > 40GeV.
• A cut on the p T of the trailing photon: p T ;2 > 25GeV.
• An isolation cut on photons: no jet is allowed in a cone of radius 0.4 around any of the two photons if it is p T > 15GeV.
We treat the Higgs boson in the zero width approximation in this article. We defer a more realistic treatment of the Higgs propagator to future work. Within this setup we show in fig. 13 the distribution of the average transverse momentum of the two photons and the distribution of the absolute of the difference in pseudorapidity between the two photons, Y * = 1 2 |y 1 − y 2 |.

Conclusions
We have presented the fully differential NNLO calculation of bb → H, a process of prime phenomenological importance for the LHC in all models with enhanced bottom Yukawa couplings. This is the first independent cross-check of the inclusive NNLO calculation performed in [19]. We have presented a variety of differential distributions for Higgs production that can only be obtained with a fully differential calculation and are useful for assessing the quality of the perturbative expansion and the level under which several features are under control at a fully differential level. We have also presented predictions for fully exclusive observables for the bb → H → γγ process in the presence of tight cuts on the final state photons including isolation cuts, demonstrating that our calculation can fully simulate any experimental setup at the partonic level. This is the second application of our approach to treat real emission singular amplitudes at NNLO [70]. It is the first application for the more complicated case of a hadron collider process. We find the approach particularly beneficial, both in terms of automatization and in terms of performance of the resulting numerical code. We find that the improvement in performance compared to the sector decomposition approach is significant. We intend to release the computer code in the near future and we defer for then any detailed comments on performance issues.
A study of significantly wider scope, including the production via gluon fusion in models with enhanced bottom Yukawa couplings, as well as the decay of Higgs to bottom quarks or tau leptons would vastly benefit the experimental searches. We defer such a study for a future publication.

A.2 Real-virtual
We decompose the real-virtual soft contribution as where only ∆ (2) RV and ∆ (4) RV are non vanishing and given by ∆ (2)

B Scale separation
The renormalization and factorization scales, µ R and µ F , can be conveniently separated by first setting µ = µ F and then applying the following relations