Associated production of a Higgs boson decaying into bottom quarks and a weak vector boson decaying leptonically at NNLO in QCD

We present the calculation of next-to-next-to-leading order (NNLO) corrections in perturbative QCD for the production of a Higgs boson decaying into a pair of bottom quarks in association with a leptonically decaying weak vector boson: pp→VH+X→ℓℓ¯bb¯+X.\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathrm{pp}\to V\mathrm{H}+X\to \mathrm{\ell}\overline{\mathrm{\ell}}\mathrm{b}\overline{\mathrm{b}}+X. $$\end{document} We consider the corrections to both the production and decay sub-processes, retaining a fully differential description of the final state including off-shell propagators of the Higgs and vector boson. The calculation is carried out using the antenna subtraction formalism and is implemented in the NNLOjet framework. Clustering and identification of b-jets is performed with the flavour-kt algorithm and results for fiducial cross sections and distributions are presented for the LHC at s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s} $$\end{document} = 13 TeV. We assess the residual theory uncertainty by varying the production and decay scales independently and provide scale uncertainty bands in our results, yielding percent-level accurate predictions for observables in this Higgs production mode computed at NNLO. Confronting a na¨ıve perturbative expansion of the cross section against the customary re-scaling procedure to a fixed branching ratio reveals that starting from NNLO, the latter could be inadequate in estimating missing higher-order effects through scale variations.


Introduction
One of the highest priorities of the LHC physics programme is the detailed exploration of the mechanism of electroweak symmetry breaking that predicts the existence of the Higgs boson and its interactions with the fermions and gauge bosons of the Standard Model (SM). In July 2012, the ATLAS and CMS Collaborations at the LHC reported the discovery of a resonance with a mass close to 125 GeV [1,2]. At the current level of accuracy, the discovered particle proves to be consistent with the Higgs boson predicted by the SM but the limited precision of some of the measurements still leaves room for possible alternative interpretations beyond the SM. Measurements of various properties of the Higgs boson have been carried out since then. One of the main goals of the completed Run II at √ s = 13 TeV and the future Run III at √ s = 14 TeV of the LHC is to test the coupling strength of the discovered Higgs-like particle to known SM particles through the study of a variety of processes at these increased luminosity and collisions energies.

JHEP10(2019)002
The production of a Higgs boson (H) in association with either a W ± or a Z boson and possible hadronic jets -also known as Higgs Strahlung -is among the most promising class of channels that can lead to the accurate determination of the Higgs-boson couplings. These were also the channels that were mainly probed during the search for a light Higgs boson at the Tevatron; and the observation of excess events at the Tevatron turned out to be consistent with the observed Higgs boson at the LHC [3].
At LHC energies the V H processes are the third (V = W ± ) and fourth (V = Z) largest production channels after the dominant gluon-gluon and vector-boson-fusion ones. These classes of Higgs production modes provide the opportunity to probe the gauge-boson-Higgs vertex (V V H) separately for V = W ± and V = Z. Moreover, a second and particularly relevant feature of the pp → V H process is the possibility to study the decay of a Higgs boson into a pair of bottom-antibottom quarks. This decay is extremely important to measure since it provides a direct measurement of the Higgs coupling to fermions, thereby testing the mechanism of fermion mass generation in the SM. Furthermore, since this decay mode dominates the total width of the Higgs boson, the uncertainty on this branching ratio enters into other studies as well, for instance in measurements of the decay of the Higgs boson to invisible final states, which are relevant for dark matter searches [4]. Such a decay is hard to measure in inclusive Higgs production through the leading production modes like the gluon-gluon or vector-boson-fusion channels due to the presence of enormous QCD backgrounds. In the Higgs Strahlung process the presence of a vector boson decaying leptonically provides a clean experimental signature and means experimental analyses related to V H production have a manageable background.
Direct searches for the SM Higgs boson through V H production and H → bb decay has been carried out at the LHC at centre-of-mass energies of 7 TeV, 8 TeV, and 13 TeV. While the use of Run I data at √ s = 7 TeV and 8 TeV by the ATLAS and CMS Collaborations was not able to firmly establish the discovery of the SM-like Higgs boson through this channel [5,6], the use of Run II data at √ s = 13 TeV enabled to do so. In 2017, The LHC experiments [7,8] announced the observation of a SM Higgs-like particle decaying to a pair of bottom-antibottom quarks precisely through this Higgs Strahlung production channel with a significance of 5.6 and 5.3 standard deviations for CMS and ATLAS respectively.
In view of prospective measurements of Higgs Strahlung final states including data from Run II and III at the LHC, it is of crucial importance to have precise theoretical predictions for cross sections and differential distributions in the kinematic regions probed by the experiments. This includes in particular QCD effects in both the production and in the decay of the Higgs boson into a bottom-quark pair.
The present status of theoretical predictions for observables related to V H production with a vector boson decaying leptonically and a Higgs boson decaying into a bottomantibottom quark pair can be summarised as follows.
The total inclusive cross section for associated V H production is known at NNLO QCD precision. It is available through the numerical program VH@NNLO [9] whose ingredients have been reported in refs. [10,11]. The electroweak corrections to the total cross section are known at NLO [12,13]. Differential distributions have also been computed at NNLO QCD, including the computation of H → bb decay at different orders. In refs. [14][15][16], the

