Higgs boson pair production at NLO in the Powheg approach and the top quark mass uncertainties

We present a new Monte Carlo code for Higgs boson pair production at next-to-leading order in the Powheg-Box Monte Carlo framework. The code is based on analytic results for the two loop virtual corrections which include the full top quark mass dependence. This feature allows to freely assign the value of all input parameters, including the trilinear Higgs boson self coupling, as well as to vary the renormalization scheme employed for the top quark mass. We study the uncertainties due to the top-mass renormalization scheme allowing the trilinear Higgs boson self coupling to vary around its Standard Model value including parton shower effects. Results are presented for both inclusive and differential observables.


Introduction
With the discovery of the Higgs boson [1,2] the study of its potential and self-interactions has become of great interest to the scientific community as an ultimate probe of the mechanism of electroweak symmetry breaking.
In the Standard Model (SM), the Higgs potential in the unitary gauge reads where the Higgs mass (m H ) and the trilinear (λ 3 ) and quartic (λ 4 ) interactions are linked by the relations λ SM 4 is the vacuum expectation value, G µ the Fermi constant and λ is the coefficient of the (Φ † Φ) 2 interaction, Φ being the Higgs doublet field.
The Higgs mass has now been measured at the level of one per mille precision [3,4].The trilinear Higgs self-coupling is hence predicted within the SM and its determination at the Large Hadron Collider (LHC), where it is accessible via the production of Higgs boson pairs, thus provides a probe of the SM.However, the main production mode, gluon fusion, has a very small SM cross-section [5] and sensitivity to the λ 3 SM value has not yet been reached, so that only constraints on λ 3 can be derived for now.However, during the recent years, the interval of allowed values for λ 3 has shrunk significantly and this trend will continue during Run 3 and the High-Luminosity (HL) phase of the LHC.Presently, from the analyses of the decay channels, HH → b bγγ, HH → b bτ τ and HH → b bb b, the ATLAS Collaboration excluded values outside the interval −0.6 < κ λ < 6.6, where κ λ = λ 3 /λ SM  3 , at 95% confidence level (CL) under the assumption that all the other couplings have SM values [6]. 1 A slightly stronger limit can be obtained by combining the information from double Higgs production with the one coming from other processes that are sensitive to λ 3 via nextto-leading order (NLO) electroweak (EW) corrections [8][9][10][11][12][13][14][15][16][17][18][19].Indeed the combination of the information coming from double-Higgs and single-Higgs production yields −0.4 < κ λ < 6.3 at 95% CL [6].
In view of future improvements in the experimental analyses of the Higgs pair production process it is interesting to reappraise the uncertainties in the theoretical prediction of this process.The process is mediated by heavy quarks loops that appear in different topologies: i) the "signal" topology given by triangular diagrams which represent the gluon-fusion production of an off-shell Higgs that then subsequently decays via the trilinear coupling into two on-shell bosons and are therefore sensitive to λ 3 ; ii) the "background" topology that at the leading order (LO) is given by box-like diagrams that do not depend on λ 3 .The LO calculation of the Higgs boson pair production has been performed more than thirty years ago [20][21][22].The NLO QCD corrections were first computed in the infinite top mass, m t , limit (HTL) [23], reweighted by the exact Born amplitude.Later they were supplemented by the inclusion of (1/m 2 t ) n corrections up to n = 6 [24][25][26].Then the HTL evaluation was retained only in the virtual two-loop corrections, while the real ones were computed including the full top mass dependence [27,28].Finally, the full top mass dependence in the virtual correction was obtained by numerical methods [29][30][31][32].At the same time analytic results for the two-loop virtual corrections, valid in specific regions of the phase space, were presented [33][34][35][36].The analytic evaluation of the corrections via a highenergy (HE) expansion was later used to replace the numerical evaluation of Refs.[29,30] in the high-energy region in order to better cover that energy range [37].Later the HE evaluation was also merged with the analytic evaluation of the corrections via a Higgs transverse momomentum, p T , in order to obtain an analytic result that covers the entire phase space [38].
The NLO fixed order calculation with the full top mass dependence was also matched to parton shower programs [39,40] and combined with next-to-next-to-leading order (NNLO) QCD corrections in the HTL limit [41][42][43][44] while keeping the double real corrections in full top mass dependence [45].Results at N 3 LO are available in the HTL [46,47].
The aim of this paper is twofold.On the one side, we present a new Monte Carlo (MC) code for Higgs pair production in the Powheg-Box approach [48,49].The code retains the full top quark mass dependence at NLO and is flexible in the input parameters.In our code it is possible to vary at the same time both the top mass scheme used and the value of λ 3 .This is achieved employing an analytic result for the virtual two-loop corrections instead of a numerical grid as in the MC code of Refs.[39,50].On the other side, we study the uncertainty related to the renormalization scheme used for the top mass for arbitrary values of the trilinear coupling including parton shower effects.
With respect to previous works in the literature where similar analyses were presented our work contains several improvements.In Ref. [50] the dependence on λ 3 of the Higgs pair total cross section and differential distributions was studied including parton shower effects.We improve this analysis by addressing also the dependence on the top mass renormalization scheme.The uncertainties in double Higgs production related to the top mass scheme were analyzed in Refs.[31,32,51].With respect to these studies, that are based on a NLO fixed order calculation, we improved the analysis by including parton shower effects.Furthermore, these previous studies were concentrating on the investigation of the total cross section or the differential cross section as a function of the Higgs-pair invariant mass, M HH , while we have the possibility to study other differential distributions.
In this paper we do not discuss the technical details of our MC code.They can be found in the instruction manual in the Docs directory of the source code tree.The code will be available on the Powheg-Box repository at https://powhegbox.mib.infn.it.
The paper is organized as follows.In section 2 we describe the basic feature of our Powheg implementation of the gg → HH process.In section 3 we present our results for the inclusive cross section and differential distributions using different top-quark-mass renormalization schemes and for several values of λ 3 .We will also compare our results with those of previous analyses.Finally we present our conclusions.

