Full NLO QCD predictions for Higgs-pair production in the 2-Higgs-doublet model

After the discovery of the Higgs boson in 2012 at the CERN Large Hadron Collider (LHC), the study of its properties still leaves room for an extended Higgs sector with more than one Higgs boson. 2-Higgs doublet models (2HDMs) are well-motivated extensions of the Standard Model (SM) with five physical Higgs bosons: two CP-even states h and H, one CP-odd state A, and two charged states H±\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H^\pm _{}$$\end{document}. In this letter, we present the calculation of the full next-to-leading order (NLO) QCD corrections to hH and AA production at the LHC in the 2HDM at small values of the ratio of the vacuum expectation values, tanβ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tan \beta $$\end{document}, including the exact top-mass dependence everywhere in the calculation. Using techniques applied in the NLO QCD SM Higgs pair production calculation, we present results for the total cross section as well as for the Higgs-pair-mass distribution at the LHC. We also provide the top-quark scale and scheme uncertainties which are found to be sizeable.

1 Introduction [1,2] are well motivated extensions of the SM.They belong to the simplest Higgs sector extensions of the SM that, taking into account all relevant theoretical and experimental constraints, are testable at the LHC.In their type II version they contain the Higgs sector of the Minimal Supersymmetric extension of the SM (MSSM) as a special case.Featuring five physical Higgs bosons after electroweak symmetry breaking (EWSB), they represent an ideal benchmark framework for the investigation of various possible new physics effects to be expected at the LHC in multi-Higgs boson sectors.

2-Higgs Doublet Models
The neutral Higgs boson pairs of the 2HDM are dominantly produced via the loop-induced gluon-fusion process gg → φ 1 φ 2 , where φ 1/2 denote scalar or pseudoscalar Higgs bosons of the 2HDM.Only for mixed scalar+pseudoscalar Higgs production the Drell-Yan-type process q q → Z * → A + h/H takes over the dominant role in large regions of the parameter space [3].The topic of our paper is the calculation of the full NLO QCD corrections to scalar Higgs-pair and pseudoscalar Higgs-pair production via gluon fusion within the 2HDM.
In the past the NLO QCD corrections to the gluon-fusion process gg → HH have been calculated within the SM and the MSSM in the heavy-top limit (HTL) [3].This calculation has been extended to the NNLO QCD corrections in the HTL [4][5][6].Quite recently, this level has been extended to the N 3 LO order in the HTL [7][8][9][10].On the other hand finite top mass effects beyond the HTL have turned out to be sizeable [11][12][13][14][15].The inclusion of the related uncertainties due to the scheme and scale dependence of the virtual top mass has been shown to be mandatory, since they dominate the intrinsic theoretical uncertainties [13][14][15].For BSM scenarios, the NLO QCD corrections to all production modes involving scalar and pseudoscalar Higgs bosons are known in the HTL [3], while partial results for the virtual corrections to pseudoscalar Higgs-pair production are known beyond NLO QCD within the HTL [16].
The paper is organised as follows.In Section 2 we introduce the 2HDM and the benchmark point we have selected to obtain our numerical results, then we give a short descrip-tion of the details of our calculation in Section 3. Our results for hH and AA production are presented in Section 4. The theoretical uncertainties are discussed in Section 5, in particular the top-quark scale and scheme uncertainties in Section 5.2.A short conclusion is given in Section 6.