JHEP10(2019)002
Higgs decay has been included at NLO while it is included up to NNLO in refs. [17,18]. In addition, for massless bottom quarks, the fully differential decay rate for H → bb known so far at NNLO QCD [19,20] has recently been computed at N 3 LO accuracy in ref. [21], although jet-flavour is not identified in this calculation. Furthermore, a differential computation at NNLO QCD of this Higgs decay rate using massive bottom quarks has been performed in [22]. The combination of fixed-order QCD computations with parton showers have also been the subject of phenomenological studies [23][24][25].
Fully differential NNLO predictions for V H observables obtained via the combination of Higgs production and decay to bottom-antibottom processes have been presented in ref. [17] (for V = Z, W + ) and in ref. [18] (for V = W − ). These computations have essential features in common: at parton level, both consider massless b-quarks except in the Higgs Yukawa coupling and use the same flavour-k t algorithm [26] to define b-jets. Furthermore, the Higgs decay is treated in the narrow-width approximation and the Higgs Yukawa coupling y b is computed at fixed scale µ = m H . Scale variations are only considered in the production sub-process using the central scale choice µ = M V H .
The aforementioned computations differ instead in the theoretical framework employed to regulate infrared divergences at NNLO level: in ref. [17] the q T -subtraction formalism [27] is used for the V H production cross section combined with the CoLoRFulNNLO subtraction method [28] for the H → bb decay. In ref. [18] the nested soft-collinear subtraction scheme [29,30] is used (an extension of the residue subtraction scheme [31]) in both production and decay sub-processes.
It is the aim of this paper to present a computation of V H observables for all three processes (V = Z, W ± ) including NNLO corrections to both production and decay subprocesses. Our goal is to yield a fully differential description of the final states, i.e. including the decays of the vector boson into leptons and the Higgs boson into bottom quarks with off-shell propagators of the vector-and Higgs-boson. The NNLO corrections to both production and decay sub-processes are calculated using the antenna subtraction formalism [32][33][34][35][36][37][38][39][40] implemented within the NNLOJET framework [41].
The structure of this paper is as follows: in section 2, we provide an overview explaining how flavour-dependent observables are computed at fixed-order accuracy within the parton level event generator NNLOJET. A detailed description of the jet-algorithm used to achieve this goal, as well as its application to the V H process are also specified. In section 3, we present the details of the V H calculation giving explicitly the different ingredients appearing in production and decay sub-processes up to NNLO level. Section 4 contains our results for the fiducial cross sections and differential observables related to V H production in pp collisions at √ s = 13 TeV. Those include, for the first time, scale uncertainty estimations related to the separate variation of production and decay scales at each order in α s . Two appendices are enclosed: in appendix A the impact of different criteria for defining the net flavour of jets is studied for a number of relevant NNLO distributions in V H production. Appendix B is dedicated to a comparison of results obtained using two different expressions of the cross section, including either a fixed branching ratio Br(H → bb) as used previously in refs. [16][17][18], or not, as in this paper.

JHEP10(2019)002 2 Flavour tagging of jets
The goal of this work is to provide fixed-order predictions for the hadron-level process pp → ¯ + 2 b-jets + X, i.e. yielding a final state which contains flavour-tagged bottomquark jets (b-jets) and (charged) leptons. The presence of two identified b-jets with a combined invariant mass consistent with m H allows us to associate this final state with the underlying process pp → V H + X → ¯ bb + X. The identification of jet flavour is an essential component of any experimental analysis of this process, which is required to reduce otherwise overwhelming background processes. It is therefore also an integral part of the requirements needed to obtain the corresponding theoretical predictions.
The computation of such observables at fixed order requires the application of a flavoursensitive jet algorithm that -besides reconstructing flavour-insensitive properties such as four-momenta -identifies the flavour of the reconstructed jets based on some welldefined (infrared-safe) criteria [26]. The application of such an algorithm requires a tracking of the flavour of individual partons, which appear in the partonic cross section at each perturbative order.
In the following, we provide a description of how this is achieved within the partonlevel event generator NNLOJET. The discussion is however not specific to the use of the antenna subtraction formalism to regulate infrared divergences occurring in partonic subprocesses beyond LO. In addition, as the application of a flavour-sensitive jet algorithm is not standard (although required from the point of view of massless fixed-order computations) for either theory or experimental communities, we also give a brief overview of the algorithm used for these computations. This section is concluded by providing specific details of the jet clustering implementation relevant for the results presented here regarding the computation of NNLO observables for V H production.