Powheg implementation of gg → HH
In this section we briefly discuss the main characteristics of our implementation of the gluon-fusion Higgs pair production process in the Powheg-Box framework.We first briefly recall the basic features of the Powheg formalism, and the required elements to implement a process in the Powheg-Box framework.We then discuss in more detail our implementation of the virtual two-loop corrections and of the real radiation contributions.
The way how these two contributions are implemented is the main difference between our MC code and the one presented in Refs.[39,50].

The Powheg approach
The Powheg formula to match NLO-QCD accurate calculations with parton showers can be written in a sufficiently general way as The Born squared matrix element is represented as B(Φ B ), with Φ B being the phase space at the leading order.The squared matrix elements of the real emission, i.e. channels with an additional parton with respect to the Born process, can be divided in two sets, according to whether they feature soft/collinear divergences, in which case we denote them as R div (Φ), or not, R reg (Φ).Φ represents the product of the Born and the real emission phase spaces Φ = Φ B Φ r .In the Powheg-Box framework it is possible to split the contribution from the divergent processes R div (Φ) into two terms, R s (Φ) and R f (Φ).The term R s should contain the singular terms, and it is matched with the parton shower using the Powheg Sudakov form factor ∆ s t .The definition of the latter is with t being the shower ordering variable.The t 0 appearing in eq. ( 2) is a lower-scale cutoff.On the other hand, R f , which should be finite, is simply added as shown in eq. ( 2), without any further treatment.The arbitrariness in choosing R f can be used to study the theoretical uncertainties linked to the matching procedure, since different definitions differ by higher-order terms.The non-divergent channels, R reg , are also added without being multiplied by the Powheg Sudakov form factor. Finally, the Bs (Φ B ) is the NLO normalization factor In this formula, Vfin represents IR-and UV-regularized two-loop virtual contribution, while Rs is the IR-subtracted real contribution as defined above.
In our MC code B(Φ B ) is obtained from Ref. [20] where the LO amplitude is presented in terms of the Passarino-Veltman functions [52].The evaluation of the latters is performed using the COLLIER code [53].