The 2-Higgs Doublet Model
The 2HDM is obtained by extending the SM by a second Higgs doublet with the same hypercharge.We work within the 2HDM version with a softly broken Z 2 symmetry under which the two Higgs doublets Φ 1,2 behave as Φ 1 → −Φ 1 and Φ 2 → Φ 2 .In terms of the two SU(2) L Higgs doublets with hypercharge Y = +1 the most general scalar potential that is invariant under the SU(2) L ×U(1) Y gauge symmetry and that has a softly broken Z 2 symmetry is given by Working in the CP-conserving 2HDM, the three mass parameters, m 11 , m 22 and m 12 , and the five coupling parameters λ 1 -λ 5 are real.The discrete Z 2 symmetry (softly broken by the term proportional to m 2  12 ) has been introduced to ensure the absence of tree-level flavour-changing neutral currents (FCNC).Extending the Z 2 symmetry to the fermion sector, all families of same-charge fermions will be forced to couple to a single doublet so that tree-level FCNCs will be eliminated [2,17].This implies four different types of doublet couplings to the fermions that are listed in Table 1 together with the transformation properties of the fermions.The corresponding 2HDM types are named type I, type II, lepton-specific and flipped.The resulting couplings of the fermions normalised to the SM couplings can be found in [2].After EWSB, the Higgs doublets Φ i (i = 1, 2) can be expressed in terms of their vacuum expectation values (VEV) v i , the charged complex fields φ + i , and the real neutral CP- Table 1 Classification of the Yukawa types of the Z 2 symmetric 2HDM.2nd-4th columns: allowed coupling combinations of Higgs doublet and fermion types; last five columns: Z 2 assignments for the quark doublet Q, the up-type quark singlet u R , the down-type quark singlet d R , the lepton doublet L, and the lepton singlet l R .
even and CP-odd fields ρ i and η i , respectively, as The mass matrices are obtained from the terms bilinear in the Higgs fields in the potential.Due to charge and CP conservation they decompose into 2 × 2 matrices M S , M P and M C for the neutral CP-even, neutral CP-odd and charged Higgs sector.They are diagonalised by the following orthogonal transformations This leads to the physical Higgs states, a neutral light CPeven, h, a neutral heavy CP-even, H, a neutral CP-odd, A, and two charged Higgs bosons, H ± .By definition, m h < m H .The massless pseudo-Nambu-Goldstone bosons G ± and G 0 are absorbed by the longitudinal components of the massive gauge bosons, the charged W ± and the Z boson, respectively.The rotation matrices are given in terms of the mixing angles ϑ = α and β , respectively, and read The mixing angle β is related to the two VEVs as The mixing angle α is given by where (M S ) i j (i, j = 1, 2) denote the matrix elements of the neutral CP-even scalar mass matrix M S .Introducing we obtain [18] in terms of the abbreviation λ 345 ≡ λ 3 + λ 4 + λ 5 (11) and using the short-hand notation s x ≡ sin x etc.
In the minimum of the potential, the following conditions have to be fulfilled, where the brackets denote the vacuum expectation values.This results in the two equations Exploiting the minimum conditions of the potential, we use the following set of independent input parameters of the model, In this work we choose a benchmark point of the 2HDM type I, in which the couplings of the two Higgs doublets to the up-and down-type fermions are equal.The benchmark point of the 2HDM type I that we use in our numerical analysis is given by the following set of input parameters It fulfils all relevant theoretical and experimental constraints.For a description of the constraints, see Ref. [19].

Partonic leading order cross section
As we work in the 2HDM type I, we are dominated by the top-quark loop contributions so that we neglect the bottomquark loops as well as light-quark loops.Note that while we work in the 2HDM type I, we could apply our approximation to the 2HDM (with natural flavour conservation) of any type as long as we work at low tan β values, as the top-quark Yukawa coupling is the same in all 2HDM types.In particular we could apply our approximation to the 2HDM type II and even to the MSSM as long as the squark contributions can be suppressed, which is the case for squark mass above 400 GeV [3].This is typically the case in current MSSM fits to data [20][21][22][23][24].The leading-order (LO) diagrams for hH and AA production, as depicted in Fig. 1 include triangle diagrams, involving a light and heavy CP-even Higgs h, H propagator coupled to the final-state Higgs bosons with various triple Higgs couplings, and box diagrams with two Yukawa couplings.Note, that we focus here on the production of a mixed CP-even and a pure CP-odd Higgs pair.
The analytical results and the numerical method for LO and NLO QCD hh and HH production can be derived from the SM results [11-14, 25, 26] by simple adjustments of the involved Yukawa and trilinear Higgs self-couplings as well as the sum over Higgs-boson propagators.It should be noted that for larger Higgs masses, as e.g. for HH production, the top-mass effects and the associated mass and scheme uncertainties will be larger than for an SM Higgs mass of 125 GeV.
We follow the conventions of Ref. [3] and decompose the cross section into scalar form factors after the application of two tensor projectors on the matrix elements.The partonic cross section σ (gg → φ 1 φ 2 ), with φ 1 φ 2 = hH or AA, can be written as where ) is the strong coupling constant evaluated at the renormalisation scale µ R , and the Mandelstam variables ŝ and t are given by with the scattering angle θ in the partonic c.m. system and where m φ 1 and m φ 2 are the Higgs boson masses, i.e. either The variable m φ 1 φ 2 denotes the invariant Higgs-pair mass.The factor S is a symmetry factor, S = 1/2 for AA production and S = 1 for hH production.The Källen function λ is given by The integrations limits read as The coefficients C h/H △ contain the triple Higgs couplings λ φ 1 φ 2 h/H and the reduced Yukawa couplings g t h/H , which are given by the 2HDM Yukawa coupling modification w.r.t. to the SM top-Yukawa coupling, as well as the CP-even Higgs boson propagators 1 , The coefficient C □ contains only reduced Yukawa couplings to the final-state Higgs bosons, For the various φ 1,2 they are given by Fig. 1 Generic one-loop diagrams for LO Higgs-boson pair production via gluon fusion, gg → φ 1 φ 2 , in the 2HDM type I.The contribution from triple Higgs couplings is marked in red.Note that φ 1 φ 2 = hH or AA.
In the heavy top-limit (HTL) approximation, the form factors reduce to with a = −1 for hH production and a = 1 for AA production.
The full m t -dependence at LO can be found in Refs.[25,26].

