Associated production of a Higgs boson decaying into bottom quarks at the LHC in full NNLO QCD

TL;DR

This paper presents a comprehensive NNLO QCD calculation for Higgs production with bottom quarks and vector bosons at the LHC, including decay effects, which improves the precision of theoretical predictions for experimental comparisons.

Contribution

It provides the first fully exclusive NNLO QCD calculation for Higgs production with bottom quarks in association with W/Z bosons, including decay and spin correlations.

Findings

NNLO QCD corrections significantly reduce the accepted cross section.

Corrections substantially alter the shape of key distributions.

Inclusion of decay and spin effects is crucial for accurate predictions.

Abstract

We consider the production of a Standard Model Higgs boson decaying to bottom quarks in association with a vector boson W/Z in hadron collisions. We present a fully exclusive calculation of QCD radiative corrections both for the production cross section and for the Higgs boson decay rate up to next-to-next-to-leading order (NNLO) accuracy. Our calculation also includes the leptonic decay of the vector boson with finite-width effects and spin correlations. We consider typical kinematical cuts applied in the experimental analyses at the Large Hadron Collider (LHC) and we find that the full NNLO QCD corrections significantly decrease the accepted cross section and have a substantial impact on the shape of distributions. We point out that these additional effects are essential to obtain precise theoretical predictions to be compared with the LHC data.

Tables

Table 1. Table 1: Cross sections and their scale uncertainties for p ​ p → V ​ H + X → l 1 ​ l 2 ​ b ​ b ¯ + X → 𝑝 𝑝 𝑉 𝐻 𝑋 → subscript 𝑙 1 subscript 𝑙 2 𝑏 ¯ 𝑏 𝑋 pp\rightarrow VH+X\rightarrow l_{1}l_{2}b{\bar{b}}+X at LHC with s = 13 𝑠 13 \sqrt{s}=13  TeV. The applied kinematical cuts are described in the text.

$\sigma$ (fb)	NNLO(prod)+NLO(dec)	full NNLO }
$pp\to W^{+}H+X\to l{\nu }_{l}b\overline{b}+X$	${3.94}_{-1.5\%}^{+1\%}$	${3.70}_{-1.5\%}^{+1.5\%}$ }
$pp\to ZH+X\to \nu \nu b\overline{b}+X$	${8.65}_{-3.5\%}^{+4.5\%}$	${8.24}_{-3.5\%}^{+4.5\%}$ }

Equations

$$d{\sigma }_{h_1{h}_2\to VH\to Vb\bar{b}}=d{\sigma }_{h_1{h}_2\to VH}\times \frac{d{\Gamma }_{H\to b\bar{b}}}{{\Gamma }_H}=\left[\sum _{k=0}^{\infty }d{\sigma }_{h_1{h}_2\to VH}^{(k)}\right]\times \left[\frac{{\sum }_{k=0}^{\infty }d{\Gamma }_{H\to b\bar{b}}^{(k)}}{{\sum }_{k=0}^{\infty }{\Gamma }_{H\to b\bar{b}}^{(k)}}\right]\times \text{mbox}Br(H\to b\bar{b}),$$

$$d{\sigma }_{h_1{h}_2\to VH\to Vb\bar{b}}^{\mathrm{NNLO}}$$

$$d{\sigma }_{h_1{h}_2\to VH\to Vb\bar{b}}^{\mathrm{NNLO}(\mathrm{prod})+\mathrm{NLO}(\mathrm{dec})}$$

Full text

TIF-UNIMI-2017-5

**Associated production of a Higgs boson decaying **

into bottom quarks at the LHC in full NNLO QCD

Giancarlo Ferrera(a), Gábor Somogyi(b)  and  Francesco Tramontano(c)

**

(a) TIF Lab, Dipartimento di Fisica, Università di Milano and

INFN, Sezione di Milano, I-20133 Milan, Italy

(b) MTA-DE Particle Physics Research Group,

H-4010 Debrecen, PO Box 105, Hungary

(c) Dipartimento di Fisica, Università di Napoli and

INFN, Sezione di Napoli, I-80126 Naples, Italy

Abstract