Virtual two-loop contribution
The virtual two-loop diagrams, that enter in Vfin (Φ B ), can be assigned to "signal" and "background" topologies as in the LO case.The contribution of the "signal" diagrams (triangular topology) is known analytically including the full top mass dependence adapting results for the production of a single Higgs with virtuality M HH [54,55].In our code we implement the expressions of Ref. [54].The possibility of varying the trilinear coupling is introduced via an additional parameter that rescales the "signal" contribution.
The "background" topologies are of two types: double-triangle and box diagrams.Concerning the former, the diagram topology is the product of two one-loop triangle diagrams and we implement in our code the analytic results derived in Ref. [26] that retain the full top mass dependence.The box diagrams are the most difficult contribution to evaluate.The box diagrams depend on four energy scales, namely ŝ, t, m t , m H , where ŝ, t, and û are the Mandelstam variables which satisfy the condition Alternatively, the scale t can be replaced by the transverse momentum of the Higgs particle, p T .Exact analytic results for two-loop box diagrams with several energy scales is at the verge of what can be obtained with the present computational technology.However, using the method of the expansion of the diagrams in terms of ratios of small energy scales vs. large energy scales, it is possible to obtain an analytic evaluation of these diagrams that is valid in specific regions of the phase space where an hierarchy among the various energy scales present in the diagrams is realized.The method of the expansion in terms of the transverse momentum of the Higgs boson [33] is valid in phase-space regions where | t|/(4m 2 t ) ≲ 1 while the high-energy (HE) expansion method [35] covers the complementary regions of the phase space where | t|/(4m 2 t ) ≳ 1.In the p T expansion it is assumed that the scales associated to m H and to p T are small compared to the scales set by ŝ and m t .Under this assumption, the box integrals are expanded in ratios of small over large scales, and the resulting simplified integrals are written as linear combinations of 52 master integrals (MI) using Integration-by-Parts (IBP) identities obtained with LiteRed [56,57].Among the 52 MI, fifty can be expressed in terms of generalised harmonic polylogarithms while two are elliptic integrals [58].On the other hand, in the HE expansion the two-loop box integrals are first expanded in terms of small m H , then an IBP reduction is performed on the expanded integrals; the resulting MIs are further expanded in the limit m 2 t ≪ ŝ, | t| and expressed in terms of harmonic polylogarithms.
As shown in Ref. [38], in order to merge the two analytic approximations the fixed-order results both in the p T expansion and in the HE expansion have to be extended up to or beyond their border of validity, i.e. t ≃ 4m 2 t , by constructing a [1/1] Padé approximant for the p T expanded result and a [6/6] Padé approximant for the HE-result.This procedure reproduces the numerical values [59] in the grid of Ref. [37], which is implemented in the MC code presented in Refs.[39,50], with an accuracy below the 1% level.
In our MC code the two-loop box contribution is implemented via the analytic expressions of the first three terms in the p T -expansion and of the first thirteen terms in the HE-expansion.With the former terms we construct the [1,1] p T -Padé approximant whose expression is evaluated when a point in the phase space satisfies | t|/(4m 2 t ) < 1 or |û|/(4m 2 t ) < 1.With the latter terms we construct a [6,6] HE-Padé approximant whose expression is evaluated when a point in the phase space lies in the complementary region, | t|/(4m 2 t ) ≥ 1 and |û|/(4m 2 t ) > 1.The evaluation of the (generalised) harmonic polylogarithms is done using the code handyG [60], while the elliptic integrals are evaluated using the routines of Ref. [61].
Finally, our analytic expressions allow to easily change the renormalization scheme employed for the top mass as discussed in Ref. [38].Our results are presented in the on-shell (OS) and the modified minimal subtraction (MS) top-mass scheme.In particular, for evaluating the top mass in the MS scheme, we first convert the OS mass to m MS t (µ t = m OS t ) using the two-loop relation [62], and then run it at two-loop order numerically to the indicated scale µ t [63].