Hadronic cross section
The structure of the NLO QCD corrections is very similar to the SM case presented in Refs.[13,14].They include twoloop virtual corrections to the triangle and box diagrams, one-particle-reducible diagrams involving two triangle diagrams connected by a virtual gluon exchange, and one-loop real corrections involving an extra parton in the final state.The partonic contributions are then convolved with the parton distributions functions (PDFs) f i evaluated at the factorisation scale µ F in order to obtain the hadronic cross section.The parton luminosities dL i j /dτ can be defined as with τ = Q 2 /s, s being the hadronic c.m. energy, so that the NLO hadronic differential cross section with respect to Q 2 can be written as with the LO and the virtual and real correction contributions for i j = gg, ∑ q, q qg, and ∑ q q q, z = Q 2 /τs, and the variable /s.We include five external massless quark flavours.The coefficients C virt of the virtual and C i j of the real corrections in the HTL have been obtained in Ref. [3] and are given by where C ∞, hH/AA △△ denotes the contribution of the one-particle reducible diagrams in the HTL with the transverse momentum The functions P gg (z) and P gq (z) are the related Altarelli-Parisi splitting kernels [27], given by with N F = 5 in our calculation.The cross section σLO (Q 2 ) is calculated in the full theory, i.e. taking into account the finite top-quark mass at the integrand-level.The total cross section can be obtained after a final integration over Q between the threshold m φ 1 + m φ 2 2 and the hadronic c.m. energy s.

