Higgs Pair Production at NLO QCD for CP-violating Higgs Sectors

Higgs pair production through gluon fusion is an important process at the LHC to test the dynamics underlying electroweak symmetry breaking. Higgs sectors beyond the Standard Model (SM) can substantially modify this cross section through novel couplings not present in the SM or the on-shell production of new heavy Higgs bosons that subsequently decay into Higgs pairs. CP violation in the Higgs sector is important for the explanation of the observed matter-antimatter asymmetry through electroweak baryogenesis. In this work we compute the next-to-leading order (NLO) QCD corrections in the heavy top quark limit, including the effects of CP violation in the Higgs sector. We choose the effective theory (EFT) approach, which provides a rather model-independent way to explore New Physics (NP) effects by adding dimension-6 operators, both CP-conserving and CP-violating ones, to the SM Lagrangian. Furthermore, we perform the computation within a specific UV-complete model and choose as benchmark model the general 2-Higgs-Doublet Model with CP violation, the C2HDM. Depending on the dimension-6 coefficients, the relative NLO QCD corrections are affected by several per cent through the new CP-violating operators. This is also the case for SM-like Higgs pair production in the C2HDM, while the relative QCD corrections in the production of heavier C2HDM Higgs boson pairs deviate more strongly from the SM case. The absolute cross sections both in the EFT and the C2HDM can be modified by more than an order of magnitude. In particular, in the C2HDM the resonant production of Higgs pairs can by far exceed the SM cross section.