Real radiation contribution
As already discussed in subsection 2.1, in the Powheg-Box framework the real emission channels are assigned to the R div or R reg groups, depending on whether they feature soft/collinear-divergent behavior or not.
In double Higgs production in gluon fusion, the gg → hhg and gq → hhq channels belong to R div , while the q q → hhg channel belongs to R reg .
The implementation of these channels has been achieved using the MadLoop matrix element generator [64].The code generated by MadLoop is interfaced directly with the Powheg-Box.The SM model file shipped with MadLoop has been modified in such a way to include an additional parameter that allows for a rescaling of the Higgs trilinear coupling.
To study the uncertainties in the matching of the NLO calculation with parton showers we use the possibility given by the Powheg-Box framework to split the contribution of the processes in R div into a singular and a finite contribution, respectively R s and R f .The events that contribute to the R f term are called "remnant events".
The separation of R div in R s and R f is achieved dynamically by using a damping factor, D h , via with where p HH ⊥ is the transverse momentum of the two-Higgs system, and the default value in Once this separation has been performed, another freedom present in the Powheg-Box framework is the choice of the shower scale for the remnant events (we recall that varying this scale is a higher order effect).By default this scale is set to the p T of the radiated parton.However, it has been found during the study of single Higgs production in gluon fusion [65,66], that such a choice yields, at large p T , harder tails for the NLO system with respect to the fixed order result.To recover the fixed order behavior, it is possible to choose lower scales, in order to limit the phase space available for further emissions by the shower.
We conclude this section by comparing the performance of our MC code with that of the code ggHH presented in Ref. [50].In the latter the virtual two-loop corrections are implemented via several numerical grids [59] in the Mandelstam variables ŝ and t for fixed values of the top and Higgs mass and α s .An interpolation framework is also provided in order to produce the virtual two-loop amplitude at any point in the phase space.This procedure is faster than our approach to compute the corrections at any point in the phase space via our analytic result.However, it lacks flexibility in the input parameters and, in some regions of the phase space, sometimes the interpolated result is not very accurate [67].
On the contrary, for what concern the real emission contributions our MC code is faster than ggHH.This is due to the different way these contributions are computed.We use the MadLoop code to compute the real radiation matrix elements while in ggHH the same contribution is computed using the GoSam code [68,69].
The net result of these two competing factors is that our MC has an average timing for evaluating one phase-space point slightly shorter than the corresponding timing in ggHH.

Results
In this section, we present our numerical results for a center-of-mass energy √ s = 13.6 TeV.
The values of the input parameters are chosen according to the latest recommendation of the LHC Higgs Working Group (LHCHWG): For our studies we adopt the NNPDF31 nlo as 0118 [70] parton distribution functions (PDF) in a five flavour scheme as reference PDF set for the NLO calculation.Correspondingly, the LO results are obtained using the same set extracted at LO.The value of strong coupling constant is set to be the same as the one used in the PDF, i.e. α s (M Z ) = 0.118.

