Bottom quark mass effects in associated $WH$ production with $H \to b\bar{b}$ decay through NNLO QCD

We present a computation of NNLO QCD corrections to the production of a Higgs boson in association with a $W$ boson at the LHC followed by the decay of the Higgs boson to a $b\bar{b}$ pair. At variance with previous NNLO QCD studies of the same process, we treat $b$ quarks as massive. An important advantage of working with massive $b$ quarks is that it makes the use of flavor jet algorithms unnecessary and allows us to employ conventional jet algorithms to define $b$ jets. We compare NNLO QCD descriptions of the associated $WH(b\bar{b})$ production with massive and massless $b$ quarks and also contrast them with the results provided by parton showers. We find ${\cal O}(5\%)$ differences in fiducial cross sections computed with massless and massive $b$ quarks. We also observe that much larger differences between massless and massive results, as well as between fixed-order and parton-shower results, can arise in selected kinematic distributions.


I. INTRODUCTION
Detailed investigation of Higgs boson production in association with a W boson is an important part of the LHC research program that aims at a comprehensive exploration of Higgs boson properties and electroweak symmetry breaking [1][2][3][4][5]. Indeed, associated Higgs boson production gives us direct access to the HW W coupling which is completely fixed in the Standard Model but can be modified in its extensions. Moreover, studies of the pp → W H process provide a unique way to study the Higgs coupling to b quarks since, by selecting Higgs bosons with relatively high transverse momenta, one can identify H → bb decays using substructure techniques [6,7].
Interest in associated W H production has inspired a large number of theoretical computations that provide refined descriptions of this process including QCD [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22] and electroweak radiative corrections [23,24]. The more recent theoretical efforts [17][18][19][20][21][22] focused on a comprehensive fully-differential description of associated production which consistently combines QCD corrections to the production and decay processes. All fully-differential NNLO QCD computations mentioned above have the common feature that the decay of the Higgs boson to b quarks is described under the assumption that b quarks are massless. The same approximation is employed in the production subprocesses which involve gluon splitting into a bb pair or b quarks that come directly from the proton.
Although, given the high energy of the LHC, the massless approximation should be fairly adequate, there are a few reasons that make it interesting to explore b-quark mass effects in this process.
The first reason is that the phase space of the pp → W H process is large and that there are important kinematic distributions which can be sensitive to energy scales smaller than the total (partonic) energy of the process. In those cases the dependence on the b-quark mass can become more pronounced. A comparison of fixed-order computations, including higher-order ones, for pp → W H performed with massless and massive b quarks will allow us to identify distributions and phase-space regions with enhanced sensitivity to the b-quark mass.
The second reason to employ massive b quarks in a calculation is that in this case the splitting g * → bb becomes non-singular. This feature makes it possible to employ conventional jet algorithms to define b jets. We remind the reader that in case of massless b quarks, this is not possible and that a special partonic flavor jet algorithm [25] has to be used. The possibility to apply conventional jet algorithms is an important improvement since it makes theoretical computations and experimental analyses more aligned.
The third reason is the appearance of certain contributions in the H → bb decay which cannot be properly described if b quarks are treated as massless. It was pointed out in Ref. [21] that an interference of singlet H → g(g * → bb) and non-singlet H → (b * → bg)b decay amplitudes forces us to keep the mass of the b quark different from zero throughout the computation since otherwise unconventional soft quark divergences appear. Such studies have already been carried out in Ref. [26].
Motivated by these considerations, we extended the computation of NNLO QCD radiative corrections reported in Ref. [21] to include b-quark mass effects in the theoretical description of Higgs production in association with a vector boson. To this end, we combined the recent NNLO QCD description of the Higgs boson decay into two massive b quarks [27] 1 with the computation of NNLO QCD corrections to the production process [21,29] which required small modifications because of the b-quark mass.
In addition to fixed-order computations, parton showers are widely used to provide theoretical predictions for collider experiments. In the context of associated Higgs production, they have been employed in Refs. [30][31][32][33][34][35][36]. For this reason, it is interesting to compare fixedorder and parton-shower results with each other. Although this has already been done in Ref. [21], the need to use different jet algorithms in fixed-order massless and parton-shower computations did not allow a direct comparison of the two. The NNLO QCD computation with massive b quarks described in this paper allows us to remedy this problem and compare fixed-order and parton-shower predictions using identical jet algorithms.
The remainder of the paper is organized as follows. In Section II we briefly review the NNLO QCD computation of radiative corrections to pp → W H [21] and H → bb [27] and discuss modifications needed in the computation of NNLO QCD corrections to the production process to accommodate massive b quarks. In Section III we show numerical results for NNLO QCD corrections to pp → W H(bb) with massive b quarks and compare 1 A calculation of NNLO QCD corrections to the H → bb decay with massive b quarks was also performed in Ref. [28]. them to results of the massless computation. In Section IV we compare a parton-shower description of associated W H production with fixed-order results. We conclude in Section V.
A detailed discussion of modifications required to accommodate massive b quarks in the NNLO QCD computation of Ref. [21] can be found in two appendices.