We consider the production of a Standard Model Higgs boson decaying to bottom quarks in association with a vector boson $W^{\pm}/Z$ in hadron collisions. We present a fully exclusive calculation of QCD radiative corrections both for the production cross section and for the Higgs boson decay rate up to next-to-next-to-leading order (NNLO) accuracy. Our calculation also includes the leptonic decay of the vector boson with finite-width effects and spin correlations. We consider typical kinematical cuts applied in the experimental analyses at the Large Hadron Collider (LHC) and we find that the full NNLO QCD corrections significantly decrease the accepted cross section and have a substantial impact on the shape of distributions. We point out that these additional effects are essential to obtain precise theoretical predictions to be compared with the LHC data.

May 2017

The discovery of the long sought Higgs boson ($H$) [1, 2] by the ATLAS and CMS Collaborations at the Large Hadron Collider (LHC) [3, 4] paved the way for the experimental investigation of the electroweak (EW) symmetry breaking mechanism and, in particular, for the measurements of the Higgs boson couplings to the Standard Model (SM) particles. In this respect the increasing amount of precise experimental data collected at the LHC demands a corresponding improvement of theoretical predictions.

One of the main production mechanisms of the SM Higgs boson at hadron colliders is the associated production with a vector boson ($V=W^{\pm},Z$). This process offers the unique opportunity to study both the Higgs boson coupling to massive gauge bosons and to bottom ($b$) quarks via the decay $H\rightarrow b\bar{b}$.

A direct search for the SM Higgs boson through associated *VH * production and $H\rightarrow b\bar{b}$ decay has been carried out at the LHC at a centre–of–mass energy of $\sqrt{s}=7/8$ TeV [5, 6] and at $\sqrt{s}=13$ TeV [7, 8]. The ATLAS and CMS Collaborations observed an excess of events above the expected background with a measured signal strength relative to the SM expectation of $0.90\pm 0.18\mathrm{(stat.)}^{+0.21}_{-0.19}\mathrm{(syst.)}$ [7] and $\mu=1.06^{+0.31}_{-0.29}$[8] respectively.

High precision theoretical predictions require detailed computations of radiative corrections for cross sections and corresponding distributions. The total cross section for associated *VH * production is known at next-to-next-to-leading order (NNLO) in QCD [9, 10, 11, 12] and next-to-leading order (NLO) in the electroweak theory [13, 14]. Fully differential calculations have been performed in NNLO QCD for the *VH * production cross section, together with the NLO QCD corrections for the Higgs boson decay rate into bottom quarks [15, 16, 17, 18]. The fully differential $H\rightarrow b{\bar{b}}$ decay rate has been computed up to NNLO in QCD [19, 20] while the inclusive rate is known up to ${\mathcal{O}}(\alpha_{S}^{4})$ [21] and up to NLO in the electroweak theory [22, 23]. Resummation and higher order (beyond NNLO) QCD effects have been investigated in Refs.[24, 28, 25, 27, 26, 29, 30] while the combination of fixed-order QCD calculations with parton shower Monte Carlo algorithms has been considered in Refs. [31, 32, 33]. NLO EW effects at fully differential level for *VH * production with the leptonic decay of the vector boson have been considered in Refs. [14, 34].

In this letter we present the fully differential calculation of the NNLO QCD corrections for the production cross section and for the Higgs boson decay rate to bottom quarks, exploiting the very good accuracy of the narrow width approximation for the Higgs boson ($\Gamma_{H}\ll m_{H}$). In Refs. [35, 16] it was shown that, when the set of kinematical cuts applied in the LHC analyses are considered, the effect of QCD corrections to the Higgs boson decay process can be large. Motivated by these findings, we extend existing calculations on higher order QCD predictions by considering the complete second order terms. Together with the NNLO corrections for the production cross section, we include the NNLO corrections to the $H\rightarrow b\bar{b}$ decay rate and the combination of the NLO terms for the production and decay stages.

We implemented our computation in the parton level Monte Carlo numerical program HVNNLO which allows the user to apply arbitrary kinematical cuts on final-state leptons, $b$ jets and associated QCD radiation, and to compute the corresponding distributions in the form of histograms.

