Higgs boson decay into four bottom quarks in the SM and beyond

We present predictions for the Higgs boson decay into four bottom quarks in the standard model and via light exotic scalars retaining full bottom-quark mass dependence. In the SM the decay can be induced either by the Yukawa couplings of bottom quarks and top quarks or the electroweak couplings. We calculate the partial decay width and various differential distributions up to next-to-leading order in QCD. We find large QCD corrections for decay via Yukawa couplings, as large as 90% for the partial decay width, and reduced scale variations. The results of this paper are therefore helpful for the measurement of this multi-jets final state at future Higgs factory of electron-positron colliders. We also propose several observables that can differentiate the SM decay channel and the exotic decay channel and compare their next-to-leading order predictions.


Introduction
The successful operation of the LHC and the ATLAS and CMS experiments have led to the discovery of the Higgs boson and completion of the standard model (SM) of particle physics [1,2]. Further refined study at the LHC has revealed one essential of the Higgs boson, Yukawa couplings of top quarks [3,4] and bottom quarks [5,6]. Precision test on properties of the Higgs boson including all its couplings with standard model particles becomes one primary task of particle physics at the high energy frontier. In the SM the Higgs boson decays dominantly to hadronic final states which are hard to access at hadron colliders due to huge QCD backgrounds. That includes decay channels of a pair of bottom quarks or charm quarks, a pair of gluons via top-quark loops, and four quarks via electroweak gauge bosons, adding to a total decay branching ratio of about 80% [7]. To measure the Higgs properties with higher accuracy and to probe rare decay modes of the Higgs boson, there have been a few proposals to build a future lepton collider that can serve as a Higgs factory. These include ILC [8], CEPC [9], CLIC [10] and FCC-ee [11]. The proposed LHeC program can also produce abundant Higgs boson and measure its couplings precisely [12,13]. Certainly with the future Higgs factory, the precious hadronic decay channles of the Higgs boson can be studied in details thanks to the clean environment [14].
Precision experiments require equally precision theoretical predictions. To further scrutinize the SM and to look for possible new physics beyond, it is necessary to calculate higher-order corrections to the production and decay of the Higgs boson. In this respect, there have been enormous advances in recent years. For example, the next-to-next-to-nextto-leading order (N 3 LO) QCD corrections to Higgs boson production via gluon fusion in the heavy top-quark limit [15,16] and to Higgs boson production via vector boson fusion within the structure function approach [17], the next-to-next-to-leading order (NNLO) corrections to Higgs boson production in association with a jet in the heavy top-quark limit [18][19][20][21], and the next-to-leading order (NLO) corrections to Higgs boson pair production with full top-quark mass dependence [22] have been known for some time. The two-loop mixed QCD and electroweak corrections have also been calculated recently for the associated production of Higgs boson and a Z boson at electron-positron colliders [23][24][25]. For decays of the Higgs boson, the partial width for H → bb is known up to the next-to-next-to-next-to-next-toleading order (N 4 LO), in the limit where the mass of the bottom quark is neglected [26]. The partial width for H → gg has been calculated to the N 3 LO in the heavy top-quark limit [27]. We refer the readers to [28,29] for a complete list of relevant calculations. At a more exclusive level, the fully differential cross sections for H → bb have been calculated to NNLO in [30,31] and N 3 LO in [32] for massless bottom quarks, and to NNLO in [33] with massive bottom quarks.
Besides, multi-jets final state, for instance, Higgs boson decays into three or four QCD partons can also be explored at future Higgs factory. There have been strong motivations for experimental searches of those exclusive hadronic decay channels to look for exotic decays of the Higgs boson induced by light states beyond the SM [34][35][36]. Among them there is the decay of the Higgs boson to four bottom quarks which will be the main focus of this paper. There have already been searches carried out by ATLAS collaboration at the LHC [37] setting an upper limit of about 50% on the decay branching ratio to four bottom quarks. On the theory side, the Higgs boson decays into four massless quarks via electroweak gauge bosons have been calculated to NLO in both QCD and electroweak couplings and have been implemented in the MC program Prophecy4f [38,39]. There have also been recent predictions for Higgs boson decays into three-jets final state [40][41][42] for the hadronic event shapes. That is particularly interested for the probe of light-quark Yukawa interactions [43]. Very recently there exist a NNLO calculation for Higgs boson decays into a pair of bottom quarks plus an additional jet for massless bottom quarks [44].
In this paper, we present predictions for the Higgs boson decays into four bottom quarks in the standard model and via light exotic scalars retaining full bottom-quark mass dependence. In the SM the decays can be induced either by the Yukawa couplings of bottom quarks and top quarks or the electroweak couplings. We calculate the partial decay width and various differential distributions up to next-to-leading order in QCD. We discuss details of the calculation in Section 2 and present numerical results for the SM case. We then move to discussion on exotic decays induced by new scalars including for the QCD corrections and comparisons to the SM case in Section 3. We conclude in Section 4.

