1 Introduction↩︎

Constraining the trilinear Higgs coupling is one of the primary goals of the High Luminosity phase of the LHC. Current constraints from the combined ATLAS and CMS analysis at the \(95\%\) confidence level, \(-0.71 < \lambda/\lambda_{SM} < 6.1\) [1], still allow for significant deviations from the Standard Model prediction, leaving room for beyond-the-Standard-Model physics.

Higgs boson pair production is the primary channel for measuring the trilinear Higgs coupling, as the value of \(\lambda\) directly affects both the total cross section and the shape of differential distributions. Precise theoretical predictions are therefore essential.

QCD corrections to Higgs boson pair production in gluon fusion have been computed at next-to-leading order (NLO) with full top mass dependence [2]–[6], including parton shower matching [7]–[12], and up to \(\text{N}^3\text{LO} +\text{N}^3\text{LL}\) with the highest order contribution calculated in the heavy-top limit [13]–[18]. These calculations reduce the residual QCD scale uncertainty to approximately \(\mathcal{O}(3\%)\). The top-mass renormalisation scheme uncertainty [10], [19], which remains the dominant uncertainty at the level of \(\mathcal{O}(20\%)\), is currently under investigation and has been addressed in Refs. [20], [21].

With the QCD uncertainties substantially reduced, NLO electroweak (EW) corrections become increasingly relevant. These corrections have been shown to contribute at the level of \(\mathcal{O}(5\%)\) to the total cross section and up to \(\mathcal{O}(30\%)\) for differential distributions [22]. In addition to the full EW corrections, several gauge-invariant subsets have been studied in Refs. [23]–[32].

These proceedings review Ref. [33], which presented the first study of EW effects in Higgs boson pair production initiated by a quark–antiquark pair. Although these contributions are negligible at the level of the total cross section, they were shown to lead to significant distortions of differential distributions close to the production threshold.

2 Amplitude Calculation↩︎

We consider Higgs boson pair production with an incoming quark–antiquark pair, \[\begin{align} q \overline{q} \to H H\, , \end{align}\] with massless quarks of the first two generations, \(q \in \{u, d, s, c\}\). Furthermore, we assume a diagonal CKM matrix. Under these assumptions, the Higgs bosons can only be produced through their couplings to EW gauge bosons. Representative Feynman diagrams contributing at leading order (LO) and NLO, separated into real-emission and virtual contributions, are shown in Tab. 1.

Since LO and real-emission amplitudes are generated automatically with GoSam-3 [34]–[36], the remainder of this section focuses on the analytic calculation of the two-loop virtual corrections.

We generate all relevant Feynman diagrams using QGRAF [37] and translate them into analytic expressions with FORM 4.2¹ [39] employing the Feynman rules of Ref. [40]. The amplitudes are computed in both the unitary and Feynman gauge for the EW gauge bosons, finding agreement. These choices are particularly convenient, as they lead to the same set of master integrals.

To treat the \(\gamma_5\) matrices appearing in quark–gauge-boson couplings, we employ naive dimensional regularisation (NDR) [41]. In the present calculation, all \(\gamma_5\) matrices occur on open fermion lines. Consequently, no traces containing \(\gamma_5\) need to be evaluated in \(D\) dimensions, thereby avoiding the ambiguities that may arise in NDR. This allows us to anticommute the \(\gamma_5\) matrices along the fermion lines and attach them to the external spinors, yielding \[\begin{align} \overline{v}(p_2)\, \mathcal{M}\,u(p_1) = \overline{v}_L(p_2)\, \mathcal{M}_L\, u_L(p_1) + \overline{v}_R(p_2)\, \mathcal{M}_R\, u_R(p_1)\, . \end{align}\] The remaining matrix elements \(\mathcal{M}_{L/R}\) are free of \(\gamma_5\) and can therefore be written in terms of a single form factor, \[\begin{align} \mathcal{M}_{L/R} \sim \mathcal{F}_{L/R} \cancel{p}_3\, . \end{align}\] For the evaluation of the Feynman integrals, we first use Reduze 2 [42] and Kira 3 [43]–[47] to apply integration-by-parts (IBP) relations and reduce them to a minimal set of \(51\) master integrals associated with two top-level sectors shown in Fig. 2. The \(45\) integrals related to the double-box topology, Fig. [fig:PL1], were computed by some of us in Ref. [30] using canonical differential equations [48], \[\begin{align} \mathrm{d}\mathbf{J} = \varepsilon \sum_{i = 1}^{51} \mathbb{A}_i \mathrm{dlog}(\alpha_i) \mathbf{J}\, , \end{align}\] where \(\mathbf{J}\) denotes the canonical integrals, \(\alpha_i = \alpha_i(s, t, m_H^2, m_{W/Z}^2)\) are algebraic functions of the kinematic scales, \(\mathbb{A}_i\) are rational matrices and \(\varepsilon\) is the dimensional regulator.