II. SUMMARY OF NNLO QCD COMPUTATIONS
In this section, we briefly review the computation of NNLO QCD radiative corrections to the associated production pp → W H and the H → bb decay processes. An earlier computation of NNLO QCD corrections to pp → W H was described in Ref. [21] using the formulation of the nested soft-collinear subtraction scheme presented in Ref. [37]. Since then, simple analytic formulas for the NNLO QCD corrections to the production of a color-singlet final state in hadron collisions were published in Ref. [29]. These formulas employ results for integrated double-unresolved soft and collinear subtraction terms computed in Refs. [38] and [39], respectively. To accommodate these developments, the code that allows us to compute NNLO QCD corrections to pp → W H was updated. In addition, we refined the description of the H → bb decay with massless b quarks using analytic results for NNLO QCD corrections to decays of color-singlet particles derived in Ref. [40].
To accommodate massive b quarks, we employed a recent computation [27] of the NNLO QCD corrections to H → bb that fully accounts for the b-quark mass. That computation is based on the nested soft-collinear subtraction scheme adapted to deal with massive particles.
On the production side, a consistent description of b quarks as massive particles forces us to work in a four-flavor scheme so that b quarks are excluded from parton distributions.
This feature leads to some changes to the renormalization procedure that we discuss in Appendix A. In addition, we have to modify the computation of NNLO QCD corrections to pp → W H to describe the splitting of a gluon into a massive bb pair, and the gluon vacuum polarization contributions due to massive b-quark loop.
We note that the gluon splitting contribution refers to the process q i q j → W H + (g * → bb) which is free of soft and collinear singularities thanks to the finite mass of the b quark.
The resulting logarithmically enhanced terms of the form log(s/m 2 b ) may, potentially, be large, but they do not appear to be particularly problematic from a numerical viewpoint.
Hence, to describe these contributions, we calculate helicity amplitudes for the q i q j → W Hbb process, parametrize the W Hbb phase space and perform numerical integration to compute observables of our choice.
Two-loop corrections to the q iqj W vertex caused by the gluon vacuum polarization due to a massive quark loop can be extracted from Refs. [41][42][43]. We recomputed these contributions and found full agreement with the results presented in Ref. [41]. For completeness, we provide the details of our calculation in Appendix B.