Decay via Yukawa interactions
In the full theory with n f = 6 active flavor of quarks, the relevant interactions include the Yukawa couplings of quarks and the QCD interactions, We neglect quark masses except for top quark and bottom quark. We use onshell scheme for renormalization of the quark or gluon fields and masses except for masses in Yukawa couplings, when calculating the QCD radiative corrections. The renormalization of the QCD coupling is carried out in MS scheme with top quark decoupled [45], namely the number of light and heavy flavor, n l = 5 and n h = 1 respectively. The corresponding renormalization constants at one-loop are, with S = (4π) /Γ(1 − ), and β (n l ) 0 = (11C A − 4T R n l )/6. Furthermore, since the relevant hard scale is much larger than the bottom quark mass, we use the MS running mass of bottom quark in Yukawa coupling, that can be related to the input pole mass [46].
Since the top quarks only appear as internal states, one can adopt an effective theory by integrating out top quark, the interactions can be expressed as to the leading power of inverse of the top-quark mass. The Wilson coefficients C 1 and C 2 [47][48][49][50][51][52][53] carry further logarithmic dependence on top-quark mass and are given by The renormalization of QCD coupling and gluon field are done with pure MS scheme with n l = 5 light flavors. Again we use onshell scheme for renormalization of bottom quark field and mass, apart from using MS running mass in the Yukawa coupling. That is equivalent to set n h = 0 in Eq. (2.2). Moreover, in the effective theory, there will be operator mixing between the last two terms in Eq. (2.3) requiring renormalization of the Wilson coefficients, with the one-loop renormalization constants in MS scheme given by [53], Under the effective theory the leading order (LO) Feynman diagrams for the Higgs boson decaying into four bottom quarks are shown in Fig. 1. Those diagrams can be obtained with interchanges of identical particles are not shown for simplicity. In this calculation the squared amplitudes needed for a next-to-leading order QCD calculation are generated automatically with program GoSam 2.0 [54], including for the one-loop virtual corrections and the real corrections. Reduction of loop integrals are performed with Ninja [55,56] and scalar integrals are calculated with OneLOop [57,58]. We do not show here the one-loop virtual or real radiation diagrams which are lengthy. We use the dipole subtraction method with massive quarks [59] for handling the QCD real corrections.