Introduction
With the discovery of the Higgs boson [1,2], the Standard Model (SM) is structurally complete. While the Higgs boson behaves very SM-like, the open questions that cannot be answered within the SM, call for New Physics (NP) extensions. In view of the lack of direct discoveries of particles predicted by extensions beyond the SM (BSM), the precise investigation of the Higgs sector plays an important role [3]. The Higgs self-couplings determine the shape of the Higgs potential. Although the Higgs couplings to the SM particles are very SM-like, the Higgs self-couplings can still deviate substantially from their SM values [4]. The trilinear Higgs self-coupling is directly accessible in Higgs pair production [5][6][7][8]. At the LHC gluon fusion into Higgs pairs provides the largest Higgs pair production cross section [9][10][11]. With a value of 32.91 fb at NLO QCD including the full top quark mass dependence for a Higgs mass of 125 GeV and a c.m. energy of √ s = 14 TeV [12][13][14] this process is experimentally challenging. In BSM models, however, Higgs pair production can be significantly enhanced, see e.g. [7,.
Higgs pair production through gluon fusion is mediated by top and bottom quark triangle and box diagrams already at leading order (LO). The next-to-leading order (NLO) QCD corrections are important and have first been obtained in the limit of large top quark masses [5]. Top quark mass effects are important and have been analysed in [37][38][39][40] with first results towards a fully differential NLO calculation presented in [39]. In Ref. [41] analytic results for the one-particle irreducible contributions to the virtual NLO QCD corrections were presented including finite top quark mass effects. Recently, the NLO QCD corrections have been calculated including the full mass dependence of the top quark in the loops [12][13][14]. The results confirm the relevance of the mass effects, in particular for the differential distributions (for former investigations, see [42,43]). The next-to-next-to-leading order (NNLO) QCD corrections in the heavy quark limit have been calculated in [44][45][46]. Results for differential Higgs pair production at NNLO QCD have been presented in [47]. Soft gluon resummation at next-to-leading logarithmic order within the SCET approach has been performed in [48] and extended to the next-to-next-to-leading logarithmic order in [49], including also the matching to the NNLO cross section. The NLO QCD corrections to Higgs pair production in the MSSM in the heavy top mass limit have first been evaluated in [5]. More recently, analytic results for the contributions from one-and two-loop box diagrams involving top and stop quarks have been obtained in the limit of large loop particle masses in [50]. The NLO QCD corrections to double Higgs production in the singlet-extended SM have been computed in [25] and those in the 2-Higgs Doublet Model (2HDM) in [51], both in the large top mass limit.
A model-independent way to parametrise NP effects realised at a scale well above the scale of electroweak symmetry breaking (EWSB), is given by the Effective Field Theory (EFT) approach where higher-dimensional operators are added to the SM Lagrangian and lead to modifications of the Higgs boson couplings. The impact on Higgs pair production through higher-dimensional operators was analysed in [24,29,[52][53][54][55][56][57][58][59][60][61]. The higher-dimensional operators do not only modify the Higgs couplings to the SM particles but introduce novel couplings not present in the SM, with different effects on the triangle and box diagrams, too. In [29] we computed the NLO QCD corrections to gluon fusion into Higgs pairs including higher dimensional operators in the large top mass limit. While the new operators modify the cross section by up to an order of magnitude, their effect on the relative NLO QCD corrections is only of the order of several per cent. The NLO QCD corrections to composite Higgs pair production in models without and with new heavy fermions have been provided in [30]. Recently, the NNLO QCD corrections in the heavy loop particle limit have been given in [59] for the inclusive as well as the differential cross section including the relevant dimension-6 operators. Since the leading order cross section dominantly factorises, the relative QCD corrections are found to be almost insensitive to the composite character of the Higgs boson and to the details of the heavy fermion spectrum.
In this paper we extend the NLO QCD corrections in the heavy top quark limit to Higgs sectors including CP violation. In the SM CP violation is incorporated in the Cabibbo-Kobayashi-Maskawa matrix and rather small. Models beyond the SM provide additional sources of CP violation, that can be significant while still being compatible with the constraints from electric dipole moments (EDMs), see e.g. [62,63]. CP violation is one of the three Sakharov conditions [64] necessary for baryogenesis. Its discovery in the Higgs sector provides an immediate proof of physics beyond the SM. We investigate CP violation both in a more model-independent EFT approach by adding additional CP-violating dimension-6 operators [65,66] and in a specific benchmark model given by the CP-violating 2-Higgs-Doublet Model (C2HDM) [67][68][69][70][71][72][73][74][75][76].
The organization of our paper is as follows. In section 2 we present the results for the NLO QCD corrections in the EFT approach including CP-violating dimension-6 operators. The subsequent section 3 contains our NLO results in the CP-violating 2HDM. The numerical analysis is presented in section 4. In section 5 we summarise and conclude.

Higgs Pair Production in the EFT including CP violation
Before presenting our analytic results for the NLO QCD corrections to Higgs pair production in the EFT approach including CP violation we introduce our notation.

The EFT including CP violation
By adding higher-dimensional operators to the SM, NP effects that appear at scales far above the EWSB scale, can be parametrised in a model-independent way. In case of a linearly realised SU (2) L × U (1) Y symmetry the Higgs boson is embedded in an SU (2) L doublet H. The leading BSM effects are then parametrised by dimension-6 operators. Note, that even though dimension-8 operators can become more important [56], the investigation of the involved kinematic regions is challenging so that we will neglect them in the following. Adopting the Strongly-Interacting-Light Higgs (SILH) basis the operators that are relevant for Higgs pair production are given by [77], where v is the vacuum expectation value (VEV) v ≈ 246 GeV, M h = 125.09 GeV [78] the Higgs boson mass, M W the W boson mass, y t the top Yukawa coupling constant and g s the strong coupling constant. The gluon field strength tensor G a µν in terms of the gluon fields g a µ and the SU (3) structure constants f abc is given by and its dualG a µν readsG a µν = where µναβ it the totally antisymmetric tensor in four dimensions, normalized to 0123 = 1. The effect of the first three operators in Eq. (2.1) is the modification of the top Yukawa and the trilinear Higgs self-coupling compared to their SM values. The second operator also induces a novel two-Higgs two-fermion coupling [79]. The last two operators parametrise effective gluon couplings to one and two Higgs bosons not mediated by SM quark loops. CP violation is accounted for by the complex part i c u of the Yukawa couplings and the last operator in Eq. (2.1) containing the dual gluon field strength tensor. An estimate of the size of the CP-conserving coefficientsc H ,c u ,c 6 andc g and the most important experimental bounds can be found in [65].
In case of a non-linearly realised EW symmetry with the physical Higgs boson h being a singlet of the custodial symmetry and not necessarily being part of a weak doublet, the non-linear Lagrangian [80] with the contributions relevant for Higgs pair production reads Here, the operators with the coefficientsc t ,c tt ,c g andc gg account for CP violation. A bound onc g has been given in [81]. While in the SILH parametrisation the coupling deviations from the SM are required to be small, the couplings c i in the non-linear Lagrangian can take arbitrary values. The relations between the SILH coefficients and the non-linear ones can be derived from the SILH Lagrangian in the unitary gauge after canonical normalization. They read [56] where α 2 = √ 2G F M 2 W /π, with G F denoting the Fermi constant. We will give results for the non-linear parametrisation in the following and summarise the SILH case in Appendix A.

The NLO QCD Corrections in the EFT
Top and bottom quark loops provide the dominant contributions to gluon fusion into Higgs pairs [10]. In the computation of the NLO QCD corrections in the heavy top quark limit we consistently neglect the bottom quark loops in the following. Their contribution in the SM amounts to less than 1% [5,11]. For the computation of the QCD corrections to Higgs pair production in the large top mass limit an effective Lagrangian can be used that is valid for light Higgs bosons. It contains the Higgs boson interactions derived in the low-energy limit of small Higgs four-momentum. In the case of SM single-Higgs production the K-factor derived in this limit approximates the result obtained with the full mass dependence to better than 5% [82][83][84][85][86]. In Higgs pair production, the low-energy approach works less well and induces an uncertainty of about 15% in the K-factor [12][13][14]. Note, that the top mass effects on the K-factor for models including higher-dimensional operators can also be expected to be of order 10-20%, as the NLO corrections are dominated by soft and collinear gluon effects. The Lagrangian with the required effective Higgs couplings to gluons and quarks can be derived from [5] as (2.6) The factor (1 + 11/4 α s /π) arises from the matching of the effective to the full theory at NLO QCD. Note, that neither the effective couplings to gluons nor the purely CP-odd contributions to the Lagrangian receive this factor. The Feynman rules for the effective couplings between one or two Higgs bosons and two gluons, obtained from this Lagrangian based on the low-energy theorems [87][88][89] are summarised in Fig. 1.  Figure 2 shows the generic diagrams that contribute to Higgs pair production through gluon fusion. Applying the effective Feynman rules of Fig. 1 results in the LO partonic cross section, which can be written as, where µ R denotes the renormalisation scale. The Mandelstam variables read in terms of the scattering angle θ in the partonic center-of-mass (c.m.) system with the invariant Higgs pair mass Q and the relative velocity (2.9) The integration limits are given by cos θ = ±1, i.e. (2.10) The F 1,2 ,F 1,2 , G 1 andG 1 summarise the various form factor contributions with their corresponding coupling coefficients. They can be cast into the form The form factors contain the full mass dependence and have been given in [10]. The triangle form factors for the projection on the CP-even and CP-odd Higgs component, F e ∆ and F o ∆ , are given by F ∆ in appendix A1 and by F A ∆ in A2 of [10], respectively. The box form factors corresponding to the spin-0 gluon-gluon couplings, F e 2 , F o 2 and F m 2 , projecting on a purely CP-even, purely CP-odd and a CP-mixed final state Higgs pair, respectively, are given by F 2 of appendix A1, A3 and A2. Finally, the CP-even, CP-odd and CP-mixed box form factors G e 2 , G o 2 and G m 2 corresponding to the spin-2 gluon-gluon couplings are the G 2 form factors of appendix A1, A3 and A2, respectively. In the heavy quark limit the form factors read The introduced abbreviations are 14) The trilinear self-coupling λ hhh , corresponding to the SM value modified by c 3 , is given by The first terms in F 1 , F 2 and G 1 , respectively, are the SM contributions modified by the rescaling c t of the Yukawa coupling and c 3 of the Higgs self-coupling (contained in C ∆ ). The contributions proportional to c ∆ , c 2 ,c ∆ andc 2 originate from the effective two-gluon couplings to one and two Higgs bosons. The novel 2-Higgs-2-fermion couplings induce the terms coming with c tt and c tt . The form factor contributions proportional toc t , respectivelyc 2 t ,c tt and the ones coming withc ∆ andc 2 are the new contributions due to the admission of CP violation. We recover the following limiting cases for the production of a Higgs pair The NLO QCD corrections to gluon fusion into Higgs pairs, composed of the virtual and the real corrections, are obtained with the help of the effective couplings defined in Fig. 1. Sample diagrams are depicted in Fig. 3. Applying dimensional regularization in d = 4 − 2 dimensions, the involved ultraviolet and infrared divergences appear as poles in . We renormalise the strong coupling constant in the MS scheme with five active flavours, i.e. with the top quark decoupled from the running of α s , in order to cancel the ultraviolet divergences. The sum of the virtual and real corrections cancels the infrared divergences. The remaining collinear initial state singularities are absorbed into the NLO parton densities. These are defined in the MS scheme with five light quark flavours. The finite hadronic NLO cross section can then be cast into the form The individual contributions of Eq. (2.17) read Here s denotes the hadronic c.m. energy and 19) and the Altarelli-Parisi splitting functions are given by [90], 20) with N F = 5 in our case. The factorisation scale of the parton-parton luminosities dL ij /dτ is denoted by µ F . The relative real corrections are not affected by the higher-dimensional operators. The virtual corrections, however, are changed with respect to the SM case due to the overall coupling modifications of the top Yukawa and the trilinear Higgs self-coupling and because of the additional contributions from the novel effective vertices. The coefficient C appearing in the virtual corrections is given by and the transverse momentum squared The last four lines in Eq. (2.21) arise from the third diagram in Fig. 3 (upper), containing the two effective Higgs-two-gluon couplings. The remaining terms originate from the diagrams with gluon loops in Fig. 3 (upper). In line 3, the factor 11/2 arises from the matching of the effective theory to the full theory. This induces the factor (1 + 11α s /(2π)) in the CP-even components of the effective couplings of Fig. 1 for the contributions arising from integrating out the top loops, while the effective couplings not mediated by SM quark loops and the CP-odd components of the couplings are not affected. The factor 6 arises in the virtual corrections to Higgs pair production for the final state projecting on the CP-mixed state, while the purely CP-even and CP-odd final state projections do not exhibit such factor, cf. [5]. Note that we have kept the full top quark mass dependence in the LO amplitude in the derivation of the coefficient C for the virtual corrections.