Inclusive Cross Section
In Table 1, we show the total cross section at 13.6 TeV at LO and NLO for several values of λ 3 adopting different top-quark-mass renormalization schemes, i.e.OS and MS with different scale choices µ t .The values of the renormalization and factorization scales are fixed to be µ C = M HH /2 and the scale uncertainty is estimated from the envelope of a 7-point rescaling of , 1), (2, 1), (2, 2).We find that the NLO corrections are large for each choice of the top-mass renormalization scheme.Moreover, the relative size of the scale uncertainties is essentially the same regardless of the top-mass renormalization scheme employed.We note that going from LO to NLO the relative size of the scale uncertainties is reduced by a factor of ∼ 40%.
As illustrated in the left panel of fig. 1, where we show the inclusive cross section at LO and NLO as a function of the Higgs trilinear coupling for different top-mass scheme, the OS scheme leads to the largest value of the total cross section both at LO and NLO up to κ λ ∼ 3.In the same range of κ λ the MS scheme for µ t = M HH gives the smallest cross section while as κ λ increases tends to give the largest one.The maximum difference between the schemes is obtained around the minimum of the cross-section, which corresponds to the maximal destructive interference between the "signal" and "background" diagrams.At LO, the maximum difference between the schemes amounts to about 20% (for κ λ = 2.4), while it decreases to 10% at NLO (again for κ λ = 2.4).We notice that the exact κ λ value that  gives the minimum of the cross section actually depends upon the top-mass scheme.As expected, going from LO to NLO the various minima get closer.
In the right panel of fig. 1, we show the K-factors, K = σ N LO /σ LO , as a function of the Higgs trilinear coupling for different top-mass scheme.We recall that going from a LO result to a NLO one the scheme dependence is expected to decrease.Then, schemes where the LO prediction are smaller show the largest K-factors.
Let us now comment on our findings for the inclusive cross as a function od λ 3 compared to the same study carried out with the ggHH MC code of Ref. [50] and with the results presented in Ref. [51].
The comparison with ggHH was done adopting for the top and Higgs masses and α s the values chosen in ggHH.While we found excellent agreement within the numerical errors of the MCs for the SM, during the course of this study, we have identified a discrepancy with the existing calculation of Ref. [50] when varying κ λ away from its SM value.We have traced this discrepancy to the two-loop virtual contributions.We contacted the authors of ref. [50] who, using our results, found indeed an issue in their two-loop virtual contributions for values of the trilinear coupling different from the SM one.The authors of ref. [50] provided us with a fixed version of their calculation.Using these new results, we found now agreement between the two codes.
For the comparison with the results of Ref. [51] we adopted the same set of PDF employed in that work.For the SM cross section we are in agreement with the results of Ref. [51] at various center-of-mass energies including their estimate of the scale uncertainty.Instead, concerning the dependence of the cross section upon κ λ , we find an agreement with the results of Ref. [51] at the level of per-mille for negative values of κ λ and for κ λ = 1.For larger positive values of κ λ the agreement starts to deteriorate.At the minimum of the cross section, κ λ ∼ 2.4, we find the maximal discrepancy between our and their evaluation of the inclusive cross section with a difference of several per-cent.This discrepancy is difficult to trace.On one side, the agreement with ggHH, after the fixing by the authors, and the agreement with Ref. [51] for κ λ ≤ 1 gives us confidence in our MC code.On the other side, the fact that the maximum difference is obtain at the minimum of the cross section, where   the destructive interference between "signal" and "background" contributions is more pronounced, let us suspect that the way the inclusive cross section is computed in Ref. [51] from finite size bins, could be not sufficiently accurate in regions of the parameter space where there are strong cancellations.This is also suggested by the fact that we find that the discrepancy with Ref. [51] decreases for large positive value of κ λ .

Differential Distributions
We now present the results for a selection of differential observables.As in the previous subsection we focus on the dependence of the observables upon the choice of the topquark-mass renormalization schemes.First, we consider the SM case, κ λ = 1, then we repeat the analyis for several values of κ λ in order to investigate the interplay between the renormalization scheme choice and the value of the Higgs trilinear coupling.(right) ratio between the MS predictions and the OS one.