Decay via electroweak interactions
Within the standard model, the Higgs boson can also decay into four bottom quarks through a cascade decay with an onshell Z boson, H → Z(bb)bb. The corresponding Figure 1: Feynman diagrams at leading order for the Higgs boson decaying into four bottom quarks via Yukawa interactions. Diagrams can be obtained with interchanges of identical particles are not shown for simplicity.
branching ratio turns to be comparable to the one as induced by the Yukawa couplings of bottom quarks and top quarks. At leading order the relevant Feynman diagrams in Feynman-'t Hooft gauge are shown in Fig. 2. Those diagrams mediated by Z boson and goldstone bosons χ must be considered together to form a gauge invariant set. Furthermore, we also include the diagrams mediated by Higgs bosons though their contributions are small, since we would like to keep full bottom quark mass dependence. In principle there are also Feynman diagrams mediated by photons. We do not consider them here since there at next-to-leading order in QCD one will also need to include one-loop QED corrections of decay via Yukawa couplings for consistency. As mentioned earlier dominant contributions to decay via electroweak couplings arise from the resonance region, namely one of the bb pairs lies at Z boson mass pole. Thus one must include finite width effects of the electroweak gauge bosons which may violate gauge symmetry since that will mix contributions from various orders of the EW couplings. In order to preserve gauge symmetry especially at one-loop level in QCD, we use the complex mass scheme [60]. There the masses of W and Z bosons and the electroweak couplings are complex numbers depending on the width of W and Z bosons, and the Lagrangian are manifestly gauge invariant. The squared amplitudes needed for the next-to-leading order QCD calculation are again generated automatically with GoSam 2.0 [54] and are checked against MadGraph5 aMC@NLO [61]. Good agreements are found between the two programs.

Inclusive decay rate
Similar to the total hadronic width, the inclusive decay width of the Higgs boson to four bottom quarks in the limit of infinite top-quark mass includes contributions from the bottom-quark Yukawa coupling, the gluon effective coupling, and their interferences. It can be expressed as [62] H H with and C 1 (µ), C 2 (µ) as given in Eq. (2.4). The leading-order form factors ∆ ij carry up to quadratic dependence on logarithm of the bottom quark mass. Scale dependent terms in the next-to-leading order corrections δ ij can be obtained through renormalization scale invariance, e.g., Furthermore, the contributions via Higgs boson decaying into electroweak gauge bosons can be written as In this case both ∆ ZZ (x) and a ZZ (x) depend on the bottom-quark mass weakly. We do not consider the interferences between decay induced by Yukawa couplings and via the electroweak gauge bosons.  Full mass dependence of factors ∆ ij and a ij can be complicated. We provide their numerical values in Table. 1  Negative sign of ∆ bg shown in Table. 1 indicates a constructive interference between the Feynman diagrams due to bottom-quark Yukawa coupling and those induced by effective coupling with gluons. We found large NLO QCD corrections for ∆ bb , ∆ gg , and the interference ∆ bg contributions, with mild dependence on the mass of the bottom quark. We further calculate the decay branching ratio of the four bottom quark channels assuming a total width of the Higgs boson Γ tot of 4 MeV [64]. We plot the decay branching ratio as a function of the renormalization scale in Fig. 3   Fractional contributions from different terms to the total decay branching ratio to four bottom quarks are shown in Fig. 4 as a function of the bottom quark mass. Decay via electroweak couplings is sub-dominant but has larger contribution than interference of the bottom-quark Yukawa coupling and the gluon effective coupling. Decay due to pure gluon effective coupling is at the level of a few percents of the total decay branching ratio.