Flavour dressing
The first step towards computing flavour-dependent jet observables is to ensure that the jet algorithm has access to both momentum and flavour information when evaluating the contributions from matrix elements and subtraction terms. To address this issue within NNLOJET, an additional "flavour-dressing" layer that tracks the flavours of all amplitudes as well as reduced matrix elements appearing in subtraction terms has been implemented.
To illustrate how this proceeds, we consider the construction of a generic NLO-type cross section for an n-body final state initiated by the two partons i and j. Following the notation of ref. [40], we may write the contribution to the partonic cross section as where the superscripts R, S, V , T indicate the real, real-subtraction, virtual, and virtualsubtraction terms, respectively.
As an example of the flavour-dressing procedure for the amplitudes, we consider the real-emission cross section (omitting the sum over potential colour orderings) which takes JHEP10(2019)002 We denote the final-state symmetry factor by S n+1 , the normalisation factor (which includes constants, couplings, colour factors) by N R NLO , the 2 → n + 1 particle phase space by dΦ n+1 , and the momentum of the partons i, j by p 1,2 . The partial squared amplitude M 0 n+3 is evaluated with the momentum set {p n+3 } and a corresponding flavour set {f n+3 }. The flavour-sensitive jet algorithm J (n+1) n builds n jets from n + 1 final-state partons which carry momentum and flavour labelled by the sets {p n+1 } and {f n+1 } respectively. The real subtraction cross section can be written in a similar fashion: where the index k runs over all possible unresolved partons in dσ R ij,NLO and X 0 3 (·, k, ·) denotes the three-parton antenna function that factorises from the associated reduced squared matrix-element M 0 n+2 . In this case, the jet algorithm acts upon mapped finalstate momentum and flavour sets {p n } and {f n } associated with the reduced squared matrix element M 0 n+2 . As the total subtraction cross section must take into account all possible unresolved limits of parton k, this cross section may be composed of multiple flavour structures. The subtraction method is only effective if the evaluation of flavourdependent observables in both the real and real-subtraction cross sections match in all possible unresolved limits. This is only ensured if an infrared-safe flavour-sensitive jet algorithm is applied.
To construct the NLO cross section according to eq. (2.1), a similar procedure must also be applied to both virtual and virtual-subtraction (in the antenna formalism, these include integrated subtraction and mass-factorisation contributions) cross sections. This construction is obtained in a similar fashion, by tracking both the momentum and flavour sets associated to all partial squared amplitudes and reduced squared matrix elements appearing in these contributions and then applying the flavour-sensitive jet algorithm to the subset of final-state particles within these sets. To allow the computation of flavourdependent jet observables at NNLO, the same ideas extend to one order higher and this flavour-dressing procedure is applied to all NNLO-type parton level contributions and their corresponding subtraction terms.

Flavoured-jet algorithm
Throughout this work jets are reconstructed with the flavour-k t algorithm, which provides an infrared-safe definition of jet flavour. The main difference with respect to a native jet algorithm is that the clustering of particles relies on both momentum and flavour information of the input pseudo-jets. For completeness, we summarise the main steps of JHEP10(2019)002 the algorithm for hadron-hadron collisions originally presented in ref. [26] (also summarised in ref. [42]).
The algorithm proceeds by assigning a net flavour to all pseudo-jets or jets based on their quark flavour content, attributing +1 (−1) if a quark (antiquark) of the flavour under consideration is present. In an experimental context, the presence of a quark flavour could be inferred from a fully/partially reconstructed hadron. A criterion is then applied to these objects to determine if they carry flavour, possible examples being: the net flavour (sum of quarks and antiquarks); or the net flavour modulo two. Objects are considered to carry flavour if they carry non-zero values of this criterion. The algorithm then proceeds by constructing distance measures for pairs of all final-state pseudo-jets i and j (d ij ) as well as beam distances (d iB and d iB ). These (flavour-dependent) distances are defined as In these definitions, k ti and k tj are the transverse momentum of the pseudo-jets i and j, and the rapidity difference and azimuthal angular separation between these pseudo-jets is given by ∆y ij and ∆φ ij , respectively. The parameters R and α define a class of measures for the algorithm. The (rapidity-dependent) transverse momentum of the beam B at positive rapidity k tB , and beamB at negative rapidity k tB , are defined as: with Θ(0) = 1/2 and the index i going over all pseudo-jets. While this flavour-aware jet algorithm is substantially more complex than the flavourblind anti-k t algorithm [43], its use is unavoidable in fixed-order computations based on massless quarks. At NLO, the flavour criterion of a pseudo-jet ensures that a collinear splitting of the form g → qq is indistinguishable from a gluon (or flavourless) jet. The subtraction formalism presented in eq. (2.1) would already be spoiled without this criterion. At NNLO, the flavour-dependent distance measure in eq. (2.4) ensures that a pair of flavoured quarks originating from a wide-angle gluon splitting is clustered into a pseudo-jet before being combined with any other (harder) pseudo-jets. This avoids the situation where one of these soft quarks may be clustered with a hard pseudo-jet that carries zero flavour, which would lead to a definition of jet flavour sensitive to soft physics. These are issues which are otherwise insurmountable for fixed-order computations involving massless quarks.
The flavour-k t algorithm described above is available in the NNLOJET framework and has been validated against an independent implementation using FastJet [44,45].

