Fully Differential Higgs Pair Production in Association With a $Z$ Boson at Next-to-Next-to-Leading Order in QCD

We present a fully differential next-to-next-to-leading order calculation of the Higgs pair production in association with a $Z$ boson, which is important for probing the trilinear Higgs self-coupling. The next-to-next-to-leading-order corrections enhance the next-to-leading order total cross sections by a factor of $1.2\sim 1.5$, depending on the collider energy, and change the shape of next-to-leading order kinematic distributions. We discuss how to determine the trilinear Higgs self-coupling using our results.

The trilinear Higgs self-coupling can be measured in two ways.One is the indirect detection from its effect on the single Higgs prductions and decays [27][28][29][30][31][32][33] or the electroweak precision observables [34,35].The other is the direct measurement of the Higgs pair productions at high-energy colliders [6-18, 21-26, 36-39].The dominant production channel at a hadron collider is the gluongluon fusion gg → hh process which involves a top-quark loop.The other channels, including the vector boson fusion, the vector boson associated production and the top quark pair associated production, have relatively smaller cross sections.However, the additional particles in the final state make it easier to distinguish the signal and background events.Actually, the different channels have different characteristics, and thus are complementary to each other and deserve discussion on the same footing.
The vector boson associated production channel, as shown in Fig. 1, is of special interest for several reasons.First, the final-state vector boson provides a characteristic tag of the event so that large quantum chromody- FIG. 1: Sample Feynman diagrams for pp → Zhh production.
namics (QCD) backgrounds can be suppressed.As a result, one can select the events with the Higgs boson pair decaying to b bb b, which has the largest branching ratio, and thus the cross section of this channel is comparable to that of gg → hh with the Higgs boson pair decaying to γγb b [38].Second, all the involved Higgs couplings in this channel are not loop-induced up to next-to-leading order (NLO) (which is the case in the gluon-gluon fusion channel), avoiding the unknown effects of virtual heavy particles.It is more direct and clear to interpret the cross section as a function of the Higgs couplings.Finally, it depends on the value of the Higgs self-coupling in a different form, compared to the gluon-gluon fusion gg → hh channel.In particular, it is sensitive to the Higgs selfcoupling which is larger than the SM value [36,38].Precise theoretical predictions of the vector boson associated Higgs pair productions are crucial for a proper interpretation of experimental data.The total cross sections have been calculated up to next-to-next-to-leading order (NNLO) in QCD [12].However, after experimental cuts are imposed on the final state, it is not clear whether the NNLO QCD corrections are the same over the full phase space.Therefore it is desirable to provide fully differential NNLO QCD calculations.We have presented NNLO differential cross sections of the Higgs pair production in association with a W boson at hadron colliders in a previous work [40].In this Letter, we re-port the first fully differential NNLO QCD calculation of the Higgs pair production in association with a Z boson, which is important for probing the trilinear Higgs selfcoupling at the LHC or future high-energy hadron colliders.In contrast to the pp → W hh process, the cross section of pp → Zhh receives a large contribution from the loop-induced process gg → Zhh starting from NNLO, as shown in Fig. 1(b).This additional contribution has a significant impact on the total and differential cross sections.Our calculation shows that the NNLO corrections indeed change the shape of NLO kinematic distributions.In the peak region of some differential distributions, the NNLO corrections reach up to 80%, compared to NLO results.Therefore, our result is an important ingredient for extracting information on the Higgs self-couplings.