Jet cross sections
We consider final state with at least 4 b-tagged jets to separate from other multi-parton hadronic decay modes, for instance, Higgs boson decaying into a bottom quark pair plus two gluons or two light quarks. We use the k T jet algorithm [65] with a resolution parameter y varied between 10 −3 and 0.5. The separation of any two clusters are calculated as with Q 2 = M 2 H /4. Constituents of jets are combined with E scheme by directly summing the four momentums. Flavor of jets are defined by counting net b-quark numbers within the jet, namely a jet with a b quark and a b anti-quark is considered as light-flavor jet [66].
We show the decay branching ratio to 4 b-jets as a function of the jet resolution parameter in Fig. 5 via either Yukawa couplings or electroweak couplings with a bottomquark mass of 4.8 GeV. In the upper panel we plot the branching ratio at LO and NLO.
In the lower inset we show several ratios including the NLO prediction to the LO one for the nominal scale choice and the scale variations of the LO or NLO predictions. The jet rate approaches the inclusive rate of the decay when y goes to 0 since then the four bottom quarks are always fully resolved. It decreases rapidly as the increasing of y especially for the case of decay via Yukawa couplings where the dominant contributions are from quasicollinear region of bb in the phase space. The effects of QCD corrections are similar as in the inclusive decay rate and exhibit mild dependence on the resolution parameter except when close to the endpoint region of the phase space. We found sizable NLO corrections and reduced scale variations for decay via Yukawa couplings. We further summarize the fractional contributions from different channels to the jet rate in Fig. 6. Contributions from ∆ bb are always dominant while the contributions from ∆ ZZ increase with the increasing of the resolution parameter. E bb,inc , and the pair energy asymmetry ∆E bb . In case that the 4 b-jets arise from a cascade decay of two scalars with identical masses, the mass or energy asymmetry peaks at zero values. We show representative results for a choice of the jet resolution parameter y = 0.02 and m b = 4.8 GeV. That corresponds to an angular separation of ∆θ ∼ 0.3 for two jets with the lower energy equals M H /4. We first plot the energy distributions of each individual jet in Figs. 7-10. Energies are always normalized to the mass of the Higgs boson. We show the results for decay via Yukawa couplings and electroweak couplings in parallel with upper panel gives the differential decay branching ratio and lower inset gives ratio of NLO predictions to LO ones and their scale variations respectively. The distributions from two decay channels show clearly differences. For instance, the (sub-)leading jet peaks at slightly lower(higher) values for decay via electroweak couplings. The energy spectrum of the third and fourth jets are broader for decay via Yukawa couplings. Beside of the normalization, the QCD corrections induce changes on shape of the distributions. The spectrum is generally pushed to lower energy region due to the recoiled gluon in QCD real corrections, except for energy of the third and fourth jets in decay via electroweak couplings, where enhancements are found right after the peak region.
We now move to distributions of invariant mass of jet pairs. In Fig. 11 we plot distribution of the highest invariant mass among all jet pairs in the event, M bb,H . The decay via Yukawa couplings show a broad peak around 0.6M H . On another hand decay via electroweak couplings has a narrow peak close to 0.7M H due to the dominant contributions from onshell Z boson. Concerning shape of the distributions, the QCD corrections shift the spectrum to lower mass region for decay via Yukawa couplings. Both sides off the narrow peak are enhanced for decay via electroweak couplings because of the QCD real radiations. For the distribution of the lowest invariant mass M bb,L in Fig. 12, the decay via Yukawa couplings shows a sharp peak right after the mass threshold. In both decay channels the QCD corrections turn to induce a change of spectrum towards high mass regions. We plot the inclusive invariant mass distribution of jet pairs in Fig. 13. They show a much broader  We present various distributions on energy of the jet pairs in Figs. 15-18. The distribution on the highest energy of all jet pairs in Fig. 15 tends to be very similar to the case of highest invariant mass as discussed in Fig. 11 since they are mostly from the same jet pair. The impact of QCD corrections are also similar, namely push the spectrum to lower energy region for the decay through Yukawa couplings and broaden the Z mass peak in the case of electroweak couplings. For the distribution on the lowest energy of all jet pairs in Fig. 16, at LO it is simply a reflection of the distribution on the highest energy of all jet pairs since the two energies add up to the Higgs boson mass. The QCD radiations change the shape of the distribution in a non-trivial way. Due to the same reason the inclusive energy distribution of all jet pairs in Fig. 17 are symmetric with respect to an energy of M H /2 at LO. Thus the decay via electroweak couplings exhibits a triple-peak structure. For their shape the QCD corrections push the distribution to lower energy region in general due to energies carried away by gluon radiations. Energy of jet pairs can also show large asymmetry as presented in Fig. 18. The QCD corrections are almost constant for decay via Yukawa couplings and are largely enhanced towards the tail of the distribution for decay via electroweak couplings.

Exotic decay
There have been proposals to study various exotic decay channels of the Higgs boson at both the LHC [35] and future electron-positron colliders [36]. Higgs boson decays into four bottom quarks is among them and has been employed to search for possible new scalars with mass m H /2 at the LHC [37]. The ATLAS collaboration recently reported an upper limit of about 50% on the decay branching ratio of the Higgs boson to four bottom quarks via two light scalars [37], depending on the mass and lifetime of the new scalar. In this section we first investigate the QCD corrections to this exotic decay channel as we did