Jet clustering for the V H process
The discussion of flavour dressing and the jet algorithm presented in this section are quite general and are applicable to all processes implemented within NNLOJET. Here we discuss a few specific points related to the application of the jet algorithm to the V H process, which will be relevant to the results presented in later sections of this paper.
The first point is that we wish to apply the flavour-k t algorithm to the partonic process qq → V H → ¯ bb, including NNLO QCD corrections which will be discussed in section 3. When higher-order corrections are considered, additional light or b-quark partons can be emitted from both production and the decay sides, as illustrated in figure 1(b). The jet clustering is performed by considering b-quarks to be flavoured (all other partons carrying zero flavour) and fully accounting for emissions from both production and decay during the jet clustering process. While our calculation focusses on the decay sub-process H → bb, it has been implemented in such a way that predictions for the hadronic process pp → ¯ + 2 c-jets + X can also be produced. This may be interesting in view of possible future measurements by the LHCb Collaboration [46].
The second point is that the definition of the transverse momentum of the beam is altered to account for the presence of a leptonically decaying gauge-boson. This is done by modifying eq. (2.6) according tõ where E t,V and y V are the transverse energy and rapidity of the reconstructed gauge-boson. A similar modification to the beam transverse momentum at negative rapidity (2.7) is assumed. This modification is introduced to provide a better estimate of the hardness of the beam, which can affect the clustering outcome. One could alternatively modify the beam hardness by including the charged leptons, which may be necessary in experimental situations where the gauge-boson cannot be fully reconstructed. The final point is related to our choice of flavour criterion during the clustering process. We have chosen to define the flavour of pseudo-jets to be the net-flavour of its constituents JHEP10(2019)002 modulo two, which means that all pseudo-jets which contain an even flavour content are considered to have zero net-flavour. The motivation for this choice is that, in our opinion, it is the most feasible realisation of the flavour-k t algorithm experimentally. Focussing on the case of b-jets, the main consideration is that most experimental approaches to flavour tagging are sensitive only to the absolute flavour [47][48][49] (and do not additionally charge tag the jets). All implementations of the algorithm must consider the combination of a bb-quark pair (or equivalently a BB-hadron pair) as carrying zero flavour, as required to guarantee its infrared safety as discussed above. Therefore, in the absence of charge tagging, any (pseudo)-jet which contains the presence of an even number of b (B) and/or b (B) quarks (hadrons) should also be considered to carry zero flavour, as experimentally these signatures are indistinguishable.
The charge-tagging of flavoured jets is also possible [50], for example in the presence of semi-leptonic B-hadron decays. However, the drawback is a large reduction in event statistics (roughly an order of magnitude for each b-jet, as the branching fraction Br(B → + X) ≈ 10%) with little informational gain. Accordingly, to present our results for NNLO observables related to V H production in section 4, we shall use the version of the flavour-k t algorithm where all even-tagged (pseudo)-jets carry zero flavour. We further provide an examination of the impact of the even-tag exclusion in the shape and normalisation of flavour sensitive observables in appendix A.

Details of the calculation
In this section we present the main ingredients that enter the calculation of the Higgs Strahlung process at NNLO. We establish how those building blocks are assembled to express the cross section in a factorised form in terms of production and decay sub-processes.