This form was achieved by first constructing a canonical basis using DlogBasis [49] for sectors up to six propagators, and a loop-by-loop approach [50] for the top sectors containing seven propagators. The corresponding differential equations were derived with the help of LiteRed [51] and FiniteFlow [52] and subsequently brought into dlog form with the help of Efforteless [53]. The boundary constants were fixed by evaluating the integrals in the large mass expansion, which corresponds to the limit \(s, t, u, m_H^2 \ll m_{W/Z}^2\).

The remaining six integrals of the second topology, which is shown in Fig. [fig:PL2], involve the same set of integration kernels and were brought into dlog form within the same framework.

Having expressed all master integrals in terms of iterated integrals over logarithmic kernels, we perform additional basis rotations to a set of graded transcendental functions [54], [55]. These functions are constructed by determining linear relations among the integrals order by order in the dimensional regulator \(\varepsilon\) using FiniteFlow. Applying these relations allows us to obtain a minimal set of independent functions at each order in \(\varepsilon\) .

Expressing the amplitude in terms of independent functions removes all artificial poles. The remaining UV poles are treated by renormalising the strong coupling constant \(\alpha_s\) in the \(\overline{MS}\) scheme, and the IR poles are removed using the Catani subtraction scheme [56].

To evaluate the independent functions numerically in the physical region, we solve their differential equations using DiffExp [57]. As an integration path, we choose two straight segments, as schematically shown in Fig. 3. The first segment connects the large mass limit to a reference point in the physical region, \[\begin{align} \left\{s_0, t_0, u_0, m_{H,0}^2, m_{W/Z,0}^2\right\} = \left\{\frac{3125}{256}, - \frac{1875}{512}, - \frac{1875}{512}, \frac{625}{256}, 1\right\}\, . \end{align}\] This part of the path is evaluated once and cross-checked against AMFlow [58], [59]. The reference point is then used as a starting point for the second segment to explore the physical region.

This split of the integration path ensures that all kinematic thresholds, which lie outside the physical region and may otherwise slow down the numerical integration, are crossed only once during the evaluation of the reference point.

The amplitude can be evaluated in this way with a precision of \(16\) significant digits in approximately \(\mathcal{O}(1')\) per phase space point. As this is not sufficiently fast for a direct implementation into \(\texttt{POWHEG-BOX-V2}\) [60]–[62], we construct an interpolation grid comprising \(\mathcal{O}(20\, 000)\) phase-space points. Therefore, we compute a rectilinear grid in the parameters \[\begin{align} \beta = \sqrt{1 - 4\, \frac{m_H^2}{s}} \in [0,1[ \quad\text{and}\quad\cos\theta = \frac{t - u}{s\, \beta} \in [-1, 1]\,, \end{align}\] with a sampling density weighted according to the amplitude squared.

The LO amplitudes as well as the real-emission contributions are evaluated on the fly using GoSam-3.

3 Results↩︎

We calculate differential distributions for a center-of-mass energy of \(\sqrt{S} = 13.6\,\mathrm{TeV}\) using the EW input parameters listed in Tab. 1. The renormalisation and factorisation scales are chosen dynamically as \(m_{HH}/2\), where \(m_{HH}\) denotes the invariant mass of the Higgs boson pair. The scale uncertainties are estimated through 7-point scale variations.

Fig. 4 (a) shows the invariant mass distribution of the Higgs boson pair. The LO and NLO quark–antiquark contributions are shown as dashed blue and dotted orange lines, respectively, together with the LO gluon-fusion channel, represented by the solid green line. The NLO quark–antiquark channel exhibits a significant enhancement close to the production threshold, reaching a relative correction of approximately \(+10\%\) compared to LO gluon-fusion. This effect is mainly driven by real-emissions, since the emitted parton softens the spectrum, which leads to a higher population of the threshold region. At larger invariant masses, the NLO corrections remain positive and reach approximately \(+3\%\).

Fig. 4 (b) shows the transverse-momentum distribution of one of the Higgs bosons. The corrections exhibit a pattern similar to that observed for the invariant-mass distribution: An enhancement at low transverse momentum, predominantly driven by real-emission contributions, and a moderate increase in the high-\(p_T\) tail.

4 Conclusion↩︎

These proceedings reviewed the calculation of the quark–antiquark channel for Higgs boson pair production up to NLO QCD. The amplitudes were computed analytically, using differential equations to express Feynman integrals in terms of iterated integrals over logarithmic kernels. Boundary conditions were obtained from a large-mass expansion, and additional basis rotations were performed to construct a basis of graded transcendental functions. The resulting amplitudes were evaluated on a grid using DiffExp and implemented into POWHEG-BOX for phenomenological studies.

We observe that the quark–antiquark channel has a sizeable impact on differential distributions, even though its contribution to the total cross section is negligible. In particular, the invariant-mass distribution of the Higgs boson pair shows an enhancement of nearly \(+10\%\) in the threshold region, as well as noticeable Sudakov-type effects at high energies.

Future work includes the public release of the POWHEG implementation and a detailed study of the interplay between additional EW and QCD contributions to Higgs boson pair production.

Acknowledgments↩︎

This work is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant no. - TRR 257, and by the Leverhulme Trust, LIP-2021-014.

References↩︎