Electroweak corrections to double Higgs production at the LHC

We present the results for the complete next-to-leading order electroweak corrections to $pp \to HH$ at the Large Hadron Collider, focusing on the dominant gluon-gluon fusion process. While the corrections at the total cross-section level are approximately $-4\%$, those near the energy of $HH$ production threshold exceed $+15\%$, and corrections at the high-energy region are around $-10\%$, leading to a shape distortion for the differential distributions. Our findings substantially diminish the theoretical uncertainties associated with this pivotal process, providing valuable input for understanding the shape of the Higgs boson potential upon comparison with experimental measurements.

-After the groundbreaking discovery of the Higgs boson [1,2], investigating the properties of the Higgs particle has emerged as a paramount avenue for unraveling the intricacies of electroweak (EW) symmetry breaking and the standard model (SM).One of the main focuses nowadays at the Large Hadron Collider (LHC) is the meticulous examination of Higgs selfinteractions, as this endeavor is pivotal in elucidating the elusive shape of the Higgs boson potential, which remains unknown.
The exploration of Higgs boson pair production, a process intricately tied to the Higgs trilinear coupling λ HHH , presents a unique opportunity for advancing our comprehension in this domain.With the continuous accumulation of data at the LHC and the potential enhancements in experimental techniques, there is a heightened expectation that λ HHH will be constrained to a reasonable range [3][4][5].The necessity for the modification of the Higgs boson potential arises if the measured value deviates from the SM value.
Theoretically, gluon-gluon fusion [6][7][8] stands as the predominant production mode for Higgs pairs at the LHC, with rates for alternative production modes being at least an order of magnitude smaller.The leading-order (LO) amplitude for gg → HH is loop induced [9][10][11], owing to the absence of interactions between the Higgs boson and gluons in the SM.This lack of interaction makes fixed-order perturbative calculations highly challenging.The comprehensive next-to-leading order (NLO) QCD calculation was accomplished only a few years ago [12][13][14][15], and subsequent to this, the resummation of soft gluons and parton shower effects was incorporated [16][17][18][19].In the heavy top-quark mass limit approximation, the next-to-next-to-next-to-leading order (N 3 LO) QCD corrections have been obtained [20], coupled with soft-gluon threshold resummation [21].This achievement yields a scale uncertainty of approximately 1.7%, with additional uncertainties stemming from our imprecise knowledge of parton distribution functions (PDFs) and α s , accounting for about 3% of the total uncertainty [22].
NLO EW corrections are notably significant in the high-energy region, primarily attributed to the Sudakov effect [23][24][25], which has been demonstrated by many processes.For example, −10% to −30% corrections are observed in the high invariant mass or high transverse momentum region for vector-boson pair production at the LHC [26].Understanding the impact of NLO EW corrections on Higgs pair production, which represents the most substantial source of uncertainties on the theoretical front, is imperative.Additionally, the Higgs quartic coupling, another crucial quantity defining the Higgs boson potential, only emerges at the NLO EW level.Consequently, computing NLO EW corrections for gg → HH has been a focal point in the 2015, 2017, 2019, and 2021 Les Houches precision wish lists [27][28][29][30].
Calculating the NLO EW corrections to gg → HH is inherently more intricate than computing NLO QCD corrections, primarily due to the challenging two-loop Feynman integrals involving a multitude of mass scales.Several attempts to address the EW corrections have been made in the literature [31][32][33][34][35][36], though they are either partial results or applicable only in specific regions.For instance, NLO EW corrections are computed using heavy top-quark mass expansion in Ref. [34], revealing suboptimal convergent behavior in the physical region and corrections as substantial as 65% when the partonic centerof-mass energy is around 260 GeV.
These findings underscore the critical need for a comprehensive computation of the NLO EW corrections to gg → HH to ensure a reliable prediction.This has been addressed in this Letter, encompassing all Feynman diagrams at the two-loop level and incorporating all mass effects.
Amplitudes. -NLO EW corrections typically encompass both real-emission and virtual contributions.However, in the context of the current problem, the process involving one additional photon emission, gg → HHγ, is prohibited by the Furry theorem.Additionally, heavy EW particle emissions can be experimentally distinguished.Consequently, NLO EW corrections for gg → HH exclusively involve virtual contributions, specifically at the two-loop level.The generation of Feynman diagrams and amplitudes is facilitated by FeynArt [37], with some representative Feynman diagrams illus-trated in Fig. 1.The amplitude for g(p 1 )g(p 2 ) → H(p 3 )H(p 4 ), M = ε 1µ ε 2ν M µν with ε iµ denoting gluon polarization vector, satisfies the current conservation relations p 1µ M µν = p 2ν M µν = 0. We have the following decomposition with each term satisfying the same relations: where ∆ µν 5 is linearly dependent on the Levi-Civita tensor, ∆ µν 0 depends on either p µ 1 or p ν 2 , T µν i contain other tensor structures which only have two degrees of freedom, and F i are gauge invariant form factors.As ∆ µν 5 appears only from the two-loop order onward, its contribution vanishes for NLO EW calculations when multiplied by the LO amplitude.The contribution from ∆ µν 0 vanishes due to current conservation relations.The remained tensors T µν 1 and T µν 2 can be chosen as [11] T where p T = ( tû − m 4 H )/ŝ represents the transverse momentum of one of the Higgs with where D represents the generic spacetime dimension used to regularize potential divergences, we define the projectors such that Because ∆ µν 5 contributes nothing to NLO EW corrections, an even number of γ 5 terms is required in each contribution to the form factors. Additionally, we observe that terms containing γ 5 make no contribution in diagrams involving two fermion loops.Consequently, the even number of γ 5 terms can only arise in a single fermion loop, making the choice of the simplest naive γ 5 scheme, adopting the anticommutativity relation {γ 5 , γ µ } = 0 naively, unambiguous in this context.
Utilizing the CalcLoop package [38], the form factors are expressed as linear combinations of scalar Feynman integrals, which are categorized into 3 (116) integral families based on the type of propagators at the one-loop (two-loop) level.The loop integrals in each family are then further reduced to a more manageable basis, referred to as master integrals, through the utilization of Blade [39].Blade harnesses the power of FiniteFlow [40] to solve integration-by-parts relations [41], and it is equipped with the block-triangular form [42,43] to further enhance computational efficiency (for other reduction packages on the market, see Refs.[44][45][46][47][48][49][50][51][52][53][54][55][56]).As the final physical result is finite and insensitive to a small dimensional regulator ϵ with D = 4 − 2ϵ, we set ϵ = 1/1000 from the beginning to the end of the calculation, aligning with the approach proposed in Refs.[57][58][59].This choice avoids dealing with Laurent expansion with respect to ϵ in intermediate steps, thereby significantly reducing computational resources.Through computations utilizing an alternative point ϵ = −1/1000, we observe that the outcomes deviate by a maximum of 0.4% from those obtained with ϵ = 1/1000.This serves as a confirmation of the effective cancellation of divergences.By combining the finite components of both sets of results, we can further mitigate the error to O(ϵ 2 ).
Master integrals are evaluated by numerically solving differential equations with respect to kinematical variables ŝ and t [60][61][62][63][64][65][66][67][68][69] with boundary conditions provided by the AMFlow method [57,58,[70][71][72][73], as described in Ref. [74].Singularities in the master integrals occur when intermediate particles go on shell, corresponding to √ ŝ = 2m H , m W + m t , 2m t , 2m t + m Z or 2m t + m H in the physical region.Analytical continuation is executed by introducing an infinitesimal positive imaginary part to ŝ across these singularities.Armed with this information, we can compute master integrals at any phase space point starting from the nearest computed point by solving the differential equations with the assistance of the differential solver in AMFlow [57].To validate our approach, we compared the results of all master integrals obtained in this manner against direct evaluations using the AMFlow method at several phase space points.We observed at least 20 digits of agreement, ensuring the robustness and accuracy of our calculations.
For the renormalization of masses and fields, we employ the on-shell scheme, while the renormalization of the electromagnetic coupling α is conducted in the G µ scheme [75].Following the renormalization process, we obtain finite results for F 1 and F 2 , and we have explicitly verified the cancellation of divergences at several phase space points.
Among the partial results reported in the literature, those presented in Refs.[34,35] can be readily compared with our own findings.In the case of Ref. [34], by choosing a nonphysical top-quark mass of 244.22 GeV to ensure the convergence of the heavy top-quark mass expansion, our results exhibit agreement to within approximately 0.1%.However, when considering the physical top-quark mass, our comprehensive calculations reveal corrections of about 34% at a chosen phase space point, as compared to the 65% corrections reported in Ref. [34] for the partial result.Regarding the partial results obtained in Ref. [35], our corresponding findings agree to within approximately 0.5%.This minor discrepancy is primarily attributed to differences in the input parameters used in each study.
Total cross sections.-The cross sections of gluon-gluon fusion channel for pp → HH at LO or NLO can be expressed as where t± = m 2 H − ŝ 2 (1∓ 1 − 4m 2 H /ŝ), f g/p (x, µ) denotes the gluon distribution function inside of the proton, µ represents the factorization scale, and σLO(NLO) are given by σLO = where F (0) i and F (1) i correspond to the lowest order and the next-order terms in the α expansion of the form factors.
We use the following SM parameters, and the top-quark mass is chosen as m t = 172.69GeV [76], keeping all other masses and widths to zero.This approximation may result in a slight deviation of a few percent in the determination of the NLO EW corrections, which is, however, insignificant for our intended purposes.α is calculated from the Fermi constant which results in α = 1/133.12= 7.512 × 10 −3 with input G µ = 1.166378 × 10 −5 GeV −2 .The Cabibbo-Kobayashi-Maskawa mixing matrix is set to diagonal.We employ NNPDF3.1 [77] as our default PDF set, with NNPDF31 nlo as 0118 for both LO and NLO calculations.The running of strong coupling α s with two-loop accuracy is provided by the LHAPDF6 library [78] including five active flavors.Our default renormalization and factorization scales are chosen as µ = M HH /2, where M HH is the invariant mass of the produced Higgs pair.
The phase space integration in Eq. ( 8) is carried out and optimized using Parni [79].A total of 3 × 10 5 events are generated at the LO, and the results are cross-checked with MadGraph5 [80].The LO events are stored in the form of modified Les Houches event files [81], enabling reweighting to the NLO events.
We calculate F 1 and F 2 at NLO for 1.8×10 4 reweighted events.This enables us to compute, on one hand, the NLO cross sections and, on the other hand, the K factor.The K factor is defined as the ratio of (differential) cross sections at NLO to that at LO, representing the main focus of this Letter.To compute the K factor, we utilize these 1.8×10 4 reweighted NLO events and the corresponding LO events.We have verified that K factors computed using this method align with those computed from uncorrelated events.It is important to emphasize that, although there are alternative choices for the aforementioned parameters such as PDF or renormalization scale, they have negligible impact on the K factor.To illustrate this point, Table I presents results for µ = m H and µ = p 2 T + m 2 H in addition to µ = M HH /2.The observed differences with varying scale choices are around 20% both at LO and NLO, attributed to the µ dependence of the strong coupling α s and gluon PDF f g/p .This dependence can be mitigated after incorporating higher-order QCD corrections [22].In contrast, the K factor remains rather stable for different choices of µ, as expected.
Differential cross sections.-Differential cross sections offer a wealth of information about physics, both within the SM and in scenarios beyond it.The impact of NLO EW corrections can vary across specific regions of the phase space compared to the full phase space.
As indicated in Table I, the statistical uncertainty for the K factor is smaller than that of the NLO cross section.This discrepancy arises because the differential K factor exhibits a much flatter behavior compared to the differential cross section, enabling the former to get a controllable error with far fewer events for numerical integration.Given that the computation at LO is significantly more economical, we proceed to compute NLO differential cross sections using the following relation: where ∆K is the K factor calculated in a specific phase space region using the same events at LO and NLO, and ∆σ LO is the LO result computed in the same region but using a significantly larger number of events.With the 1.8 × 10 4 reweighted events, we can compute the K factor quite accurately for most bins, except for those with very large M HH or p T .For each of these bins, we compute an additional 400 reweighted events and use them to determine the corresponding K factor.In Fig. 2, we present the M HH distribution.A significant positive correction of approximately +15% is observed in the first bin.In fact, we find that the EW correction for phase space points near the HH production threshold can exceed +70%.A similar result has also been obtained in Ref. [32] where top-Yukawa corrections have been considered, partially in the heavy top-quark mass limit.This can be understood by examining the EW corrections using heavy top-quark mass expansion.As shown in Ref. [34], the leading term in the expansion at NLO is larger than that at LO by m 4 t , which explains the substantial increase near threshold.However, Same as Fig. 2, but for transverse momentum distribution of one of the two Higgs bosons.
above the threshold, the expansion becomes unreliable, and consequently, the enhancement should no longer exist.Indeed, as M HH increases, the K factor decreases dramatically initially and then slows down as it moves away from the threshold.The pattern is similar for the p T distribution in Fig. 3, where the correction is positive initially and subsequently becomes negative.In regions of either large M HH or large p T , we find the NLO EW correction to be approximately −10%.We explicitly checked phase space points with √ ŝ close to 14 TeV and found the corrections to be as substantial as −30% at the matrix element squared level.However, the gluon luminosity is highly suppressed in this region, and thus, it does not contribute significantly to (differential) cross sections.In Fig. 4, we display the rapidity distribution of one of the two Higgs bosons.A nearly flat K factor is observed, approximately 0.96, similar to the total cross section.
Summary.-Double Higgs production is considered the golden channel to probe the Higgs self-coupling, with NLO EW corrections representing the most substantial source of theoretical uncertainties.In this Letter, we, for the first time, compute the complete NLO EW corrections for the production of double Higgs bosons in the dominant gluon-gluon fusion channel at the LHC.The NLO EW corrections are approximately −4% at the total cross-section level, proving to be insensitive to the choice of parameters.For the differential cross sections, the EW corrections can be significant in certain phase space regions.Specifically, for dimensionful observables, EW corrections reach up to +15% at the beginning of the spectrum and −10% in the tail.For dimensionless observable, the rapidity distribution, the NLO EW corrections are flat and approximately −4%.The complete NLO EW corrections are found to be close to the size of N 3 LO QCD corrections in heavy top-quark mass limit [20], which amounts to 8%.
As emphasized in the Introduction, the size of NLO EW corrections constituted the primary source of theoretical uncertainties in this process before this Letter.Our Letter has substantially reduced the theoretical uncertainties by addressing this critical aspect.The remaining theoretical uncertainties at present primarily stem from high-order QCD effects and QCD parameters, accounting for approximately 3% [22].Upon future comparisons of such precise predictions with experimental measurements, the standard model can be further confirmed if agreement is observed, or potential beyond-thestandard-model effects may be indicated if any tension arises.

FIG. 2 .
FIG. 2. Invariant mass distribution of the Higgs pair with √ s = 14 TeV.The upper plot shows absolute predictions, and the lower panel displays the differential K factor with error bars representing statistical errors.
FIG. 3.Same as Fig.2, but for transverse momentum distribution of one of the two Higgs bosons.

FIG. 4 .
FIG.4.Same as Fig.2, but for rapidity distribution distribution of one of the two Higgs bosons.