General framework
We consider the process pp → V H + X → ¯ bb + X where the vector boson (V ) decays leptonically and the Higgs boson (H) decays into a pair of bottom quarks bb. We compute NNLO QCD observables related to these reactions by including corrections up to order α 2 s in both production and decay sub-processes. This enables us to express the fully differential cross section at the kth order in a factorised form given as The term dσ V H , which corresponds to the production part, comprises the vector propagator and the leptonic decay V → ¯ , including spin correlations between the initial-state partons and the final-state leptons. The term denoted by dσ (j) H→bb corresponds to the decay part and includes the Higgs propagator and its subsequent decay to a bottom-antibottom quark pair. We treat all light quarks as massless including the bottom quark with the exception of the Yukawa coupling mediating the H → bb decay. In the decay the bottom quark Yukawa coupling to the Higgs boson is renormalised in the MS scheme at the scale µ dec. , taken to be JHEP10(2019)002 proportional to the Higgs-boson mass m H . 1 Note that the factorised form of the associated Higgs production cross section (3.1) does not include interferences between production and decay sub-processes. This is a valid approximation owing to the smallness of the Higgs decay width, which further formed the basis of the narrow-width approximation (NWA) as used in previous calculations.
Up to O(α 2 s ), the cross section may then be written as We note that this formulation of the NNLO cross section does not contain the Higgs-boson branching ratio into b quarks given as Br(H → bb) = Γ H→bb /Γ H . This is in contrast to previous calculations for the V H process at NNLO, either considering the decay sub-process at NLO [14][15][16] or NNLO [17,18], which all employed a scaled variant of the cross section incorporating the Higgs-boson branching ratio at a fixed value and thus not subject to an α s expansion. It is worth mentioning that this scaled variant of the cross section was essential in describing the data using fixed-order predictions at LO and NLO [14,15]. With this formulation, the LO predictions have the correct normalisation; NLO corrections are kept small and have a small residual theoretical uncertainty. If computed up to order α 2 s , we here argue that the need of such scaling factors in the formulation of the cross section becomes questionable.
In appendix B, we will further elaborate on this matter and compare the results obtained with both approaches for the fiducial cross sections up to order α 2 s . We find that a consistent treatment of theoretical uncertainties outweighs the precision gain that one might obtain by scaling the cross section to a fixed branching ratio, if the cross section is computed including NNLO corrections in both production and decay parts. This further motivates the simpler formulation of the cross section given above in eq. (3.2) where no scaling factors are applied. This will be our default setup throughout this work and for the results presented in section 4.
As a validation of our calculation, we performed a comparison to the results of ref. [18] by adopting their calculational setup and found perfect agreement with the reported values for the total cross sections in table I at each perturbative order in α s .

Production and decay parts up to O(α 2 s )
Based on our master formula (3.2) for the V H process at NNLO, we specify the individual ingredients of the production and decay parts in the following and describe how they are assembled to the final prediction for the Higgs Strahlung process.

Production parts
Up to order α s , only one type of contribution enters the associated Higgs production cross section, which is given by Drell-Yan-like diagrams with a subsequent Higgs emission from the gauge-boson leg. Starting from O(α 2 s ), additional quark-loop induced contributions arise. These can be treated independently from the aforementioned Drell-Yan-type ones as the relevant Feynman amplitudes are individually gauge invariant. In the following, we describe these two production modes one after the other.
Drell-Yan-type: these contributions arise from the Drell-Yan-like production of a virtual W ± or Z boson, which then splits into a real vector boson and a Higgs particle. In our calculation we include them up to O(α 2 s ) using off-shell amplitudes that effectively treat both the directly produced vector boson and the vector boson that decays leptonically as virtual particles. Representative Feynman diagrams for this production mode are illustrated in figure 2 at LO (a) and NNLO (b).
These contributions only involve the square of Drell-Yan-like amplitudes and the infrared singularities are dealt with using the NNLO antenna subtraction formalism [32]. The subtraction terms can be readily constructed based on the NNLO calculation for the Drell-Yan processes, which are available within the NNLOJET framework.
Top-quark-loop induced: starting from O(α 2 s ), new types of diagrams induced by quark loops must be taken into account for the V H production process. Depending on whether the gauge boson and/or the Higgs boson couple to the quark loop, these contributions can be classified into three categories:   Figure (a) depicts an R I -type amplitude, which is present for both ZH and WH production channels. Figures (b), (c) illustrate representative gluon-gluon induced heavy quark loop amplitudes, which are only present for ZH production. In case (c) on the other hand, the Higgs boson does not couple directly to the quark loop and we have to consider all quark flavours inside the loop. For the quarks of the first two generations: q = u, d, s, c, the corresponding amplitudes cancel due to the fact that an equal count of up-and down-type quark flavours are evaluated. This cancellation is spoiled in the case of the third generation due to the non-vanishing mass of the top quark. As a result, both the top-and bottom-loop components are included.
The complete O(α 2 s ) top-loop-induced contributions were computed for on-shell vector bosons in ref. [11], relying in some cases on the infinite-top-mass approximation. In our NNLO calculation we include those that are known in the exact theory and numerically sizeable but omit those which are only known in the infinite-top-mass limit. Specifically, for the NNLO contributions associated with diagrams (a), we include diagrams with topquark-loop insertions onto an external gluon line. The related amplitudes are referred to JHEP10(2019)002 as R I in ref. [11] and have been included in the previous computations [17,18]. 2 Regarding the amplitudes of class (b) and (c), which are exclusively present in ZH production, we only include the gluon-gluon-induced channels shown in figures 3(b),(c). Phenomenologically, they represent the dominant component among the top-loop-induced contributions due to the large gluon luminosity at the LHC and were also considered in the previous calculations at NNLO. 3 The heavy-quark-loop-induced contributions included in our calculation have been either independently rederived or implemented using known results, in particular those given in ref. [16]. A validation of the implementation was performed against OpenLoops amplitudes [51] and full agreement was found in all cases.