Renormalization scheme dependence in the SM
In fig. 2 we show the M HH distribution adopting different renormalization schemes for the top mass.This and the following figures are obtained at the NLO+PS level using our Powheg MC code2 in conjunction with the Pythia 8 shower [71,72].In the top left figure, A, we present the absolute distribution that shows the peak at the opening of the 2 m t threshold.Although not clearly visible in the figure, the position of the peak actually depends on the top mass scheme.In fig. 2 B, the ratio between the MS predictions and the OS one is shown.We notice that for large values of M HH , (M HH ≥ 600 GeV) this ratio is approximately constant with the m t (M HH ) scheme giving the smallest cross section.The MS m t (M HH /4) values are up to M HH ∼ 700 GeV the closest to the OS ones.In this M HH range the m t (M HH /4) scheme is the one where the running of the top mass is the least.Below M HH ∼ 400 GeV the other MS predictions show significant deviations from the OS one especially in the region around the opening of the 2 m t threshold.As already said, the opening of the 2 m t threshold occurs at different values of M HH in the various top mass schemes and therefore the ratio shown is influenced by the position of the peak.The bottom part of fig. 2 shows the K-factor in the various schemes.In fig. 2 C the LO result is evaluated using the LO PDF while in the NLO result the NLO PDF is used.In fig. 2 D instead, the NLO PDF is used both at LO and NLO.As already said going from a LO result to a NLO one the scheme dependence is expected to decrease.Then if an MS scheme has a LO prediction smaller than the corresponding OS one its K-factor is expected to be larger than the OS K-factor and vice versa.Taking into account the information shown in fig. 2 B, one sees that the behavior of the K-factors in the bottom left figure shows exactly this feature.In fig. 2 D, the effects due to the variation of the PDFs are eliminated so that the pure scheme-dependent effects are more manifest.The effect induced by the variation of the PDF is quite mild as can be seen comparing figs.consider two observable which are sensitive to the recoil against jet activity, namely the transverse momentum of the two Higgs system, p ⊥ HH , and the transverse momentum of the leading Higgs, p ⊥ H 1 , that is identified as the final state boson with the largest transverse momentum.
In fig. 3 we show the transverse momentum distribution of the two Higgs system for different top mass renormalization schemes.This observable is sensitive to soft gluon radiation.In the left figure the absolute distributions at NLO+PS are presented which show the suppression in the region p ⊥ HH → 0 where the fixed order NLO results become unreliable.The right plot shows the ratio between the MS predictions for the p ⊥ HH distribution and the OS one.We notice that the MS results are always smaller than the OS one.However, in the small p ⊥ HH region they are all quite close to the OS result while, as the p ⊥ HH increases the MS results tend to be more spread among themselves and with respect to the OS one.The small scheme dependence as p ⊥ HH → 0 is likely related to the shower, that in this region has a relevant effect.scheme dependence as the corresponding plots in fig. 2. In fig. 4 C all the K-factors are quite close up to p ⊥ H 1 ∼ 100 GeV showing that the p ⊥ H 1 predictions in this region are little scheme-dependent.For higher values of p ⊥ H 1 the scheme dependence in more pronounced and the K-factors can be very large.However, fig. 4 D, which is the same as plot C but with LO and NLO distributions computed with NLO PDFs, shows that the choice of the PDF plays an important role and, once the same PDF is used both at LO and NLO, the K-factors become smaller and the scheme dependence is approximately the same for all values of p ⊥ H 1 .