III. THE PROCESS pp → W H(bb)
In this section we present results for the associated production pp → W H(bb) including b-quark mass effects. We begin by specifying how corrections to production and decay processes are combined. Since the Higgs boson is a scalar particle, these corrections can be put together in a straightforward manner. The only subtlety worth discussing is how to treat the total decay width of the Higgs boson that appears in the differential cross section for pp → W H(bb) when it is computed in the narrow-width approximation. We begin by writing the cross section as follows [18] dσ W H(bb We treat Br(H → bb) as an input parameter and do not expand it in a series in α s 2 . For numerical computations we take Br(H → bb) = 0.5824, as recommended by the Higgs Cross Section Working Group [44].
Keeping the branching fraction fixed, we compute an expansion of Eq. (1) in α s by first expanding the W H cross section and the decay rate then introducing normalized quantities to describe the decays and, finally, defining physical cross sections computed through different orders in QCD perturbation theory We note that dγ (i) = 1 provided that the integration is performed over unrestricted phase space. An identical definition of the cross section was used in an earlier massless computation reported in Ref. [21].
To present the results of our computation, we focus on the associated production process We treat both decay processes W + → ν eē and H → bb in the narrow-width approximation.
We we approximate the CKM matrix by an identity matrix. 4 We also need to specify the selection criteria that are used to define the W (ν eē ) H(bb) final state. To this end, we require that an event contains at least two b jets that are defined with the anti-k t jet algorithm [45,46]. For the sake of comparison, we also calculate W H(bb) cross sections for massless b quarks. In that case, we employ the flavor-k t jet algorithm [25] to describe massless b jets. In both cases, we choose the jet radius R = 0.4. Moreover, we 3 At NNLO a residual dependence on y b remains in the dΓ bb /Γ bb ratio because of the singlet-non-singlet interference which depends on the product of top and bottom Yukawa couplings. 4 We have checked through NLO QCD that in case of the associated production, this approximation is accurate to about a percent. impose the following cuts on pseudo-rapidities and transverse momenta of leptons and b jets |η l | < 2.5 , p t,l > 15 GeV , Finally, following experimental analyses, we may employ an additional requirement that the vector boson has a transverse momentum of p t,W > 150 GeV. We always state explicitly when this cut is used.
To present numerical results we use the five-and four-flavor parton distribution function sets (p W + p H ) 2 , whereas the renormalization scale for the decay process is set to the Higgs boson mass, µ r,dec = M H . The uncertainty of the cross sections is obtained by varying the scale in the production process by a factor of two and, independently, by changing the decay scale by a factor of two as well. The total uncertainty is taken to be an envelope of these nine numbers.
We begin by presenting fiducial cross sections for the process pp → W + H(bb) at the 13 TeV LHC in Table I. Comparing these results with massless predictions, we observe that the massive cross sections are systematically larger than the massless ones. The difference is very small at LO but increases when radiative corrections are included. At NLO, the differences range from about four percent, in case of the basic fiducial cuts, to six percent, if the additional p t,W > 150 GeV cut is applied. At NNLO, the differences between massive and massless results increase further and reach 6 − 7 percent.
We note that these differences may be obscured by the scale variation uncertainties. This is indeed what happens at leading and, to some extent, also at next-to-leading order, whereas  We emphasize that the NNLO scale variation uncertainties shown in Table I are likely to be too conservative [22]. Indeed, it was shown in Ref. The O(5%) differences between massive and massless fiducial cross sections can be traced back to gluon radiation in H → bb decays. Indeed, it is well-known that the collinear radiation pattern of massive and massless b quarks differs significantly. In case of massless b quarks, we expect a logarithmic enhancement of the collinear gluon emission probability dP ∼ dθ 2 /θ 2 , where θ is the relative angle between the b quark and the gluon momenta. This feature leads to a logarithmic dependence of the fiducial cross section on the clustering radius R. At the same time, when massive b quarks radiate, the probability distribution becomes , where E b is the energy of the radiating quark. This probability distribution implies that the collinear singularity is screened by the b-quark mass and that the cross-section dependence on the jet radius changes if ∆R < m b /p t,b . We have checked that for the chosen value of the jet radius, the amount of radiation included in a b jet is different for massless and massive quarks. This means that, in the case of radiative decay of the Higgs boson H → bbg, the acceptance of events with massless b quarks is smaller, by about a factor of two, than the acceptance computed with massive b quarks, when fiducial cuts described above are applied. We also observe that this difference is reduced if we consider larger jet radii or reduce the transverse momentum cut on b jets.
Finally, it turns out that the O(y t y b ) interference of singlet H → g(g * → bb) and non- , is a minor effect. For the fiducial cuts discussed above, it contributes to cross sections only at a sub-percent level and is, therefore, below the scale uncertainty and much smaller than the differences between massive and massless computations.
We will now proceed with the discussion of kinematic distributions. Since in an experimental analysis the Higgs boson can only be observed through its decay products, we will study kinematic distributions of the bandb-jet pair whose invariant mass is closest to the Higgs boson mass. Throughout this paper, we refer to such pairs of jets with the subscript H(bb), e.g. their four-momenta are written as p H(bb) and their invariant masses as M H(bb) .
We begin by presenting the rapidity distribution of pairs of b jets in Fig. 1. We observe that the distributions computed with massive and massless b quarks are very similar and differ, to a good approximation, by an overall rescaling factor that can be inferred from the results for the cross sections reported in Table I. Such behavior is expected given the well-known inclusiveness of rapidity distributions.
We proceed with the invariant mass distribution of the two b jets, M H(bb) , which is presented in Fig. 2. At leading order this distribution is described by a δ-function, δ(M 2 H(bb) − M 2 H ), but the situation becomes more complex when higher-order corrections are considered. In particular, if a b quark from the decay is clustered with a gluon emitted in the production process, the invariant mass of two b jets can exceed M H and, conversely, a three-body decay H → bbg leads to two b jets with an invariant mass that is smaller than M H . Hence, already at NLO, the M H(bb) distribution is non-vanishing both below and above M H . We present the M H(bb) distributions obtained at NLO and NNLO in Fig. 2  is not applied, we observe that below the Higgs peak, the massless results are larger than the massive ones except at very low invariant masses. In the region above the peak, which is populated by events with radiative corrections to the production process, the two results are very similar to each other. In the most populated bin, adjacent to the Higgs boson mass, M H(bb) = M H , the massive result is larger than the massless one; this feature drives the observed behavior for fiducial cross sections discussed earlier, c.f. Table I. When the additional p t,W > 150 GeV cut is applied, the massless result stays below the massive one; we observe an O(15%) difference at very low invariant masses which decreases when getting closer to the peak. Above the Higgs mass, we see a constant difference of about 10%.
Next, we consider the transverse momentum distribution of those b-jet pairs whose invariant Such events are then rejected by fiducial cuts since (at least) two b jets are required. Since As we already pointed out, some differences in kinematic distributions computed with massive and massless quarks arise already at leading order. If we assume that radiative effects are similar in massive and massless cases, one can construct approximate NNLO distributions from massive NLO computations and massless differential K-factors defined as dK = dσ NNLO /dσ NLO . We compare the (so constructed) approximate and exact NNLO distributions for M H(bb) and p t,H(bb) in Fig. 6. We observe that such an approximation is only partially successful; it provides a decent description of the true p t,H(bb) distribution but does not capture all the details of the M H(bb) spectrum.

IV. COMPARISON WITH PARTON SHOWER
Having discussed fixed-order calculations with massive and massless b quarks, we turn to a comparison of these calculations with parton showers. Such a comparison is important because experimental analyses often rely on parton showers and one needs to understand their reliability by comparing them to fixed-order computations.
The second result shown in Eq. (7) is obtained by requiring that, in addition to standard fiducial cuts, the transverse momentum of the W boson p t,W exceeds 150 GeV. The uncertainties shown in Eq. (7) correspond to numerical integration errors.
The parton-shower cross sections Eq. (7) differ from NNLO cross sections computed with massive b quarks by about 2% in the full fiducial region and by about 4% if the additional p t,W cut is applied (c.f. Table I). These differences are only natural given that the POWHEG+Pythia8+MiNLO setup is different compared to what we use to obtain fixed-order predictions, see Ref. [33] for further details.
We proceed with the comparison of fixed-order and the parton-shower descriptions of selected kinematic distributions for a pair of b jets whose invariant mass is closest to the mass of the Higgs boson. We present the transverse momentum distribution of such b-jet pairs in Fig. 7, and their invariant mass distribution in Fig. 8. In the case of the transverse momentum distribution, both with and without the additional p t,W cut, we see that in the region p t,H(bb) 100 GeV the parton-shower result is smaller than the massive NNLO result by about five percent, whereas for transverse momenta below the peak of the distribution, p t,H(bb) 50 GeV, the parton-shower prediction exceeds the fixed-order result by about five percent. We note that such behavior is expected since additional QCD radiation, simulated by a parton shower, reduces energies of the b jets leading to a softer spectrum.
Differences between parton-shower predictions and the massive fixed-order NNLO result for the invariant mass of the bb-system are more significant than in case of the transverse  are O(25%) smaller than fixed-order results. We note that the parton-shower and the fixedorder distributions can be made well aligned provided that the fixed-order distribution is shifted along the x axis by δM bb ∼ −4 GeV.

V. CONCLUSIONS
In this paper, we discussed the associated production of the Higgs boson, pp → W H, and the decay of the Higgs boson to bb pairs at the LHC. We included the NNLO QCD corrections to the production and decay processes, retaining the dependence on the b-quark mass. The inclusion of the b-quark mass in the calculation is important as it allows us to use realistic jet algorithms to describe b jets, making theoretical and experimental analyses more aligned.
We compared theoretical predictions for the associated production that are obtained with massive and massless b quarks. We observed O(6%) differences between the two results once fiducial cuts are applied. Such relatively large differences can be traced back to different acceptances in radiative decays of the Higgs boson H → bbg when they are computed in the massive and in the massless approximations for a standard set of fiducial volume cuts. We also found that radiative corrections to the production process are less sensitive to b-quark mass effects.
Interestingly, mass effects can become much more pronounced in kinematic distributions.
For example, we observed large differences between massive and massless predictions in kinematic regions where b jets have large transverse momenta. In these cases, differences in clustering algorithms employed with massive and massless partons, needed to unambiguously define a jet's flavor, combine with rapidly changing distributions and lead to O(20%) discrepancies between the theoretical predictions.
We note that in some cases such large discrepancies are driven by differences in lower-order distributions while massive and massless K-factors turn out to be similar. If this is the case, an approximate massive NNLO result may be constructed from massive NLO result and massless NNLO/NLO K-factor. We have identified the transverse momentum p t,H(bb) as one such observable. However, there are also other cases where the differences in NNLO distributions are driven by different (massive and massless) K-factors; if this is the case, the approximate distribution will not provide a decent description of the true result. This is the case, e.g., for the invariant mass M H(bb) .
Differences between massive NNLO QCD and parton-shower computations, discussed in Sec. IV, are easily understood if we assume that parent b quarks lose more energy in a parton-shower computation than in a fixed-order one. This implies that shapes of, at least some, distributions in both cases are similar but the distributions themselves are shifted relative to each other, e.g. dσ (P S) /dx (x) ∼ dσ (F O) /dx (x + δ x ). We have found that, in case of the invariant mass of two hardest b jets, δ x ∼ 4 GeV which appears to be a rather natural value.
In summary, we studied effects of the b-quark mass on associated production of the Higgs boson, pp → W H, followed by decay of the Higgs boson into a bb pair. Although such effects are not large, we found that they are typically larger than naively expected and that they can affect both fiducial cross sections and kinematic distributions in a somewhat unexpected way. We look forward to future studies of such effects in other processes relevant for the LHC phenomenology. The renormalization of the H → bb decay process was discussed at length in Ref. [27] and we do not repeat it here. Instead, in this appendix, we focus on the production process.
We start by discussing the renormalization of the qq → V H amplitude A with q being a massless quark, i.e. q = b. Neglecting b-quark contributions altogether and considering n f = 4 massless flavors, we write the MS-renormalized amplitude as We continue by expressing Eq. (A1) through the bare coupling constant α s,b and find where S = (4π) e − γ E is the standard MS factor and β (n f ) 0 In order to include the b-quark contribution to Eq. (A1), we need to add the gluon vacuum polarization diagram Fig. 9a and to account for additional contributions to renormalization constants that arise in the theory due to loops with massive b quarks. For the amplitude A an additional renormalization factor is the wave function renormalization constant of a massless quark Z q that receives b-quark contributions at two loops, see Fig. 9b. Another contribution that arises in the theory with massive b quarks is the gluon wave function renormalization constant Z A , see Fig. 9c.
Starting from Eq. (A2), we re-express the renormalized amplitude through the coupling constant defined in a theory with five active flavors. We find (A4) From now on, we will always work with α s renormalized in a theory with n f = 5 massless flavors at a scale µ. Therefore, unless stated otherwise, we will use the short-hand notation To proceed further, it is convenient to express A (n f =5) through two-loop contributions to the wave function renormalization constants Z q and Z A . To this end, we write We leave the discussion of the massless quark and gluon self-energies,Σ 2 (0) and Π 1 (0), to Appendix B. Here, we only remark that the difference of the two β-functions in Eq. (A4) can be expressed through Π 1 (0) and an additional constant term, cf. Eq. (B16). Hence, we with Using Eqs. (A6) and (A5) we write Eq. (A4) as where we introduced The square of the amplitude A (n f =5) expanded to second order in α s gives the following contribution to the cross section where dσ V and dσ VV are the one-and the two-loop contributions to cross sections, respectively, and dσ VV,(b,reg) is the two-loop contribution proportional to 2Re in Appendix B.
We now discuss the amplitude qq → V H +g which is needed to describe real and real-virtual contributions to NLO and NNLO cross sections. As in the previous case, we start with the amplitude computed in a theory with n f = 4 massless quarks and write In Eq. (A10) g s . Equivalently, we re-express Eq. (A10) using the bare coupling constant In this case, there are no explicit n f -dependent contributions to the unrenormalized amplitude so that all the b-quark effects only enter through the renormalization. Since Z q = 1 + O(α 2 s ), we only need to renormalize the strong coupling constant α s and to multiply the unrenormalized amplitude by the gluon renormalization factor √ Z A . We obtain where g s = g (5) s (µ) is the strong coupling constant defined in the theory with five flavors and renormalized at a scale µ. We finally write the contribution of the renormalized qq → V H +g amplitude Eq. (A12) to the cross section The last two contributions that we need to discuss are the double-real emission processes and the PDFs renormalization term. The double-real emission processes do not require any renormalization and can be obtained as a direct sum of n f = 4 contributions that we discussed earlier [29,37] and an additional finite contribution where a virtual gluon splits into a massive bb pair.
In the context of PDF renormalization, we stress that we work in a theory with four active massless flavors in the proton, but we write the result using the QCD coupling constant computed in a theory with n f = 5 flavors. Taking into account the change in the coupling constant, we find an additional contribution to the NNLO cross section that reads whereP (0) are the LO Altarelli-Parisi splitting functions and "⊗" denotes the standard convolution product, see Ref. [37] for more details.
The resulting NNLO cross section is obtained by combining Eqs. (A9), (A13) and (A15). We find that the terms proportional to K 1 assemble themselves into a finite NLO cross section.
Therefore, we write Finally, we emphasize that no modifications are required to compute leading and next-toleading order W H production cross sections.
Appendix B: Contributions of a massive b quark to a two-loop form factor of a massless quark In this appendix, we calculate the contribution of a massive b quark to the two-loop amplitude A (b,reg) 2 defined in Eq. (A8). We note that such a calculation was performed in Refs. [41][42][43]; we discuss it here for completeness.
We begin by considering A (b,bare) 2 , which corresponds to Fig. 9a. Since helicity of a massless quark is conserved and since flavor-changing currents are anomaly-free, there is no difference between the form factors of a vector and of a vector-axial current. Therefore, for simplicity we consider radiative corrections to a matrix element of a generic vector current J µ =qγ µ q between the vacuum state and a qq pair q(p 1 )q(p 2 )|J µ (0)|0 .
The gluon self-energy Π(k 2 ) satisfies the once-subtracted dispersion relation We now insert this dispersion relation into Eq. (B1). The Π(0) term gives rise to a contribution proportional to the one loop amplitude A 1 in the Landau gauge. However, since A 1 is gauge-independent, we can write The term in the second line of Eq. (B5) is proportional to the one-loop vertex correction due to an exchange of a gluon with the mass q 2 in the Landau gauge. As a consequence, it is both UV and IR finite. After simple manipulations, we cast Eq. (B5) into the following form iA (b,bare) 2 dq 2 q 2 ImΠ(q 2 ) Tri(d, q 2 , s).

(B6)
We note that in the limit q 2 → ∞, both Π(q 2 ) and Tri(d, q 2 , s) approach constants, so that the integration over q 2 diverges. To remove this divergence, we need to incorporate the wave function renormalization constant of a light quark, (Z q − 1) ∼Σ 2 (0), cf Eq. (A5), into the computation.
To computeΣ 2 (0), we evaluate the self-energy in Fig. 9b and write We note that, thanks to helicity conservation, the self-energyΣ is proportional top We extractΣ from Eq. (B7) and use dispersion relations, Eq. (B4), to arrive at Combining Eq. (B9) with Eq. (B6), we find that which implies that in a combination of the relevant vertex correction and the wave-function renormalization contribution the constant asymptotic at large q 2 cancels out and the q 2 integration becomes convergent. This allows us to take the d → 4 limit in Π(q 2 ). Following this discussion, we write the regulated b-quark amplitude in Eq. (A8) as dq 2 q 2 ImΠ(q 2 ) Tri(d = 4, q 2 , s) + Inserting Eq. (B12) into Eq. (B11) and integrating over q 2 , we obtain the final result