Decay parts
For the decay sub-process H → bb, we required the corrections up to O(α 2 s ) as indicated in our master formula (3.2). The corresponding amplitudes at one-and two-loop level were obtained from the analytic expressions of refs. [19,20] and were decomposed into different colour levels according to antenna formalism conventions. A validation of all one-loop amplitudes was performed against the OpenLoops library [51], yielding full agreement. In addition, subtraction terms capturing the infrared singularities are required. Those have been constructed for the Higgs decay up to order O(α 2 s ) for the present computation. Checks for the correct divergent behaviour in all single-and double-unresolved limits have been performed in order to ensure the proper cancellation of singularities in the realemission corrections as well as the cancellation of poles against the virtual amplitudes.
The decay sub-process up to NNLO only enters in eq. (3.2) when combined with the Drell-Yan-type production parts. For the top-quark-loop induced contributions, which are already of O(α 2 s ), the decay only needs to be considered at tree level. 2 We did not include the two-loop amplitudes from this class as they are currently not known in the full theory but only in the infinite-top-mass limit. Diagrammatically, these would be given by and are referred to as VI in ref. [11]. The numerical impact to the total NNLO cross section is estimated to be below the percent level. This contribution was omitted in ref. [17] but kept in ref. [18]. 3 The contributions that we omitted from this class are are given by diagrams of the following type: They are denoted as RII and VII respectively in ref. [10]. The one-loop amplitude RII is known in the full theory but it merely constitutes a sub-permille effect. The two-loop amplitude VII is currently only known in the large-top-mass limit but its impact is estimated to be at the sub-percent level. These contributions were also omitted in ref. [17].

JHEP10(2019)002 4 Numerical results
In this section we present phenomenological results obtained for the different V H processes using our implementation in the parton-level event generator NNLOJET. We first summarise the general setup in section 4.1 and move on to discuss the integrated fiducial cross sections obtained within this setup in section 4.2. We devote section 4.3 to validating the scale dependence of our numerical results and present differential distributions for flavour-sensitive observables in section 4.4.

General setup
We provide predictions for proton-proton collisions at √ s = 13 TeV using the parton distribution function NNPDF31_nnlo_as_0118 provided via the LHAPDF library [52]. Each event was required to contain at least two b-jets with transverse momentum p ⊥,b > 25 GeV and rapidity |y b | < 2.5. Charged leptons were required to have a transverse momentum above p ⊥, > 15 GeV and for their rapidity to satisfy |y | < 2.5. For the W ± H processes, we additionally demanded a minimum missing transverse energy of E ⊥,miss > 15 GeV to identify the neutrino in the final state. We used the flavour-k t algorithm with an even-tag exclusion to cluster b-jets as described in sections 2.2 and 2.3 with the parameters R = 0.5 and α = 2.
We employed the G µ -scheme for the electroweak input parameters and the full set of independent parameters entering the computation are given by The running of the strong coupling (α s ) was evaluated using the LHAPDF library with the associated PDF set, while the MS mass of the bottom quark (m b ) was directly computed within NNLOJET. Finally, in the case of W ± H production, we assumed a diagonal CKM matrix for the vector-boson-quark couplings. For the unphysical scales appearing in the calculation, we chose to set and vary them independently for the production and decay parts. The central factorisation and renormalisation scales of the production sub-processes were chosen to the invariant mass of the V H system M V H , whereas the central renormalisation scale of the decay was set to the Higgs-boson mass m H . We evaluate the differential cross section for a total of 21 different scale settings that are obtained from all possible combinations of with the additional constraint 1 2 ≤ µ F /µ prod. R ≤ 2 following the conventional 7-point scale variation for the production sub-process.

Fiducial cross section
The cross-section predictions including fiducial cuts for the different V H processes are summarised in table 1  Note that the reduction of scale uncertainties observed here is spoiled in all cases when a rescaling prescription is employed that incorporates a fixed branching ratio for the H → bb decay, as is commonly done in previous calculations for the V H processes. A comparison of our results in table 1 to such a rescaled cross-section prediction is presented in appendix B.

Scale variations
The dependence on the renormalisation scales µ prod./dec. R can serve as a non-trivial check of the final results obtained from the numerical computation. To this end, we ensure that the different scale settings of eq. (4.2) are correctly reproduced by the analytic renormalisationgroup running starting from the central scale choice. 4 This is of particular importance for the calculation at hand, as the independent variation of scales for the different sub-processes was for the first time implemented in the NNLOJET framework for the present work.