Interplay between a modified trilinear coupling and the renormalization scheme
In this subsection the same analyses presented in the previous subsection are repeated for different values of the Higgs trilinear coupling.For shortness we concentrate on the ratio between the MS predictions and the OS one at NLO+PS.In fig. 5 the invariant mass distribution of the two Higgs system is considered.The four plots in the figure have to be compared with fig. 2 A that correponds to the case κ λ = 1.The figure shows that the cases κ λ = −0.6 (A) and κ λ = 0 (B) are very similar  to κ λ = 1.The case κ λ = 0, where there is no "signal" contribution, indicates that the scheme dependence of the "signal" part is much milder than that of the "background" one.Instead, for κ λ = 2.4 (C) i.e. the λ 3 value where the negative interference between the "signal" and "background" diagrams in the OS scheme is maximal, one sees that the scheme dependence is more pronounced in the region around the 2 m t threshold with repect to the SM case.Finally, the case κ λ = 6.6 (D) shows a much milder scheme dependence, as expected, because in this case the "signal" contribution is very amplified.
In fig.6 we consider the transverse momentum distribution of the two Higgs system.This figure has to be compared with fig. 3 B. The cases κ λ = −0.6 (A) and κ λ = 0 (B) are similar to the SM case although with a slightly smaller spread among the MS predictions that also tend to be closer to the OS one.Instead, the cases κ λ = 2.4 (C) and κ λ = 6.6 (D) show a scheme dependence quite similar for any value of p ⊥ HH and relatively small especially in the κ λ = 6.6 case.
Finally, fig.7 presents the case of the transverse momentum distribution of the leading Higgs to be compared with fig. 4 B. As for the previous observables the cases κ λ = −0.6 (A) and κ λ = 0 (B) are very similar to the SM case.The case κ λ = 2.4 (C) shows a sensible reduction of the scheme dependence in the region p ⊥ HH ≲ 150 GeV, while, as before, for κ λ = 6.6 (D) the scheme dependence is very reduced.

Conclusions
We have presented a new MC code for Higgs boson pair production at NLO in the Powheg-Box approach.The main characteristic of this new code is the flexibility both in the inputs parameters, including the Higgs trilinear coupling, and in the choice of top-mass renormalization scheme.This is obtained employing analytic results for the two-loop virtual contributions instead of numerical grids.
Results are presented for the inclusive cross section and differential observables, including parton shower, for several values of λ 3 and different choices of the top mass.We find, as a general trend, that going from a LO result to a NLO one the top mass scheme dependence is reduced.However, for large invariant mass of the Higgs-pair system or large transverse momentum, the scheme dependence in the SM case, κ λ = 1, can be significant, up to 20%.For other values of κ λ we find similar results but also cases where there is a more pronounced reduction of the top mass scheme dependence.
We have compared our results with those of similar investigations present in the literature [50,51].We find a good agreement for the SM case with all the results present in the literature.For values of the Higgs trilinear coupling different from the SM one our results for the inclusive cross section are in agreement with those of the MC ggHH, after the authors of that code corrected their evaluation of the two-loop virtual contributions.With respect to the results of Ref. [51] we found good agreement only for κ λ ≤ 1.Given the fact that in our MC the cases κ λ ̸ = 1 are obtained just assigning a value different from one to the parameter that multiplies the "signal" contribution we are quite confident in our result.
Finally we notice that our MC code is structured in such a way that can be easily extended to include other beyond the SM effects besides the rescaling of the Higgs trilinear coupling.

Figure 2 :
Figure 2: The invariant mass distribution of the two-Higgs system for different choices of the top-mass renormalization scheme: (A) absolute distributions at NLO+PS; (B) ratio between the MS predictions and the OS one; (C) ratio between the distributions computed at NLO+PS and their LO counterpart (K-factors); (D) same as C but with the LO distributions computed with NLO PDFs.

Figure 3 :
Figure 3: The transverse momentum distribution of the two-Higgs system for different choices of the top-mass renormalization scheme: (left) absolute distributions at NLO-PS; (right) ratio between the MS predictions and the OS one.

Figure 4 :
Figure 4: The transverse momentum of the leading Higgs for different choices of the topmass renormalization scheme: (A) absolute distributions at NLO+PS; (B) ratio between the MS predictions and the OS one; (C) ratio between the distributions computed at NLO+PS and their LO counterpart (K-factors); ( D) same as bottom left, but with the LO distributions computed with NLO PDFs.
2 C with 2 D. The M HH observable is quite insensitive to shower effects.To appreciate the latter we