Higgs Pair Production in the C2HDM
In this section we present the NLO QCD corrections to Higgs pair production in a specific UV complete model. Investigations in well-defined UV complete models complement the EFT approach, as the latter cannot account for NP effects arising from light resonances. Here we resort to the CP-violating 2-Higgs-Doublet Model, the C2HDM. The extension of the SM Higgs sector by a complex Higgs doublet naturally fulfills the constraints from the ρ parameter. In the type II 2HDM, furthermore, the two Higgs doublets couple in the same way to the fermions as in the Minimal Supersymmetric Extensions of the SM (MSSM). The 2HDM Higgs couplings, however, are not constrained by supersymmetric relations and thus entail more substantial deviations from the SM that are still compatible with the data. In this sense, the 2HDM [91,92] is an important benchmark model for the experimental study of the effects of extended Higgs sectors. We briefly summarise the basics relevant for our process and refer to the literature for more details, cf. e.g. [63,76].

The C2HDM
The Higgs potential of a general 2HDM with two SU (2) L doublets Φ 1 and Φ 2 and a softly broken discrete Z 2 symmetry reads The absence of tree-level Flavor Changing Neutral Currents (FCNC) is ensured by the required invariance under the Z 2 transformations Φ 1 → −Φ 1 and Φ 2 → Φ 2 . By hermiticity all parameters in V are real except for the soft Z 2 breaking mass parameter m 2 12 and the quartic coupling λ 5 . For arg(m 2 12 ) = arg(λ 5 ), the complex phases of m 2 12 and λ 5 can be absorbed by a basis transformation leading to the real or CP-conserving 2HDM, in case the VEVs of both Higgs doublets are assumed to be real. Otherwise we are in the CP-violating 2HDM, which depends on 10 real parameters. In the following, we will adopt the conventions of [76] for the C2HDM. The VEVs of the neutral components of the Higgs doublets developed after EWSB can in principle be complex if CP violation is allowed. The relative phase between the VEVs can, however, be rotated away by a global phase transformation in the field Φ 2 [67] so that without loss of generality it can be set to zero. After EWSB the two doublets Φ i (i = 1, 2) are expanded about the real VEVs v 1 and v 2 and we have where ρ i and η i denote the real neutral CP-even and CP-odd fields, respectively, and φ + i the charged complex fields. Requiring the minimum of the potential to be located at induces the minimum conditions where Equations (3.27) and (3.28) can be used to trade the parameters m 2 11 and m 2 22 for v 1 and v 2 , and Eq. (3.29) yields a relation between the two sources of CP violation in the scalar potential, thus fixing one of the ten C2HDM parameters. We introduce the mixing angle β given by which rotates the two Higgs doublets into the Higgs basis [93,94]. Defining the CP-odd field ρ 3 ≡ −η 1 sin β + η 2 cos β (the orthogonal field corresponds to the massless Goldstone boson), the neutral mass eigenstates H i (i = 1, 2, 3) are obtained from the C2HDM basis ρ 1 , ρ 2 and ρ 3 through the rotation The neutral mass matrix is diagonalised by the orthogonal matrix R through We order the Higgs bosons by ascending mass as the mixing matrix R can be parametrised as (3.36) For the 9 independent parameters of the C2HDM we then choose [70] v ≈ 246 GeV , t β , α 1,2,3 , m H i , m H j , m H ± and Re(m 2 12 ) . The m H i and m H j denote any of the masses of two among the three neutral Higgs bosons. The mass of the third Higgs boson is obtained from the other parameters [70]. The analytic relations between the above parameter set and the coupling parameters λ i of the 2HDM Higgs potential can be found in [76]. For α 2 = α 3 = 0 and α 1 = α + π/2 the CP-conserving 2HDM is obtained [68]. The mass matrix Eq. (3.33) then becomes block diagonal, ρ 3 is identified with the pure pseudoscalar Higgs boson A and the CP-even mass eigenstates h and H result from the gauge eigenstates through the rotation parametrised in terms of the angle α, i.e.
For the computation of the Higgs pair production process we need the Higgs couplings to two fermions, the Z couplings to two Higgs bosons and the self-couplings among three Higgs bosons.  We allow one type of fermions to couple only to one Higgs doublet, in order to avoid tree-level FCNC. This is achieved by imposing a global Z 2 symmetry under which Φ 1,2 → ∓Φ 1,2 . There are four phenomenologically different 2HDM types as shown in Table 1. The Yukawa Lagrangian from which the Higgs couplings to fermions are derived, reads with Ψ denoting the fermion fields of mass m f . The CP-even and CP-odd Yukawa coupling coefficients c e (H i f f ) and c o (H i f f ) were derived in [76] and we summarise them in Table 2. In the CP-conserving 2HDM the Z boson would only couple to the CP-mixed Higgs pair combination hA or HA. In the case of CP violation it can couple to any pair of Higgs bosons H i H j . In terms of the SU (2) L gauge coupling g and the cosine of the Weinberg angle θ W the Feynman rule for the coupling Z µ H i H j is given by where the four-momenta of both Higgs bosons, p H i,j , are taken as incoming. The coupling coefficients c(ZH i H j ), parametrised by the mixing matrix elements R ij and the ratio of the two VEVs, tan β, read Note, that the coupling coefficient c(ZH i H j ) becomes zero for i = j. The trilinear Higgs self-couplings λ H i H j H k are quite lengthy and we will not list them here explicitly. In the CPconserving limit for a SM-like Higgs boson h, the trilinear self-coupling approaches 3M 2 h /M 2 Z .