JHEP10(2019)002
The comparison between the analytic evolution and the 21 points obtained from the numerical computation using NNLOJET are shown in figure 4 for the case of the W + H process at NLO (a)-(c) and NNLO (d)-(f). We performed a scan in the two-dimensional (µ prod. R , µ dec. R ) space by choosing three different slices that cover the combinations in eq. (4.2) where the three choices in the factorisation scale µ F = M WH 1, 1 2 , 2 are shown as separate curves: (a,d) We keep the decay renormalisation scale fixed to µ dec. R = m H and vary the scale in the production sub-process according to We keep the production renormalisation scale fixed to µ prod. R = M WH and vary the scale in the decay sub-process according to (c,f) We choose a diagonal slice in the (µ prod. R , µ dec. R ) plane setting K prod.
R ≡ K R with the individual scales given as Note that the invariant mass M WH constitutes a dynamical quantity that varies on an event-by-event basis. The curves in figure 4 are obtained by picking a specific bin M WH ∈ [220, 230] GeV to assign a value to the production scale, where the width of the bands in the smooth curves correspond to the uncertainty that arises from the finite bin width. We observe an excellent agreement between the numerical results from NNLOJET and the curves predicted from the renormalisation group equations. The dramatic reduction in scale uncertainties can be further appreciated by contrasting the vertical ranges in the figures at NLO (left) and NNLO (right). We carried out the same tests also for the W − H and the ZH processes as well as for other individual M V H bins in the distributions and found that the scale variation of the numerical results match the analytical formulae in all cases.

Distributions
In figures 5-7 we present differential distributions of flavour-sensitive observables for the three different associated Higgs boson production processes W + H, W − H, and ZH: 5 where in (b-d) the two b-jets are selected whose invariant mass is closest to m H in order to identify the candidate pair that is most likely to originate from the Higgs decay.

JHEP10(2019)002
G GS E >IE*H9@    contributed at NNLO, which also explains the widening of the theoretical uncertainty bands around the threshold regions of these distributions. Concerning the invariant mass distribution of all three production modes shown in figures 5(d)-7(d), the features previously noted in refs. [17,18] can be confirmed by our predictions as well: due to the very narrow width of the Higgs boson, the m bb distribution has a natural kinematic threshold at m H = 125.09 GeV and the phase space away from this value is barely populated at LO. Consequently, NNLO corrections are effectively NLOaccurate for most of the bins, which explains the large corrections and relatively larger uncertainty bands for these distributions. The left shoulder below m H is mainly the result of radiation from the decay, whereas the shoulder above m H is due to radiative corrections to the production. Fixed-order predictions at the threshold region of m bb ≈ m H , however, should not be trusted as they are prone to Sudakov-type instabilities. A proper treatment of this region would require the inclusion of resummation effects. In our case, the binning is sufficiently coarse so that such instabilities only manifest in larger uncertainty bands for the m bb = m H bin and not as an explicit divergence.

Summary and conclusions
We reported on the calculation of NNLO corrections to the Higgs Strahlung processes W + H, W − H, and ZH including the off-shell leptonic decay of the gauge boson as well as the Higgs decaying into a bottom-antibottom pair. The calculation consistently takes into account NNLO corrections to the production and decay sub-processes and fully retains the differential information on the final state.
The study of V H (H → bb) processes critically relies on the tagging of bottom jets in order to isolate the candidate pairs associated to the Higgs boson. We described our independent implementation of the infrared-safe flavour-k t algorithm in the NNLOJET partonlevel event generator and the necessary modifications this entails in the framework of the antenna subtraction formalism.
A detailed account was given on the residual theory uncertainties by allowing the scales in the production and decay sub-processes to vary independently. This conservative approach resulted into taking the envelope of 21 scale variations for the full process but allowed for a more comprehensive view into the impact of higher orders on the reduction of scale uncertainties. The NNLO corrections to the fiducial cross section were found to exhibit a good perturbative convergence with residual uncertainties at the percent level. We contrasted our naïve perturbative expansion of the cross section with a more commonly employed rescaling procedure using the branching ratio BR(H → bb), where we observed the latter to overestimate the residual scale uncertainties. This was attributed to a miscancellation in the scale dependence among the terms that receive different rescaling factors and lead us to advocate the simpler prescription to be more reliable beyond NLO.
Flavour-sensitive observables were studied by investigating differential distributions where a similar stabilisation of the perturbative series was found as in the cross sections. Larger effects from higher-order corrections were seen in the invariant mass distributions of two b-jets, which can be attributed to this observables being only NLO-accurate away JHEP10(2019)002 Even-tag exclusion 20.6828 ± 0.0055 Original flavour-k t 20.7093 ± 0.0063 Ratio 99.87% from m bb ∼ m H . A comparison between the W ± H and ZH processes showed a qualitatively similar behaviour but also emphasised the phenomenologically sizeable impact that arise from the gluon-gluon-induced top-quark loop amplitudes.
The study of flavour-sensitive jet observables with fixed-order predictions, such as those associated to b-jets in the present work, must be performed in an infrared-safe way. For calculations based on massless QCD this can only be achieved with a flavour-aware jet algorithm (such as flavour-k t ), while for a massive calculation this is achievable with a flavour-blind algorithm (such as anti-k t ). In many cases the corresponding massive calculation may not be available, or the massless calculation may actually be preferred due to the presence of large logarithmic corrections which may be easily resummed via PDF evolution. Future comparisons to measurements are only viable if a similar prescription is also employed in the experiment, and the application of the even-tag exclusion here was mainly motivated to facilitate the experimental implementation. The use of flavour-sensitive jet algorithms is not only important to the V H process class but we expect it to be of relevance for any flavour-sensitive jet observable, such as the associated production of the flavoured jet with a gauge boson. Such studies will be left for future work. with the criterion that the flavour of (pseudo)-jets is assigned as the net flavour of its constituents modulo two, which we believe is more motivated from an experimental point of view as discussed in section 2.3.
To investigate the impact of this "Even-tag exclusion" on the fixed-order predictions, we have re-computed the fiducial cross-section and distributions reported in section 4.2 and 4.4 without the additional "modulo two" criterion -we refer to these results as "Original flavour-k t ". This impact of the choice of this criterion is visualized in the case of the W + H process in figure 8 for the p ⊥,b , p ⊥,bb , and ∆R bb distributions. In that figure, the ratio of the two NNLO central values are divided bin-by-bin, demonstrating that this choice has no overall effect on the shape of these distributions. The small variation between bins can be attributed to statistical fluctuations. This behaviour is also confirmed at the level of the fiducial cross section as reported in table 2, where the results are consistent within statistical uncertainties. This supports our claim that no significant portion of the events are discarded by switching to the even-tag-excluded version of flavour-k t in our fixed-order predictions.