The main result of our study is that for a typical set of kinematical cuts applied in the LHC analyses we observe a substantial decrease of the complete NNLO QCD prediction with respect to lower order calculations. Therefore the inclusion of the QCD effects we have calculated could be relevant to improve the agreement between the SM predictions and the current LHC data. In this letter we discuss the main ingredients of our computation, a more comprehensive phenomenological analysis will appear elsewhere.

We consider the inclusive hard scattering reaction $\,h_{1}+h_{2}\rightarrow VH+X\rightarrow l_{1}l_{2}b\bar{b}+X\,$, where the collision of the hadrons $h_{1}$ and $h_{2}$ produces the ${\it VH\,}$ system ($V=W^{\pm},Z$) which subsequently decays into the lepton pair $l_{1}l_{2}$ ($l_{1}l_{2}\!\equiv\!l\nu_{l}$ in the case of $W^{\pm}$ decay) and the bottom quark-antiquark pair $b\bar{b}$, while $X$ denotes the accompanying QCD radiation. We consider a high value of the ${\it VH\,}$ invariant mass ($M_{{\it VH\,}}$), which sets the hard-scattering scale of the process, and we treat the colliding hadrons, the leptons and the $b$ quarks in the massless approximation.

By using the narrow width approximation for the Higgs boson, the perturbative QCD expansion of the fully differential cross section can be written in the following factorized form †††In order to simplify the notation the leptonic decay $V\rightarrow l_{1}l_{2}$ of the $V$ boson (including spin correlations) has been understood, since it has no effect from the point of view of QCD corrections.:

[TABLE]

where $\Gamma_{H\rightarrow b\bar{b}}$ and $\Gamma_{H}$ are the Higgs boson partial decay width to bottom quarks and the total decay width respectively, and the expansion in powers of $\alpha_{S}$ is given by the exponent $k$. Eq. (1) is arranged in a form such that we can exploit the precise prediction of the Higgs boson branching ratio into $b$ quarks $\mbox{Br}(H\rightarrow b\bar{b})=\Gamma_{H\rightarrow b\bar{b}}/\Gamma_{H}$ (see for instance Ref. [36]), by which we normalize the contributions to the differential decay rate of the Higgs boson ‡‡‡Indeed, by considering observables that are inclusive over the Higgs boson decay products, we obtain the production cross section times the branching ratio $d\sigma_{h_{1}h_{2}\rightarrow V\!H}\times{\mbox{Br}}(H\rightarrow b\bar{b})$..

By expanding Eq. (1) up to the second order in $\alpha_{S}$ we have:

[TABLE]

Eq. (S0.Ex2) contains the complete NNLO contributions, which include the $\alpha_{S}^{2}$ corrections: (i) to the decay rate (included in the first term on the right hand side), (ii) from the combination of the NLO contributions for production and decay (included in the second term on the right hand side), and (iii) to the production cross section (the third term on the right hand side). The novel result of this letter compared with previous approximations concerns the full computation of the terms (i) and (ii).

In order to compare with the partial NNLO calculations considered so far [16, 17, 18], we also consider the truncation of Eq. (1) defined as:

[TABLE]

which contains the NNLO corrections for the production cross section together with the NLO corrections for the $H\rightarrow b{\bar{b}}$ decay rate.

Each term in the Eqs. (S0.Ex2) and (S0.Ex3) includes all the relevant contributions from (double-) real, real-virtual and (double-) virtual corrections. In our implementation we have employed the $q_{T}$ subtraction method for the *VH * production cross section [15, 17] and the CoLoRFulNNLO method for the $H\rightarrow b{\bar{b}}$ decay rate [20]. Details of these formalisms can be found in Refs.[37, 38] and [39, 40, 41] respectively.