Virtual corrections
Three generic types of diagrams contribute to the virtual corrections cf.Fig. 2: (i) two-loop triangle diagrams involving the light and heavy scalar Higgs bosons in the s-channel propagators, (ii) one-particle reducible diagrams emerging from two triangular top loops coupling to a single external Higgs boson that are connected by t-channel gluon exchange and (iii) two-loop box diagrams.The diagrams of class (i) consist of off-shell single scalar Higgs production dressed with the trilinear Higgs vertex.The relative QCD corrections coincide with the NLO QCD corrections to scalar Higgs boson production with mass Q and can thus be adopted from the single-Higgs calculation [28][29][30][31][32].The diagrams of class (ii) define the coefficients c 1 , c 2 in Eq. ( 28).The analytical expressions of the coefficients c 1 , c 2 of the one-particle reducible contributions can be obtained from the corresponding Higgs decay widths of φ → Zγ (φ = h, H, A) [33][34][35] with the corresponding adjustments of the involved couplings.
The full top-mass dependence of c 1 , c 2 is given by 2 with τ φ = 4m 2 t /m 2 φ (φ = h, H, A) and λ t = 4m 2 t /t.The generic loop functions are given by In the case of different pseudoscalar Higgs bosons as in more extended Higgs sectors, the coefficient reads These expressions approach the HTL values given in Eq. ( 28).
The involved part of our calculation is the two-loop box diagrams of type (iii).We have used the same method as in Refs.[13][14][15], i.e. we have performed a Feynman parametrisation, end-point subtractions and the subtraction of special infrared terms to allow for a clean separation of the ultraviolet and infrared singularities.For the stabilisation of the 6-dimensional Feynman integrals we have applied integrations by parts to reduce the powers of the singular denominators and performed the integrations with a small imaginary part of the virtual top mass.In order to arrive at the narrowwidth approximation for the virtual top mass, we have used Richardson extrapolations [36] along the lines of our SM calculation of Refs.[13,14].However, here we needed to extend the calculation for scalar Higgs-boson pairs to the case of different final-state Higgs masses.For the calculation of pseudoscalar Higgs-boson pairs, we have used a naive anti-commuting γ 5 matrix at the pseudoscalar vertices, since only even numbers of γ 5 contribute to the (C P-even) virtual corrections diagram by diagram.For this case, we have used the same projectors as in the double-scalar case, since the contributing tensor structures are the same.Since each individual two-loop box diagram is singular for the t integration, we have applied a technical cut at the integration boundaries and included a suitable substitution to stabilise this integration for each diagram.We have checked explicitly that our results do not depend on this technical cut.
The top mass has been renormalised in both the on-shell scheme and in the MS scheme.The on-shell scheme predictions are our default central predictions while the MS scheme predictions are used to calculate the top-quark scale and scheme uncertainties, see below.The strong coupling constant is renormalised in the MS scheme with 5 active flavours.We have obtained finite results for the virtual corrections by subtracting the HTL results as in the SM case so that we end up effectively calculating the NLO mass ef-fects only.To obtain the final hadronic differential cross section, we have added back the HTL results calculated with HPAIR 3 .The calculation of each two-loop box diagram has been performed independently at least twice with different Feynman parametrisations and we have obtained full agreement within the numerical precision.

Real corrections
The calculation of the finite mass effects in the real corrections, ∆ σ mass i j = ∆ σ i j − ∆ σ HTL i j , follows closely the method described in Refs.[13,14] for the SM case.The HTL contributions are calculated again with the program HPAIR while the partonic mass effects are obtained as where the exact four-momenta p i are mapped onto LO subspace four-momenta pi following Ref.[37].
The HTL matrix elements have been calculated analytically, while the full one-loop matrix elements have been obtained by two different methods.They have been generated using FeynArts [38] and FormCalc [39] on the one hand, and obtained analytically using FeynCalc [40] on the other hand.The scalar one-loop integrals have then been calculated numerically using the library COLLIER 1.2 [41].The phase-space has also been parameterised in two different ways.The two methods agree within the numerical precision.