THE METHOD
We perform the (N)NLO QCD differential calculations using the q T subtraction method [41], where q T denotes the transverse momentum of the final-state colorless system.This method is based on the understanding of the cross section near the infrared divergent regions, i.e., q T → 0. When q T is small, the cross section can be factorized as a combination of a hard function, a soft function and transverse-position dependent parton distribution functions (PDFs).The soft function and transverseposition dependent PDFs describe the low-energy dynamics near the infrared divergent region, which are independent of the high-energy dynamics except for the dependence due to collinear anomaly [42], and thus can be considered as universal.Most of them have been obtained up to NNLO, based on which a large number of processes have been calculated at (N)NLO [41,[43][44][45][46][47][48][49][50][51].
In our calculation, we divide the (N)NLO cross section into two parts by a cutoff parameter q cut T .In one part with q T < q cut T , the cross section can be obtained by expanding the transverse momentum resummation formula up to NNLO 1 .The other part with q T > q cut T is just the NLO correction to pp → Zhhj, which can be calculated using standard NLO subtraction method, such as the FKS [52] or the dipole subtraction [53].
Let us first discuss the part with q T < q cut T .We make use of the transverse momentum resummation based on soft-collinear effective theory (SCET) [54][55][56].Since the process of q q → Zhh can be considered as a production of an off-shell Z * boson and its decay to Zhh, the cross section of q q → Zhh production in the small q T region 1 Notice that the gg → Zhh channel does not depend on q cut T and is computed apart from the q T subtraction method.As a result, we neglect this channel in all description and discussion about the q T subtraction method, but include its contribution in the numerical results.can be written as [42] dσ Zhh dq 2 where q T , M and y are the transverse momentum, invariant mass and rapidity of the Zhh system, dΦ 3 the threebody phase space, and f i/N (x, µ) the standard PDF.
The integration lower limits ζ 1 and ζ 2 are determined by S, M, q T , y.The hard function H q q (M, µ) contains the contribution from hard-scale interactions and is extracted by matching the (axial) vector currents in QCD onto effective currents built out of fields in SCET.The two-loop results can be obtained from the hard function of a Drell-Yan process [57].All the q T -dependent terms are contained in the collinear kernel [42] . The function F q q (x 2 T , µ) arises from the effect of collinear anomaly [42].The function I q←i describes the evolution of a parton i to q at fixed x ⊥ , of which the two-loop results can be found in Refs.[58,59] .With all the NNLO ingredients available it is straight forward to perform the integration of q T from 0 to q cut T in Eq. (1).Next, we turn to the other part of the cross section with large q T .Due to the constraint q T > q cut T , there must be at least one parton in the final state.The (N)LO cross section of pp → Zhhj contributes to the (N)NLO cross section of pp → Zhh.Therefore, we need to calculate only the NLO corrections to pp → Zhhj production.This is one of the advantages by using q T subtraction, i.e., the present tools to perform NLO calculations can be utilized without any substantial change.Notice that it is even unnecessary to apply any jet algorithms in the final state and all the infrared singularities are either regularized by the cut q T > q cut T or cancelled between the virtual and real corrections, since the constraint q T > q cut T prevents the momentum of the harder parton in the final state from being arbitrarily soft or collinear to the initial-state momenta.The combination of phase spaces of pp → Zhhj at NLO with large q T and pp → Zhh at NNLO with small q T recovers the whole phase space of pp → Zhh at NNLO.In this work, we use the modified MadGraph5 aMC@NLO [60] to calculate the (N)LO corrections to pp → Zhhj.In the bottom panels, the deviation is defined as σ(q cut T )/σ(q cut T = 10 GeV).The curves are drawn using the linear interpolation method.
Combining the two parts together, we obtain the (N)NLO QCD differential cross section of the process pp → Zhh where q max T is fixed by the partonic center-of-mass energy and the phase space constraints.The cross section of the loop-induced process gg → Zhh is both ultraviolet and infrared finite, and thus there is no need to introduce a new subtraction method to calculate this process.We also use MadGraph5 aMC@NLO [60,62] to compute this contribution.We have compared our results with Ref. [61] and found agreement within uncertainties.

NUMERICAL RESULTS
In the following of this paper, we present numerical results for Zhh production at the proton-proton colliders with √ S = 14 TeV and 100 TeV.The CT14 PDF sets [63] and the associated strong coupling are used throughout our calculation.In particular, we use the LO, NLO, and NNLO PDF sets for the corresponding LO, NLO, and NNLO calculations of the cross section, respectively.The default factorization scale µ F and renormalization scale µ R are both set equal to M to avoid possible large logarithms.The scale uncertainties are estimated by varying the default value by a factor of two up and down.The

SM parameters are given by
where v = 246 GeV is the vacuum expectation value of the Higgs field.When using the q T subtraction method in Eq. ( 3), the first priority is to check that the total cross section is independent of the cut-off parameter q cut T .Figure 2 shows the predictions from the two parts with q T < q cut T and q T > q cut T individually, as well as their sum, at both NLO and NNLO.Here, the loop-induced gg fusion channel is not included in the NNLO result because it is independent of q cut T .The two parts at NLO change dramatically when varying q cut T from 2 GeV to 20 GeV, while the NLO total correction is independent of the cut-off parameter.At NNLO, the two parts change slightly because the O(α 2 s ) q cut T -dependent corrections are almost equal to the O(α s ) ones, but with a relative minus sign.In Fig. 2 the statistical uncertainties of the total cross section are less than 0.2%.Notice that the part with q T < q cut T is only accurate at the leading power of q cut T /M .The power corrections are about (q cut T /M ) 2 ∼ 0.04% for the choice of q cut T = 10 GeV and a typical M ∼ 500 GeV.In the following discussion we choose q cut T = 10 GeV.As a cross check, we have compared the NLO cross section of pp → Zhh calculated by Eq.( 3) with that by the standard NLO program MadGraph5 aMC@NLO and found good agreement.
In Fig. 3, we present the total cross section at different collision energies as well as the K-factors of higherorder corrections.With the increasing of collision energy, the cross section increases significantly.Compared to the leading order (LO) results, the NLO cross sections have much smaller scale uncertainties, especially when the collision energy is large.The NLO K-factors are around 1.26 for 14 TeV ≤ √ S ≤ 100 TeV.The NNLO corrections can enhance the NLO cross section further by a factor of 1.2 ∼ 1.5, depending on the collision energy, but have larger scale uncertainties, about ±5%.This is mainly due to the very large contribution from the loop-induced gg → Zhh process, which is about 13% (14 TeV) ∼ 34% (100 TeV) of the NNLO total cross section as shown in Fig. 3.In order to reduce the theoretical uncertainty, it is desired to include higher-order corrections to this process, which means that one needs to calculate two-loop Feynman diagrams of a 2 → 3 process with three different masses.As far as we know, this is beyond current techniques, and we leave it to future work.By adopting the same input parameters, our calculations of the total cross sections are consistent with those in Ref. [12].
In Fig. 4, we show the transverse momentum p T distributions of the leading Higgs (the one with larger transverse momentum) at the 14 TeV LHC, which have not be obtained in previous calculations.It can be seen that the shape of the p T distribution is nearly unchanged from LO to NLO, but at NNLO the peak region increases more significantly than the tail region; i.e., the shape of the kinematic distribution is changed.Similarly to the total cross section, the NNLO corrections are very large, and the scale uncertainties are also larger than NLO ones because of the contribution from the loop-induced pro-  cess.Figure 5 shows the same distributions at a 100 TeV proton-proton collider.The kinematic features are the same as Fig. 4 except for larger NNLO corrections.In particular, the differential NNLO/NLO K-factor in the peak region is as large as 1.8.
Regarding that a Z boson associated Higgs pair production can be used to probe the trilinear Higgs selfcoupling, we investigate the dependence of the total cross section on the self-coupling at a 100 TeV hadron collider in Fig. 6, where the factor κ is defined as It can be seen that the total cross section changes by about −40% to +100% in the range −5 < κ < 3, compared to the SM prediction, and thus provides a potential method to measure the trilinear Higgs self-coupling.However, there are, in general, two values of κ derived from the same total cross section.For example, the NNLO SM total cross section corresponds to κ = 1 and κ = −3.67.Therefore, we need other observables, e.g., the differential distributions, to ascertain this selfcoupling.Figure 7 shows the normalized NNLO transverse momentum distributions of the Z boson at a 100 TeV collider with κ = 1 and κ = −3.67,respectively.The SM prediction (κ = 1) has a larger peak but a smaller tail compared with the case of κ = −3.67.

CONCLUSIONS
In this paper, we have presented complete NNLO QCD predictions for the total and differential cross sections of Zhh production at hadron colliders using the transverse momentum subtraction method.The NNLO corrections enhance the NLO total cross sections by a factor of 1.2 ∼ 1.5, depending on the collider energy, and change the shape of NLO kinematic distributions.We also investigate the sensitivities of the total and differential cross sections to the trilinear Higgs self-coupling, and find that both of them are needed in order to fix this self-coupling.Our precise predictions will be helpful for the extraction of information on the Higgs trilinear self-coupling in future experimental analyses.

FIG. 2 :
FIG.2:The total cross sections of pp → Zhh production at NLO (left) and NNLO (right) without contribution from loopinduced gg fusion channel.In the bottom panels, the deviation is defined as σ(q cut T )/σ(q cut

FIG. 3 :
FIG.3:The total cross section of pp → Zhh production as a function of the collision energy.The bands denote the scale uncertainties when varying µ = µF = µR by a factor of two.The NNLO total cross sections include the loop-induced gg channel, and the contribution from this channel is also shown in the upper panel individually.

FIG. 4 :
FIG.4:The kinematic distributions of pp → Zhh production at the 14 TeV LHC.h1 denotes the Higgs boson with larger transverse momentum.

FIG. 5 :
FIG. 5: The kinematic distributions of pp → Zhh production at a future 100 TeV hadron collider.h1 denotes the Higgs boson with larger transverse momentum.

FIG. 6 :FIG. 7 :
FIG.6:The total cross section as a function of κ at a hadron collider with √ S = 100 TeV.