While the full NLO and part of the NNLO QCD corrections to *VH * production are the same as those of the Drell–Yan (DY) process [45], with the Higgs boson radiated by the $V$ boson, additional contributions appear at NNLO, with the Higgs boson coupled to a heavy-quark loop. In the case of *ZH * production at the LHC, the impact of the gluon-gluon initiated subprocess involving a heavy-quark loop is substantial, due to the large gluon luminosity. We have taken into account these corrections with the full dependence on the top and bottom heavy-quark masses [17]. At NNLO there is yet another set of non DY like contributions involving quark induced heavy-quark loops both for *ZH * and *WH * production. These corrections have been computed in Ref. [11], relying in some cases on the large-$m_{t}$ approximation, and have been shown to have an impact on the *VH * cross section at the LHC at the $1\%$ level (for $m_{H}\sim 125$ GeV). However since the validity of the large-$m_{t}$ approximation is challenged in the high invariant mass region probed by the *VH * kinematics, we considered in our computation only the terms which can be presently calculated retaining the full $m_{t}$ dependence. In particular we included the NNLO terms obtained by radiating the Higgs boson off a top-quark bubble-insertion into an external gluon line. These terms, called $R_{I}$ in Ref. [11], contribute both to *ZH * and *WH * production. On the other hand the $R_{II}$ terms of Ref. [11], which are present only for *ZH * production, have been shown to contribute at the sub-per-mille level and have been thus neglected in this paper.

We are interested in the identification of the $b$-quark jets which originate from the Higgs boson decay. Besides the $b$-quark pair directly produced in the Higgs boson decay, we consistently include the effect of $b$-quark emissions from initial and final state partons §§§ Therefore, within our NNLO calculation, we have up to four $b$ quarks in the final state.. However the standard jet clustering algorithms [42] do not provide an infrared and collinear safe definition of flavoured jets with massless quarks. In the present case, at NNLO, the splitting of a gluon in a soft or collinear (massless) $b\bar{b}$ pair may affect the flavour of a jet. While the collinear unsafety can be removed by defining as a “$b$-jet” a jet containing a number of $b$ quarks different from the number of $\bar{b}$ quarks, the definition of infrared safe $b$-jets using standard jet clustering algorithms is less trivial. In order to deal with an infrared and collinear safe $b$-jet definition, we consider the so called flavour $k_{T}$ algorithm [43]. According to this algorithm, the definition of the $k_{T}$-distance measure in the presence of flavoured partons (particles) is modified in such a way that the flavour of a jet is insensitive to soft parton emissions.

We now present numerical results for $pp$ collisions at a center–of–mass energy of $\sqrt{s}=13$ TeV. For the electroweak couplings, we use the $G_{\mu}$ scheme and the following input parameters: $G_{F}=1.1663787\times 10^{-5}$ GeV*-2*, $m_{Z}=91.1876$ GeV, $m_{W}=80.385$ GeV, $\Gamma_{Z}=2.4952$ GeV, $\Gamma_{W}=2.085$ GeV, $m_{t}=172$ GeV and $m_{b}=4.18$ GeV ¶¶¶We consider the pole mass for the top quark ($m_{t}$) and the $\overline{MS}$ scheme for the bottom quark mass $m_{b}=\overline{m}_{b}(\overline{m}_{b})$.. The mass and the width of the SM Higgs boson are set to $m_{H}=125$ GeV and $\Gamma_{H}=4.070$ MeV respectively, while the $H\rightarrow b{\bar{b}}$ branching ratio is set to ${\mbox{Br}}(H\rightarrow b{\bar{b}})=0.578$ [36].

As for the parton distribution functions (PDFs), we use the NNLO PDF4LHC set [44] with $\alpha_{\mathrm{S}}(m_{Z})=0.118$. We set the renormalization and factorization scales to the dynamical value $\mu_{R}=\mu_{F}=M_{{\it VH\,}}$ (i.e. the invariant mass of the *VH * system) and the renormalization scale for the $H\rightarrow b{\bar{b}}$ coupling to the value $\mu_{r}=m_{H}$. To assess the impact of scale variation, we fix $\mu_{r}=m_{H}$ varying $\mu_{R}$ and $\mu_{F}$ independently in the range $M_{VH}/2\leq\{\mu_{R},\mu_{F}\}\leq 2M_{VH}$, with the constraint $1/2\leq\mu_{R}/\mu_{F}\leq 2$. We then fix $\mu_{R}=\mu_{F}=M_{VH}$ and vary the decay renormalisation scale $\mu_{r}$ between $m_{H}/2$ and $2m_{H}$. The final uncertainty is obtained by taking the envelope of the two (production and decay) scale uncertainties. Jets are reconstructed with the flavour-$k_{T}$ algorithm with $R=0.5$ [43]. We define a $b$-jet as a jet which contains a number of $b$ quarks different from the number of anti-$b$ quarks ($N(b)\neq N(\bar{b})$).