Kinematic distributions and QCD corrections
In the calculation of the exotic decay, we assume the light scalar to be CP-even similar to the SM Higgs boson. We renormalize the Yukawa coupling between the bottom quark and the new scalar with the MS scheme. We also assume the new scalar has a small width thus the narrow width approximation can be applied. We are mostly interested in two kinematic distributions, i.e., the inclusive invariant mass distribution and the inclusive energy distribution of b-jet pairs. They are relevant for the exotic decay via new scalars which will show up either as a sharp peak at the scalar mass or a sharp peak at half of the Higgs mass respectively. We plot the normalized distribution of M bb,inc in Fig. 19 for a scalar mass of 20 and 40 GeV. We compare the shape at LO and NLO. For a lower mass of the new scalar, the QCD corrections barely modify the shape except for the tail regions. We observe a slightly broader peak with the QCD corrections for a larger scalar mass. The impact of QCD corrections are similar for E bb,inc shown in Fig. 20.

Comparison to the SM case
For a future electron-positron collider, for instance, CEPC, the Higgs boson can be produced abundantly and collected with a high efficiency. With a center of mass energy of 250 GeV and an integrated luminosity of 5.6 ab −1 , we expect a total number of the Higgs boson of about 1.1 × 10 6 in ZH production. Thus in the SM it predicts a total of about 4000 events for the Higgs boson decaying into four bottom quarks based on the decay branching ratio in Sect. 2. Such a rate can be measured with a precision of 2% if taking into account only statistical uncertainties. The SM decay of course contributes as an important background for searches of the exotic decay via two new scalars, not mentioning other non-resonant processes, e.g., e + e − → Zbbbb. In Fig. 21 we compare the normalized distribution of M bb,inc for the decay to four bottom quarks in the SM and via the new scalars at NLO in QCD. The SM result includes both decays via Yukawa couplings and electroweak couplings. We consider new scalars with a mass of 20, 40 and 60 GeV respectively. We further apply a Gaussian smearing on the distributions since the signal over background ratio strongly depends on the detector energy resolution. In the left and right plots we assume an energy resolution of 1% and 5% respectively. One can see clearly distortions of the signal peak while the background is less modified. In Fig. 22 we present similar comparison for normalized distribution of E bb,inc . Impact of the energy smearing are similar as shown in Fig. 21. The inclusive energy of b-jet pairs show less discrimination power as comparing to the inclusive invariant mass of b-jet pairs especially after taking into account the realistic jet energy resolution.

Conclusions
We study in details decays of the standard model Higgs boson to two bottom quark and anti-quark pairs. The hadronic decay can be triggered either by Yukawa couplings of bottom and top quarks or the electroweak couplings. Both channels are calculated to nextto-leading order in QCD with the former one utilizing effective theory with top quarks integrated out. We found large NLO QCD corrections for decay via Yukawa couplings that lead to a 90% increase of the inclusive partial decay width and residual QCD scale variations of about 20%. On another hand we found moderate NLO QCD corrections of about 12% for inclusive partial decay width of decay via electroweak couplings. We also calculated various jet cross sections and kinematic distributions. The QCD corrections can result in significant changes not only on normalizations but also on shape of various distributions. At future Higgs factory, all hadronic decay channels of the Higgs boson including the one to four bottom quarks can be explored and measured. Especially they can be employed to directly search for new physics beyond the standard model, e.g., possible new light scalars that can induce exotic decay of the Higgs boson to four bottom quarks. Thus we further compare predictions on various kinematic distributions of the decay in the SM and induced by those new scalars. To complete the study on searches for new physics in the four bottom quark decay channel it will be desirable to carry out a refined study with all SM backgrounds included, e.g., those from non-Higgs boson production processes, and with realistic detector coverage and resolution included. Further matching of the NLO calculations with parton showering and QCD hadronizations will also be required. We leave those topics for a future study.