The NLO QCD Corrections in the C2HDM
The diagrams contributing to the LO production of a C2HDM Higgs pair H i H j are depicted in Fig. 4. In contrast to the EFT approach, the cross section does not receive contributions from the effective couplings, obtained from integrating out heavy states. Furthermore, as we have now three CP-violating Higgs states H i , we can have different combinations of Higgs pairs in the final state, and in the first diagram of Fig. 4 we have to sum over all three possible Higgs boson exchanges H k (k = 1, 2, 3). Finally, we have an additional diagram contributing to Higgs pair production where a virtual Z boson couples to the triangle and subsequently decays into a Higgs pair, cf. second diagram in Fig. 4. This diagram does not contribute for equal Higgs bosons in the final state, as the coupling coefficient c(ZH i H j ) vanishes in this case. The LO partonic cross section for the production of the Higgs pair H i H j (i, j = 1, 2, 3) can then be cast into the form Figure 4: Generic diagrams contributing to C2HDM Higgs pair production in gluon fusion at LO.
where a t = 1 denotes the axial charge of the top quark in the loop, and where we have used the short-hand notation In Eq. (3.43) the Γ H k and Γ Z denote the total widths of the Higgs boson H k and the Z boson, respectively. 1 For a given set of input parameters we obtain the total Higgs width with a private version of HDECAY [96,97], adapted to the C2HDM, that will be published in a forthcoming paper. The form factors appearing in Eqs. (3.44) are the same as the ones given in section 2.2, apart from F Z ∆ . They can all be found in the appendix of Ref. [10]. In the computation of the NLO QCD corrections in the heavy quark limit we again consistently neglect the bottom quark loops in the following, and we use the Feynman rules for the effective couplings between one and two Higgs bosons to two gluons in the heavy top quark limit. They are given in Fig. 5. For H i = H j ≡ h, these rules can be obtained from the corresponding Hi iδ ab α s πv Hi Hj iδ ab α s πv 2 ones in the EFT approach, Fig. 1, by making the replacements replaced by the 2HDM result in Eq. (3.42) and the factor C for the virtual corrections given by In Eq. (3.47) we have implicitly assumed summation over same indices. The factors a i andã i are given in Eq. (2.22) and