We start the presentation of our results by considering *W$\;\!{}^{+}\!$H * production and decay at the LHC at $\sqrt{s}=13$ TeV. Our choice of kinematical selection cuts on the final states closely follows the fiducial setup considered in the CERN Yellow Report of the LHC Higgs Cross Section Working Group [46]. We require the charged lepton to have transverse momentum $p_{T}^{l}>15$ GeV and pseudorapidity $|\eta_{l}|<2.5$ while the missing transverse energy of the event is required to be $E_{T}^{miss}>30$ GeV. The $W$ boson is required to have a transverse momentum $p_{T}^{W}>150$ GeV. Finally we require at least two $b$-jets each with $p_{T}^{b}>25$ GeV and $|\eta_{b}|<2.5$. The corresponding cross sections in the fiducial region are reported in the first row of Table 1, where we present the full NNLO prediction (see Eq. (S0.Ex2)) compared with the partial NNLO prediction (see Eq. (S0.Ex3)) ∥∥∥The results for the case of *W$\;\!{}^{-}\!$H * production and decay are qualitative similar, with a numerical reduction of fiducial cross section around $40\%$.. We observe that the inclusion of the full NNLO corrections reduces the cross section by around $6\%$ with respect to the partial NNLO result ******In particular we note that roughly $40\%$ of the reduction is due to the combination of the NLO contributions for production and decay and $60\%$ is due to the NNLO contributions to the decay rate (see Eq.(̇S0.Ex2) and subsequent comments)..

We next consider differential distributions. In Fig. 1 (left) we present the transverse-momentum distribution $p_{T}^{bb}$ of the leading $b$-jet pair (i.e. the two $b$-jets with largest $p_{T}$). In the lower panel we show the ratio of the two theoretical predictions defined above.

We observe that the additional $\alpha_{S}^{2}$ corrections included in the full NNLO prediction have an important effect also on the shape of the $p_{T}^{bb}$ distribution. In particular the cross section is increased by around $2-5\%$ for p_{T}^{bb}\raisebox{-2.58334pt}{\hbox to0.0pt{,\sim,\hss}}\raisebox{1.72218pt}{,<,}140 GeV and it is decreased by around $6-8\%$ for p_{T}^{bb}\raisebox{-2.58334pt}{\hbox to0.0pt{,\sim,\hss}}\raisebox{1.72218pt}{,>,}140 GeV. The corresponding $K$-factor, defined as the ratio between the full NNLO prediction in Eq. (S0.Ex2) and the partial NNLO prediction in Eq. (S0.Ex3), is thus remarkably not constant (see the lower panel of Fig. 1 (left)). The qualitative behaviour of these effects is not unexpected. The additional QCD radiation in the Higgs boson decay, which is included in the full NNLO calculation, has the effect of decreasing the transverse-momentum of the leading $b$-jet pair, making the $p_{T}^{bb}$ distribution softer.

In Fig. 1 (right) we present the invariant mass distribution of the leading $b$-jet pair, $M_{bb}$. We consider again the comparison between the full NNLO QCD prediction in Eq. (S0.Ex2) and the partial NNLO prediction in Eq. (S0.Ex3) and we show the ratio of the two predictions in the lower panel. For this observable the effect of the NNLO corrections to the decay rate are even more substantial. While the position of the peak is rather stable around the value of the Higgs boson mass $M_{bb}\simeq m_{H}$, the spectrum receives large positive corrections (up to $+60\%$) for $M_{bb}<m_{H}$ and sizeable negative corrections (from $-30\%$ to $-10\%$) for M_{bb}\raisebox{-2.58334pt}{\hbox to0.0pt{,\sim,\hss}}\raisebox{1.72218pt}{,>,}m_{H}. The large impact of these corrections can be understood by noting that the leading order (LO) computation would produce an invariant mass distribution which exactly fulfills the constraint $M_{bb}=m_{H}$. Higher-order corrections to the decay decrease the invariant mass of the leading $b$-jet pair. In the $M_{bb}<m_{H}$ region the partial NNLO prediction (which contains just the NLO correction to the decay rate) is effectively a first-order calculation and the next-order term is contained only in the *full * NNLO correction. Conversely, higher-order corrections to the production cross section typically increase the invariant mass of the leading $b$-jet pair and the region $M_{bb}>m_{H}$ receives contributions only from partons emitted from the initial state. In this case the effect of the additional $\alpha_{S}^{2}$ corrections contained in the *full * NNLO calculation has a sizeable but moderate impact with respect to the partial NNLO calculation.