B Comparison with previous formulations
As mentioned in section 3.1, NNLO-accurate observables for associated Higgs production have also been presented in [16][17][18]. However, the cross section in these calculations is assembled in a different manner compared with our expression in eq. (3.2). Specifically, the Higgs decay at the different orders had been scaled up to a fixed value of the accurately  known branching ratio of the H → bb process. Up to NNLO, the cross sections in this formulation is assembled as follows: The branching ratio Br(H → bb) is kept fixed and is not a subject to an α s expansion.
In the following, we elaborate on possible drawbacks that this prescription entails, in particular concerning theory uncertainties as estimated through scale variations.
Firstly, the scaling factors effectively divide out the Yukawa coupling y b (µ dec. R ) ∝ m b (µ dec. R ) from the prediction. As a result, any running of the mass as induced by the MS scheme exactly cancel in the final result. This can lead to underestimating the uncertainties, which is especially apparent at LO where the scale dependence of the Yukawa coupling otherwise dominates the uncertainties.
Secondly, analysing the structure of the scaled cross sections at NLO (B.2) and NNLO (B.3), it is apparent that they are assembled as a sum of terms where different scaling factors K (i) accompany the different perturbative coefficients of the production cross section dσ (j) V H . This mismatch can interfere with the compensation mechanism between terms of different orders, possibly distorting the theory error estimate obtained through variations of the production scale µ prod. R .
To quantify the differences between the two approaches, in table 3 we report the fiducial cross sections obtained according to (B.1)-(B.3) using Br(H → bb) = 58.09% [54]. Comparing these predictions with those given in table 1 using the unscaled cross section JHEP10(2019)002 formulae (3.1), we observe that the central value of the LO prediction is substantially improved in the scaled predictions thanks to absorbing higher-order effects to the H → bb decay through the branching ratio. At NLO, however, the scaled prediction appears to slightly overestimate the cross section, while the associated theory uncertainties are comparable in size between the two formulations. At NNLO, both prescriptions agree well in their respective central values, however, sizeable differences can be seen in their associated uncertainties. The scaled predictions at NNLO show almost no reduction in scale uncertainties -even increasing for ZH -compared to the respective NLO number, whereas our formulation (3.1) exhibits a substantial reduction in scale uncertainties when going from NLO to NNLO. This difference can be attributed to the aforementioned compensation of scale dependences, which is spoiled by the different rescaling factors used in eq. (B.3).
The effects of dividing out the Yukawa coupling in the decay and the scaling factor mismatch during the assembly of production cross sections are apparent as the theoretical uncertainties of the NNLO cross section barely change compared to their NLO values. In our opinion, the consistent treatment of theoretical uncertainties outweighs the precision gain that one might (or might not) get by scaling to a fixed branching ratio, especially in the case of NNLO-accurate observables. This further motivates our initial and simpler formulation we presented in eq. (3.1) of section 3.1 where no scaling factors are applied.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.