Numerical Analysis
We have implemented the LO and NLO Higgs pair production cross sections both for the EFT approach including CP violation and for the C2HDM in the Fortran program HPAIR [98].

Impact of CP Violation on NLO QCD Higgs Pair Production in the EFT Approach
The impact of the new CP-violating couplings in the EFT approach on the QCD corrections can be read off Figs. 6-8. 2 They display the K-factor, which is defined as the ratio of the NLO and LO cross sections, K = σ NLO /σ LO , where the parton densities and the strong couplings α s are taken at NLO and LO, respectively. Deviations from the SM K-factor arise both in the virtual and the real corrections. In the virtual corrections they emerge from the terms in the curly brackets of the coefficient C, Eq. (2.21). In the real corrections the different weights in the τ integration due to the modified LO cross section induce deviations from the SM. In Fig. 6   in the SM. The upper plot shows that the CP-violating new interactionc g induces a variation of the K-factor between the SM-value 1.94 and 2.03 in the chosen range 3 . We define the maximal deviation of the K-factor away from the SM value K SM = 1.94, induced by the coupling c x , as The impact on the total cross section is measured by the quantity With these definitions, we find forc g = ±0.15, Both the impact on the K-factor and the total cross section is small. The effect ofc t in the numerator of C, Eq. (2.21), almost cancels against the one in the denominator. The impact of the variation ofc tt finally, is shown in Fig. 8. It is of the per-cent order on the K-factor, with δ K,ctt max = 0.018 (4.57) forc tt = ±1.5. With δ σ,ctt max = 14.23 (4.58) the change of the cross section is substantial. Forc tt = ±0.2, however, it induces with δ σ,ctt = 0.25 similar changes as the other CP-violating couplings.
In summary, the new CP-violating couplings change the K-factor by a few per cent only, while the total cross section itself is affected much more significantly. We have varied here, however, the new couplings only one by one away from the SM values. The combined effect of all dimension-6 couplings, both CP-even and CP-odd, might induce more substantial deviations in the K-factor. same renormalisation scale as here, but another pdf set (MSTW08), we found the SM K-factor Ktot = 1.71.