As for the perturbative scale variation we have found that the scale dependence is dominated by the effect of the renormalization scale of the decay process $\mu_{r}$ and is particularly small: at the $1\%$ level for the fiducial cross section. The scale variation of the “full” NNLO result is around $\pm 5\%$ in the case of $p_{T}^{bb}$ distribtution and around $\pm 10\%$ in the case of $M_{bb}$ distribution. The scale dependence of the “partial” NNLO result is quantitatively similar being significantly larger (around $\pm 17\%$) only in the region $M_{bb}<m_{H}$, where the “partial” NNLO result is a first-order calculation.

We observe that the uncertainty bands for the “partial” and “full” NNLO results fail to overlap for the fiducial cross section and in various regions of differential distributions.

We next turn to the case of *ZH * production and decay at the LHC at $\sqrt{s}=13$ TeV. We consider the invisible $Z$ decay into neutrinos ($Z\rightarrow\nu\bar{\nu}$) and we require to have at least two $b$-jets each with $p_{T}^{b}>25$ GeV and $|\eta_{b}|<2.5$ and a missing transverse energy $E_{T}^{miss}>150$ GeV. The corresponding cross sections in the fiducial region are reported in the second row of Table 1. We observe that the inclusion of the full NNLO corrections reduces the cross section by around $5\%$ with respect to the partial NNLO result.

In Fig. 2 (left) we present the transverse-momentum distribution of the leading $b$-jet pair, $p_{T}^{bb}$. As in the previous case we compare the full NNLO QCD prediction (Eq. (S0.Ex2)) with the partial NNLO prediction (Eq. (S0.Ex3)) and in the lower panel we show the ratio of the two predictions.

In this case the inclusion of the NNLO corrections to the decay rate decreases the cross section up to about 10% below the peak and around 5% above the peak. The corresponding $K$-factor is shown in the lower panel of Fig. 2 (left).

Finally in Fig. 2 (right) we consider, for the *ZH * case, the invariant mass distribution of the leading $b$-jet pair, $M_{bb}$. The effect of the full NNLO corrections is similar to the *W$\;\!{}^{+}\!$H * case. The spectrum receives large positive corrections (up to $+70\%$) for $M_{bb}<m_{H}$ and sizeable negative corrections (from $-30\%$ to $-10\%$) for 125\raisebox{-2.58334pt}{\hbox to0.0pt{,\sim,\hss}}\raisebox{1.72218pt}{,<,}M_{bb}\raisebox{-2.58334pt}{\hbox to0.0pt{,\sim,\hss}}\raisebox{1.72218pt}{,<,}150 GeV.

We finally observe that when a kinematical boundary is present at a given order in perturbation theory, higher order corrections are affected by instabilities of Sudakov type [47], which spoil the reliability of the fixed-order expansion around the boundary. This is the case for both the $p_{T}^{bb}$ distribution ($p_{T}^{V}>150$ GeV LO kinematical boundary) and the $M_{bb}$ distribution (LO condition $M_{bb}=m_{H}$). While a proper treatment of this misbehaviour requires an all order resummation of perturbatively enhanced terms, the effect of these instabilities can be mitigated by increasing the bin size of the distribution around the critical point.

In the case of $ZH$ production, due to the substantial effect of the gluon-gluon initiated subprocess involving a heavy-quark loop, scale uncertainty is dominated by the effect of the renormalization scale $\mu_{R}$ and the ensuing scale variation band turns out to be larger (at the $4\%$ level).