Numerical results
We present our numerical results at a hadron pp collider for c.m. energies of √ s = 13 and 14 TeV (LHC energies), √ s = 27 TeV (high-energy variant of the LHC, the HE-LHC), and √ s = 100 TeV (FCC energy).We use m t = 172.5 GeV for the on-shell top-quark mass.We have performed the calculation using the NLO PDF set PDF4LHC15 [42] as implemented in the LHAPDF-6 library [43].Our central scale choice is µ R = µ F = µ 0 = Q/2, and α s (M 2 Z ) is set according to the chosen PDF set, with an NLO running in the fiveflavour scheme.As done also in the SM calculation [13,14], we have used the narrow-width approximation for the top quark.We use the 2HDM benchmark scenario given in Eq. (16).
We have calculated a grid of Q-values from Q = 259.907(269.422)GeV, for hH production (for AA production), to Q = 1500 GeV, so that we obtain the invariant Higgs-pairmass distributions depicted in Fig. 3 for hH production (left) 3 The program can be downloaded at http://tiger.web.psi.ch/hpair/.and AA production (right), for the LHC at 13 TeV.The results at 14 TeV are shown in Fig. 4, while the results for the HE-LHC are shown in Fig. 5 and the results for the FCC in Fig. 6.The full NLO QCD results are displayed in red, including the numerical errors as well as a band indicating the renormalisation and factorisation scale uncertainties obtained with a standard seven-point variation around our central scale choice (cf.Subsec.5.1).The blue line shows the (Born-improved) HTL prediction, while the yellow line displays the HTL supplemented by the full mass effects in the real corrections only and the green line (including numerical errors) the HTL supplemented by the full mass effects in the virtual corrections only.
The mass effects in the real corrections increase with increasing c.m. energy both for hH and AA final states.In CPeven hH production, they reach a negative peak at around Q = 400 GeV and are of the order of −10% at 13 TeV (of the order of −20% at 100 TeV) before mildly increasing up to around -6% at Q = 1500 GeV at 13 TeV (−14% at 100 TeV).In CP-odd AA production, the behaviour of the mass effects in the real corrections is slightly different.There is also a negative peak around Q = 400 GeV, of the order of −8% at 13 TeV (−14% at 100 TeV), but then it mildly increases before reaching a plateau around Q = 1000 GeV.The mass effects are then practically constant, about −6% at 13 TeV (−11% at 100 TeV).The mass effects in the virtual corrections are negative at large Q values for both hH and AA final states, as expected by the restoration of partial-wave unitarity in the high-energy limit.Combined with the mass effects in the real corrections, the full mass effects reach about −30% (−40%) at Q ≃ 1500 GeV for hH production, at lower c.m. energies (at 100 TeV), while the mass effects in the virtual corrections are smaller for AA production, reaching about −15% (−20% for Q ≃ 1500 GeV, at lower c.m. energies (at 100 TeV).This is the same behaviour that is observed in the SM case [11][12][13][14], albeit with a smaller correction for AA production.Note that the mild increase in the mass effects in the virtual corrections at large Q values for AA production can be attributed to numerical fluctuations.The most striking difference between CP-even and CP-odd pair production can be seen around the t t threshold and below.There is a distortion of the shape that is distinctly different from the SM case and also between hH and AA productions, hence discriminating between the two production channels.
We have also obtained the total cross sections from the differential distributions, using a numerical integration of Q.For Q between 300 GeV and 1500 GeV we have used the trapezoidal method supplemented by a Richardson extrapolation [36] while we use a Simpson's 3/8 rule [44] for Q between 270 GeV and 300 GeV and a simple trapezoid for Q between the threshold and 270 GeV.For the FCC c.m. en-  ergy of 100 TeV we have also included three new Q bins between 1500 GeV and 2500 GeV and add their contribution using a Simpson's rule.Including the numerical errors on the final decimal number, we have obtained the following results for the full NLO QCD total cross sections for hH and AA production in our 2HDM benchmark scenario, using PDF4LHC15 PDF sets, The corresponding results in the (Born-improved) HTL approximation, obtained using the same numerical integration of the Q grid, are The comparison of Eq. ( 33) with Eq. ( 34) gives a ≃ −12% top-mass effect correction at NLO on the total cross section for hH production at LHC energies (≃ −21% at the 100 TeV FCC), and a ≃ −5% correction for AA production at LHC energies (≃ −11% at the 100 TeV FCC).While the mass effects are of similar size as the SM Higgs-pair production for CP-even Higgs bosons, they are smaller for CP-odd Higgs pair production.

Factorisation and renormalisation scale uncertainties
We have estimated the factorisation and renormalisation scale uncertainties using the standard seven-point method.We have varied both the factorisation scale µ F and the renormalisation scale µ R around our central scale choice µ R = µ F = Q/2, by a factor of two up and down while avoiding the choices leading to the ratio µ R /µ F being either greater than two or smaller than one-half.The maximal and minimal cross sections obtained by this procedure are then compared to the nominal cross section obtained with the central scale choice.
The scale uncertainties are similar to what is obtained for SM Higgs pair production [11][12][13][14].They are slightly larger in AA production than in hH production.We have also found the following scale dependences for the differential cross section at 13  minimal and maximal cross sections against the central OS prediction are used to calculate the scale and scheme uncertainties.This procedure has already been used for SM predictions and this gives rise to significant uncertainties that are comparable or even larger than the usual factorisation and renormalisation scale uncertainties [13][14][15].We compare the five predictions (the OS predictions and the three MS predictions) in Fig. 7 at the 13 TeV LHC, in Fig. 8 at the 14 TeV LHC, in Fig. 9 at the 27 TeV HE-LHC, and in Fig. 10 at the 100 TeV FCC.The red lines display the OS full NLO QCD Higgs-pair invariant mass distributions, the blue lines the MS full NLO QCD predictions with m t (m t ), the yellow lines the MS full NLO QCD predictions with m t (Q/4), and the green lines exhibit the MS full NLO QCD predictions with m t (Q).For Q values above Q = 400 GeV, the MS prediction with µ t = Q always leads to the smallest distribution while the maximum at large Q values is given by the OS prediction.The lower panels in each figures display the ratios of the various predictions to our central OS prediction.As in the SM case, we see large deviations at large Q values, ≃ −50% at Q = 1500 GeV for all c.m. energies.We have obtained the following uncertainties at 13 TeV for selected Q values in hH production using PDF4LHC15 parton densities, dσ (gg → hH) dQ As already seen in the SM case, the top-quark scale and scheme uncertainties turn out to be significant, as large or even larger than the factorisation and renormalisation scale uncertainties.For Q > 400 GeV, the maximum cross section is always the OS prediction.
From the differential distributions, we can obtain the topquark scale and scheme uncertainties on the total cross section.We adopt the envelope for each Q-bin individually to build up two maximal and minimal differential distributions and we integrate these distributions over Q using fits of the various distributions which are then numerically integrated.We have arrived at the following top-quark scale and scheme uncertainties for the CP-even hH total cross section, The scale and scheme uncertainties are sizeable and should be included in an uncertainty analysis of the 2HDM Higgspair production cross sections according to the procedure of Ref. [15].

Conclusions
In this work, we have calculated the full NLO QCD corrections to mixed scalar and pure pseudoscalar Higgs-boson pair production via gluon fusion gg → hH, AA within the 2HDM type I, working in our benchmark scenario that is not excluded at the LHC.We have integrated the two-loop box diagrams numerically by performing end-point and infrared subtractions of the contributing Feynman integrals.A numerical stabilisation across the virtual thresholds has been achieved by integration by parts of the integrand to reduce the power of the problematic denominators of the Feynman integrals.The results of the triangle diagrams, involving s-channel scalar Higgs propagators and the corresponding trilinear Higgs couplings, have been adopted from the single-Higgs case.The one-particle reducible contributions emerging from either two single scalar or pseudoscalar Higgs couplings to gluons can be derived from the known results for h, H, A → Zγ with appropriate replacements of the contributing couplings and masses.After renormalising the top mass and the strong coupling, we have subtracted the (Born-improved) HTL to obtain the pure virtual NLO top-mass effects.The real corrections have been computed by generating the full one-loop matrix elements with automatic tools.These have then been connected to suitable subtraction matrix elements in the HTL for the radiation part, but keeping the full LO top-mass dependence.This could be achieved by suitably projected 4-momenta inside the LO sub-matrix elements.This yields the pure NLO top-mass effects of the real corrections.Adding both subtracted virtual and real corrections, we obtain the full NLO QCD top-mass effects that have then been added to the (Born-improved) HTL results of Ref. [3] by using the code Hpair.Very similar to the corresponding SM calculation of Refs.[11][12][13][14][15], we find NLO top-mass effects of about 15-25% (depending on the collider energy) for the total cross sections if the top mass is defined as the top pole mass.For the invariant Higgs-pair mass distribution, the NLO top-mass effects can reach a level 30-40% for large invariant mass values.The larger the hadronic collider energy, the larger NLO top-mass effects emerge.The renormalisation and factorisation scale dependence induces uncertainties at the level of 10-15% for scalar Higgs pairs and 12-17% for pseudoscalar Higgs pairs at NLO, i.e. similar to the SM case.We have studied the additional theoretical uncertainties originating from the scale and scheme choice of the virtual top mass and obtained additional uncertainties of about 5-15% for scalar and about 10% for pseudoscalar Higgs-pair production that are significant and should be included in future Higgs-pair analyses.These uncertainties are larger for distributions at large invariant Higgs-pair masses.

Fig. 3 Fig. 4 Fig. 5
Fig.3Invariant Higgs-pair-mass distributions for Higgs boson pair production via gluon fusion at the 13 TeV LHC as a function of Q using the PDF4LHC15 PDF set, in the 2HDM type I. Left: CP-even hH production.Right: CP-odd AA production.In both panels, the Born-improved HTL results (in blue), HTL results including the full real corrections (in yellow), HTL results including the full virtual corrections (in green, including the numerical error), and the full NLO QCD results (in red, including the numerical error) are depicted.The inserts below display the ratio to the NLO HTL result for the different calculations.The red band indicates the renormalisation and factorisation scale uncertainties for the results including the full NLO QCD corrections.gg → hH at NLO QCD | √ s = 14 TeV | PDF4LHC15 dσ/dm hH [fb/GeV] µ R = µ F = m hH /2 NLO scale uncertainty

Fig. 7
Fig.7Higgs-pair invariant mass distribution at the 13 TeV LHC with different scale and scheme for the top-quark mass, in the 2HDM type I. Left: CP-even hH production.Right: CP-odd AA production.The lower panels display the ratio to the default OS prediction.
gg → hH at NLO QCD | √ s = 13 TeV | PDF4LHC15 dσ/dm hH [fb/GeV] µ R = µ F = m hH /2 m AA[GeV] Top-quark scale and scheme uncertaintiesThe calculation of the NLO QCD corrections has been performed in two different schemes for the renormalisation of the top-quark mass.Our central predictions use the on-shell (OS) scheme with a mass m t = 172.5 GeV both in the Yukawa couplings and in the loop propagators.The MS scheme can instead be used, with an appropriate choice of the topquark mass counterterm.On top of this scheme choice, there is also a scale choice for the renormalisation of the top-quark mass, m(µ t ).To obtain the top-quark scale and scheme uncertainties, we have compared three MS predictions to our central OS prediction, for µ t = Q/4, Q, and µ t at the MS top mass itself, m t (m t ) = 163.02GeVforour choice of the OS top-quark mass value, obtained with an N3LO evolution and conversion of the pole into the MS mass m t (m t ).The gg → hH at NLO QCD | √ s = 13 TeV | PDF4LHC15 dσ/dm hH [fb/GeV] µ R = µ F = m hH /2 MS scheme with m t (m t ) MS scheme with m t (m HH /4) MS scheme with m t (m HH ) OS scheme, m t = 172.5 GeV MS scheme with m t (m t ) MS scheme with m t (m AA /4) MS scheme with m t (m AA ) OS scheme, m t = 172.5 GeV m AA[GeV] MS scheme with m t (m t ) MS scheme with m t (m HH /4) MS scheme with m t (m HH ) OS scheme, m t = 172.5 GeV MS scheme with m t (m t ) MS scheme with m t (m AA /4) MS scheme with m t (m AA ) OS scheme, m t = 172.5 GeV m AA[GeV]Fig.8Same as in Fig.7but for √ s = 14 TeV.gg → hH at NLO QCD | √ s = 27 TeV | PDF4LHC15 dσ/dm hH [fb/GeV] µ R = µ F = m hH /2 MS scheme with m t (m t ) MS scheme with m t (m HH /4) MS scheme with m t (m HH ) OS scheme, m t = 172.5 GeV Ratio to OS m hH [GeV] gg → AA at NLO QCD | √ s = 27 TeV | PDF4LHC15 dσ/dm AA [fb/GeV] µ R = µ F = m AA /2 m A = 134.711GeV MS scheme with m t (m t ) MS scheme with m t (m AA /4) MS scheme with m t (m AA ) OS scheme, m t = 172.5 GeV Ratio to OS m AA [GeV] gg → hH at NLO QCD | √ s = 100 TeV | PDF4LHC15 dσ/dm hH [fb/GeV] µ R = µ F = m hH /2 Full NLO results for different top-quark masses 2HDM type I m h = 125.09GeV m H = 134.817GeV MS scheme with m t (m t ) MS scheme with m t (m HH /4) MS scheme with m t (m HH ) OS scheme, m t = 172.5 GeV Ratio to OS m hH [GeV] gg → AA at NLO QCD | √ s = 100 TeV | PDF4LHC15 dσ/dm AA [fb/GeV] µ R = µ F = m AA /2 A = 134.711GeV MS scheme with m t (m t ) MS scheme with m t (m AAg /4) MS scheme with m t (m AA ) OS scheme, m t = 172.5 GeV m m AA [GeV]