Impact of CP Violation in the 2HDM on NLO QCD Higgs pair production
After applying the minimisation conditions of the Higgs potential, we are left with 9 input parameters for the C2HDM, which we choose as given in Eq. (3.37). For our numerical analysis we adopt as starting point a scenario that is compatible with all the relevant constraints. It is obtained from the sample generated in Ref. [63], where we used the tool ScannerS [100,101] to perform a scan over the input parameters and check for the experimental and theoretical constraints: The potential has been required to be bounded from below, the EW vacuum has been ensured to be a global minimum and it has been checked that tree-level perturbative unitarity holds. The relevant flavour constraints have been applied and agreement with the EW precision observables has been verified. The mass of one of the neutral Higgs bosons has been set to 125 GeV, and compatibility with the Higgs exclusion bounds for the non-SM Higgs bosons and the individual signal strength fits for the 125 GeV Higgs boson has been checked. Finally, the constraint arising from the measurement of the electron EDM, which is the most constraining of the EDMs, has been taken into account. For the numerical analysis applied here, we resort to the C2HDM type II. While this choice does not affect the couplings involved in Higgs pair production, as the bottom loop contribution has consistently been neglected, this is relevant for the parameter ranges that are still allowed after applying the constraints. We included the latest bound on the charged Higgs boson mass, m H ± > 580 GeV, of Ref. [102], resulting from the recently updated analysis by the Belle collaboration of the inclusive weak radiative B-meson decays [103]. For further details on the scan and the applied constraints, we refer to Ref. [63]. The chosen scenario, finally, is given by we find a pseudoscalar admixture to H 1 of Ψ i = 1.06%. In [63]  Starting from this scenario, we vary α 2 , one of the two angles inducing CP violation, in the range −0.13 ≤ α 2 ≤ 0.15 . not necessarily compatible with all applied constraints. Also the total width of H 3 becomes very large for α 2 ≤ −0.11. As we want to investigate the impact of CP violation on NLO QCD Higgs pair production, we still allow these scenarios for illustrative purposes. With these coupling coefficients we are near the SM case for the CP-even component of the top Yukawa coupling and the variation ofc t is comparable to the one in the EFT approach, where we chose −0.15 ≤c t ≤ 0.15 in Fig. 7. Contrary to the EFT approach, however, in the C2HDM we have additional Higgs bosons. These and the Z boson contribute in the triangle diagrams of Higgs pair production, cf. diagrams 1 and 2 in Fig. 4 Figure 9 shows the masses of the next-to-lightest and heaviest neutral Higgs bosons, H 2 and H 3 , as a function of α 2 . At α 2 = 0 there is a cross-over and the Higgs bosons change their roles: the initially lighter H 2 becomes heavier than H 3 . Still, we stick to our convention and call the heaviest Higgs boson H 3 and the next heavier one H 2 . Thus, we have for −0.13 ≤ α 2 ≤ 0.  Both Higgs bosons are heavy enough to decay on-shell into an H 1 pair, so that we can expect the cross section to be larger than in the SM case. This is confirmed by Fig. 10, which shows the cross section for H 1 H 1 production at NLO QCD as a function of α 2 at a c.m. energy of 14 TeV. The smallest value of the cross section is obtained for α 2 = −0.13. With a value of 604.14 fb it exceeds by far the SM cross section of 38.19 fb. 4 At α 2 = 0.03, we observe a strong increase in the cross section. The largest value, given at α 2 = 0.15, is 28.48 pb. The strong increase at α 2 = 0.03 can be understood by inspecting the H 2 and H 3 branching ratios. They are shown in Fig. 11 for the decays into tt (dashed) and H 1 H 1 (full) for H 2 (red) and H 3 (blue). The cross-over at α 2 = 0, where H 2 becomes heavier than H 3 and they change their roles, is clearly visible by the jump in the branching ratios. As can be inferred from the plot, the H 2 branching ratio into H 1 H 1 strongly increases for α 2 ≥ 0.03. The H 2 mass value here drops below the tt threshold, so that this decay channel gets closed and the branching ratio into H 1 H 1 becomes large and even dominating, as the H 2 couplings to the gauge bosons are suppressed. This increase explains the increase in the Higgs pair production cross section. Also resonant H 3 production with subsequent decay into H 1 H 1 plays a role for positive α 2 although it is much less important. At negative α 2 only the H 2 branching ratio into H 1 H 1 is non-negligible and contributes to the resonant production.
We now turn to the investigation of the K-factor, K = σ NLO /σ LO , which is displayed in Fig. 12 together with the individual K-factors of the virtual and real corrections. Again, in the total K-factor the NLO (LO) cross section is evaluated with NLO (LO) parton densities and α s . The K-factor varies between 1.99 at α 2 = −0.13 and 2.07 at α 2 = 0.15. Between α 2 = 0.03 and 0.04, where the total cross section gets strongly enhanced, the K-factor increases a little bit. The maximum deviation from the SM K-factor is found to be δ K,α 2 max = 0.071 (4.66) for α 2 = 0.15. While the deviation in the K-factor is small, the deviation in the absolute cross section is much more substantial. For α 2 = 0. 15   This exceeds by far the deviations found in the EFT approach, and is due to the resonant production of a heavy Higgs boson, subsequently decaying into  In Table 3 we list for our initial scenario defined in Eq. (4.59) the LO and the NLO cross sections as well as the total K-factor for all final states H i H j (i, j = 1, 2, 3). With increasing mass of the final state Higgs pair the cross sections decrease as expected. For the Higgs mass values of our scenario, the only final state that can include resonant contributions, is H 1 H 2 production: A resonantly produced H 3 subsequently decays into H 1 H 2 . The branching ratio BR(H 3 → H 1 H 2 ) = 0.647 · 10 −4 is, however, very small, so that resonant production does not play a role in this case. 5 All non-resonant cross sections exhibit smaller K-factors than H 1 H 1 production with resonant contributions. They deviate significantly from the SM K-factor 1.94 and lie between 1.43 and 1.86. We furthermore observe that the heavier the final state the smaller becomes the K-factor. These findings are in accordance with previous investigations of the MSSM Higgs sector, where, depending on the final state, the K-factor ranges between about 1.73 and 1.96 [5].