The scale variation of the “full” NNLO result is around $\pm 3-5\%$ in the case of $p_{T}^{bb}$ distribtution and around $\pm 10\%$ in the case of $M_{bb}$ distribution. As in the *W$\;\!{}^{+}\!$H * case, the scale dependence of the “partial” NNLO result is significantly larger (around $\pm 17\%$) only in the $M_{bb}<m_{H}$ region and the uncertainty bands for the “partial” and “full” NNLO results fail to overlap in various regions of differential distributions.

As already pointed out in Ref. [17], we are interested in a specific “boosted” kinematical regime where the size of the NNLO corrections tends to be underestimated by the customary NLO scale uncertainty band. Therefore the NLO scale variation cannot be regarded as a reliable approximation of the “true” perturbative uncertainty and it casts some doubts also on the reliability of the NNLO scale variation band.

A hint on the reliability of the customary scale uncertainty at NNLO can be obtained considering missing higher-order contributions that can be calculated through a suitable combination of individual parts of our computation. We have therefore calculated the $\mathcal{O}(\alpha_{S}^{3})$ contributions proportional to (i) $d\sigma^{(1)}_{h_{1}h_{2}\rightarrow{\it VH\,}}\times d\Gamma^{(2)}_{H\rightarrow b\bar{b}}$ and (ii) $d\sigma^{(2)}_{h_{1}h_{2}\rightarrow{\it VH\,}}\times d\Gamma^{(1)}_{H\rightarrow b\bar{b}}$ in Eq. 1. We have found that the numerical impact to the fiducial cross-section of the N3LO terms (i) and (ii) above is respectively around $-0.4\%$ and $+0.4\%$ ($-0.3\%$ and $+1.5\%$) for $WH$ ($ZH$) production. The fact that these effects are covered by the scale variation in Tab. 1 suggests that the NNLO scale dependence (contrary to the NLO case) could be considered as a trustable estimate of the “true” perturbative uncertainty of the calculation. However a more conservative estimate of the uncertainty can be obtained by comparing the NNLO result to what is obtained at the previous order.

We briefly comment on expected PDF uncertainties and on electroweak effects. PDF uncertainty has been calculated, within a similar setup, in Ref. [46] and has been shown to be at the $\pm 1.5\%$ level for fiducial cross sections. The NLO EW effects have been calculated for $pp\rightarrow VH+X\rightarrow l_{1}l_{2}H+X$ only (i.e. without the inclusion of EW effects for $H\rightarrow b\bar{b}$ decay) [14, 34] and it has been shown to be significant ($\sim-10\%$) [46].

In conclusion, we have presented a fully differential QCD computation for the associated production of a vector boson and a Standard Model Higgs boson in hadron collisions including the QCD radiative corrections up to next-to-next-to leading order (NNLO) both for the *VH * production cross section and for the differential Higgs boson decay width into bottom quarks. Our calculation also includes the leptonic decay of the vector boson with finite-width effects and spin correlations and it is implemented in the parton level Monte Carlo numerical code HVNNLO.

We have studied the impact of the full NNLO QCD corrections to the *VH * production and decay at the LHC by focusing on the most relevant distributions, namely the transverse momentum and invariant mass of the Higgs boson candidate. We have studied the renormalization and factorization scale dependence of the results in order to estimate the perturbative uncertainty of our predictions. We have found that the additional second-order corrections included in the present calculation have a substantial effect both for *WH * and *ZH * production. Therefore the inclusion of these effects turns out to be essential in order to obtain a precise theoretical prediction for associated *VH * production and decay at the LHC.

Acknowledgements.

We gratefully acknowledge Massimiliano Grazzini for discussions and comments on the manuscript, Raoul Röntsch for comparison with the results of Ref. [48], and the computing resources provided by the SCOPE Facility at the University of Napoli and by the Theory Farm at the University of Milano and INFN Sezione di Milano. This research was supported in part by Fondazione Cariplo under the grant number 2015-0761, by the Italian Ministry of Education and Research (MIUR) under project number 2015P5SBHT and by grant K 125105 of the National Research, Development and Innovation Fund in Hungary.