Conclusions
The Higgs sector plays an important role in the search for NP. While in extensions beyond the SM, the 125 GeV Higgs boson needs to have SM-like couplings to the other SM particles, the trilinear Higgs self-coupling and consequently Higgs pair production can still deviate significantly from the SM expectations. CP-violation in the Higgs sector plays an important role to explain the observed baryon-antibaryon asymmetry. In this work we computed the NLO QCD corrections to Higgs pair production including CP violation. We worked in the large top mass limit and performed the calculation on the one hand in the effective field theory approach, where NP effects are parametrised by higher-dimensional operators. On the other hand, we resorted to a specific benchmark model, given by the CP-violating 2HDM.
The various contributions to Higgs pair production are affected differently by the QCD corrections. In the EFT approach including CP-violating effects, we found that the K-factor is changed by several per cent only in the investigated parameter regions that are compatible with the LHC Higgs data. This reflects the dominance of the soft and collinear gluon effects in the QCD corrections. The impact of the novel dimension-6 operators on the absolute value of the cross section is much more important, however, as already found previously in the case of the CP-conserving EFT.
Also in the C2HDM, the K-factor of SM-like H 1 H 1 production varies only by a few percent in the investigated parameter range. For the other possible pair production processes with heavier Higgs bosons in the final state we find smaller K-factors, in accordance with previous findings in the MSSM, representing a specific realisation of the 2HDM model. The total cross sections can, however, be much larger than in the SM. This is due to the possibility of resonant heavy Higgs production with subsequent decay into the Higgs pair final state.
With K-factors between 1.4 and 2.1 in the C2HDM and 1.9 and 2.0 in the EFT approach, the inclusion of the QCD corrections in the gluon fusion process is necessary for reliable predictions of the cross section. . (A.74) As can be inferred from Eqs. (A.70) and (A.73), in the SILH approximation inclusive Higgs pair production is not affected by CP-violating effects at LO in the coupling deviation, i.e. at the dimension-6 level.