Electroweak corrections in the 2HDM for neutral scalar Higgs-boson production through gluon fusion

We have computed the two-loop, electroweak corrections to the production of a light and a heavy neutral, scalar Higgs-boson through the important gluon fusion process in the Two-Higgs-Doublet Model. We provide our results in various renormalization schemes for different scenarios and benchmark points, which will be valuable for experimental studies at the LHC. We describe the technicalities of our two-loop calculation and augment it by a phenomenological discussion. Our results are also applicable to the gluonic neutral, scalar Higgs-boson decays.


Introduction
The Higgs boson, which has been discovered at the Large Hadron Collider (LHC) [1] with a mass of M h = 125.09 ± 0.21 ± 0.11 GeV [2], could be part of an extended Higgs sector. One of the simplest extensions of the Standard Model (SM) is the Two-Higgs-Doublet Model (2HDM), where an extra scalar doublet is added [3]. Such an extension can be an additional source of CP-violation, which may contribute to explain the observed matteranti-matter asymmetry of the Universe. It can also provide a dark matter candidate. Thus, such an extension can assist to address problems which are not solved within the SM. Phenomenological studies of the 2HDM have been performed by the ATLAS [4,5] and CMS [6] collaborations. Several benchmark scenarios for the new parameters which arise due to the addition of a second scalar doublet have been collected by the LHC Higgs Cross Section Working Group in Ref. [7].
Extensions of the SM can strongly modify Higgs-boson production and decay processes, which allows to perform exclusion studies for the new parameters of such model extensions. This requires precise theoretical predictions. For example, the addition of a sequential fourth generation of heavy fermions increases the leading order (LO) cross section of the Higgs-boson production process through gluon fusion already by about a factor of nine [8]. The next-to-leading order (NLO) electroweak corrections in this model have been computed in Refs. [9,10], which helped to exclude this model at the LHC.
In the 2HDM, there are two neutral, scalar Higgs bosons H l and H h . One of them is considered to be the SM-like Higgs boson which has already been discovered at the LHC. In addition, there are two charged Higgs bosons H ± as well as a pseudoscalar Higgs boson Various important Higgs-boson production and decay processes have already been studied at NLO in the 2HDM. For example, NLO electroweak corrections to Higgs-boson production in Higgs strahlung and through vector-boson fusion have been determined in the 2HDM in Refs. [11,12]. The NLO electroweak corrections to the decay of the light Higgs boson of the 2HDM into four fermions have been computed in Ref. [13]. The dominant Higgs-boson production mechanism at the LHC proceeds via gluon fusion. Its precise theoretical knowledge is thus a mandatory task. In the SM, the complete NLO electroweak corrections for the most recent value of the top-quark mass m t = 173.1 GeV [2] amount to 5.1% [14,15]. One can thus expect that the electroweak corrections are also sizable in extensions of the SM. In this work, we will contribute to the effort of studying the 2HDM by computing the two-loop electroweak corrections to the production of neutral, scalar Higgs bosons in gluon fusion. On the one hand, it is important to know the NLO electroweak corrections in the 2HDM to Higgs-boson production in gluon fusion for the already discovered Higgs boson. First results in the special case of the alignment limit (cos (α − β) = 0) have already been presented in Ref. [11]. In addition to the alignment limit, we consider the more complicated, general case of cos (α − β) = 0 in this work. On the other hand, it is important to know the NLO electroweak corrections for the production of the heavy, neutral, scalar Higgs boson H h of the 2HDM, which we present, too. We augment our presentation with the technicalities and obstacles which need to be overcome in order to accomplish this calculation. The results of our computations will be provided in different renormalization schemes for the mixing angles α and β.
The structure of this paper is as follows: In Section 2, we introduce the model in which we perform our calculation. Section 3 describes the different renormalization schemes which we use for the mixing angles. In Section 4, we discuss the computational techniques, which have been used in order to accomplish this aim; and in Section 5 we provide the numerical results. Finally, we close with our summary and conclusion in Section 6. In the appendices, we provide supplementary information on the scale dependence arising from the MS renormalization of the mixing angles α and β as well as on the perturbative behaviour of the coupling constants of the Higgs potential.

The model
We discuss the 2HDM extension of the SM which has two complex, scalar doublet fields Φ i , i = 1, 2. In the generic basis, they are parametrized through where the v i are the vacuum expectation values (vev). We consider the 2HDM with a discrete Z 2 symmetry under the transformation Φ 1 → −Φ 1 , Φ 2 → Φ 2 of the two Higgs doublets. This Z 2 symmetry is important, since it has also implications on the Yukawa sector of the 2HDM , where it suppresses tree-level flavour-changing neutral currents (FCNC), which are experimentally unobserved. Before coming to the Yukawa sector, let us complete the discussion of the Higgs potential. The Z 2 symmetry requirement also reduces the number of terms in the potential. However, we allow for a soft-breaking term [16], which does not induce the FCNC problem. With these constraints, the Higgs potential has the following form The five couplings λ 1 ,...,λ 5 and the soft-breaking parameter m 12 are taken to be real as well as the two masses m 1 and m 2 . For physical applications, we work in the physical basis, where the sector of the potential, which is quadratic in the scalar fields, is diagonalized. This leads to the introduction of the mass eigenstates through a change of basis Here, the symbols H l , H h , H ± ≡ H c and H a are the fields of the physical light, heavy, charged and pseudoscalar Higgs boson, which receive the masses M H l , M H h , M Hc and M Ha , while G 0 and G ± are the neutral and charged would-be Goldstone bosons.
The rotation matrix which performs the diagonalization reads The five couplings in the potential can then be expressed in terms of the Higgs-boson masses and the mixing angles. The explicit formulae are given in Eqs.
Finally, we also introduce the Higgs basis [17] in which only the neutral component of one of the two Higgs doublets, say the first one, acquires a vev. This is achieved by the linear combination Expressing the fields Φ a,b in terms of the fields in the physical basis leads to where we defined the neutral Higgs-boson fields through the following linear combinations It becomes apparent that only the neutral component of Φ a acquires a vev, while in Φ b it does not. We introduce the abbreviations c αβ = cos(α − β) and s αβ = sin(α − β) as well as the vev v = (v 2 1 + v 2 2 ) 1/2 . The choice of phases in the definition of the two doublets Φ i in Eq. (1) leads to s αβ = − 1 − c 2 αβ . An additional property of the Higgs basis is that the Goldstone bosons are all contained in one doublet, Φ a . Note, that the Higgs basis was particularly valuable in Ref. [11] in order to show that an MS renormalization of the mixing angle β can become gauge dependent when using popular tadpole renormalization schemes. In the case of the alignment limit In order to give masses to the fermions in the 2HDM, one has different options to construct Yukawa interactions between the fermionic and scalar fields. In particular, one distinguishes four different types of 2HDM Yukawa terms. In type I, only Φ 2 couples to fermions. In type II, the down-type quarks as well as the leptons couple to Φ 1 , while the up-type quarks couple to Φ 2 . Then, in type Y (or type flipped), Φ 2 couples to up-type quarks and charged leptons, while Φ 1 couples to down-type quarks, and finally, we have type X (or type lepton specific), where Φ 2 couples to quarks, while Φ 1 couples to charged leptons. Within this work we focus on type II. This corresponds to the configuration which is realized in the Minimal Supersymmetric SM. As already mentioned in the beginning of this section, we impose a discrete Z 2 symmetry in order to avoid tree-level FCNC [18]. For the type II 2HDM, down-type quarks and charged leptons need to be odd under this while all other fields remain unchanged. Here, the fields d R and R are the right-handed, down-type quarks and the right-handed charged leptons. All fermions except for the top-quark are taken to be massless such that our results are valid for all types of 2HDM Yukawa terms. This is only justified for small values of t β . Furthermore, we neglect flavour mixing in the following.
As new input parameters of the extended Higgs sector, we have the masses of the heavy, charged and pseudoscalar Higgs bosons cos β sin β (9) and the mixing angles α and β, which we express through Several limits of the new parameters can be defined. Beside the already mentioned alignment limit, the decoupling limit will be studied in this work. In this limit, not only c αβ is equal to zero, but in addition, all new mass scales of the 2HDM are much larger than the electroweak scale [19,7].

Renormalization of the mixing angles α and β
In this section, we consider different renormalization schemes for the mixing angles α and β. We distinguish the MS renormalization, which leads to scale-dependent amplitudes, and the on-shell, p * as well as two process-dependent schemes. The latter guarantee scaleindependent amplitudes. In the process considered in this work, the parameter M sb does not enter at LO. Thus, it does not require renormalization unless the running is taken into account in the MS renormalization scheme.
In Section 5, we will show numerical results for these renormalization schemes. Having different renormalization schemes at hand can be valuable in order to estimate the residual uncertainty due to unknown higher order corrections.

MS renormalization scheme
An MS renormalization of the mixing angles α and β has been defined in Ref. [11]. Here, the MS counterterm δβ MS is defined at the τ −H a −τ vertex, and the MS counterterm δα MS is defined at the τ −H l −τ vertex. The counterterms δα MS and δβ MS are required to render all amplitudes finite. The treatment of tadpoles requires special care when using an MS renormalization of the mixing angles α and β in order to guarantee gauge-independent amplitudes. The tadpoles have been treated in the FJ tadpole scheme, which has been introduced by Fleischer and Jegerlehner for the SM in Ref. [20] and was extended for a general Higgs sector, in particular to the 2HDM, in Ref. [11]. An other application of the FJ tadpole scheme in the 2HDM can be found in Ref. [21]. In the SM, the FJ tadpole scheme corresponds to the β t scheme of Ref. [22]. In this approach, the tadpole contributions are not absorbed into bare physical parameters, which intrinsically assures gauge-independent physical counterterms. A renormalization of the soft-breaking scale M sb in the MS scheme at the H c − H h − H c vertex was introduced in Ref. [11]. For an MS renormalization of M sb see also Ref. [23]. At next-to-leading order, the counterterms δα and δβ of the 2HDM can be obtained from the off-diagonal elements of the fieldrenormalization constants of the Higgs-and Goldstone-boson fields [24]. In particular, for a gauge-independent MS renormalization, this relation reads as a central renormalization scale, see Ref. [13] for a similar scale choice.
In addition to the explicit scale dependence of the amplitude, the renormalization scale dependence, i.e. the running of the parameters, can be taken into account by solving a system of coupled differential equations: The explicit expression for the functions B α , B β and B M sb can directly be obtained from the pole part of the corresponding counterterms and are lengthy in the non-alignment limit.
In the MS scheme the explicit values of the input parameter c αβ and t β need to be defined at a given scale µ. We use two different choices of the scale at which we define these values: the averaged scale µ 0 of Eq. (13) as well as the vacuum expectation value v. For the latter choice we run the values from the scale v to the scale µ 0 with the help of the renormalization group equations (RGE) (14). We perform the scale variation µ 0 /2 and 2µ 0 and run the values of c αβ and t β to these corresponding scales. All results in the MS scheme in Section 5 are presented with the parameters defined at µ 0 . In Appendix C we also compare the different MS schemes for the M * benchmark scenarios, which will be introduced in Section 5.

On-shell renormalization scheme
We have also considered the on-shell (OS) scheme, where the renormalized mixing angles are defined such that they diagonalize the radiatively corrected Higgs mass matrices introduced in Section 2. This connects the counterterms of the mixing angles α and β to the off-diagonal terms of the on-shell two-point functions of the CP-even Higgs bosons and of the CP-odd Higgs boson and the neutral Goldstone boson, respectively. The mixing two-point functions are evaluated on-shell, which guarantees scale independence of the rotation matrices in Eq. (4), see Ref. [25]. However, using finite parts of the mixing two-point functions to define the renormalization constants of α and β leads to gaugedependent renormalization constants in addition to the gauge-dependent tadpole terms discussed in Ref. [11]. The problem has been overcome in Ref. [21] using the pinch technique [26,27,25]. Recently, it has been treated in a more general framework in Ref. [12] within the background field method (BFM), leading to the same result. In the pinch technique the two-point function Σ SS (p 2 ) of two scalars S and S for the R ξ -gauge with ξ = 1 is given by where Σ 1PI SS (p 2 ) is the one-particle irreducible (1PI) mixing two-point function. In the 2HDM the additional contribution Σ add SS (p 2 ) reads [21] Σ add with the cosine of the weak mixing angle c W , the gauge coupling g and the one-loop, scalar, two-point integral B 0 [28]. The counterterms to the mixing angles are then given by where the tadpole counterterms t H h H l and t G 0 Ha are given in the Appendix of Ref. [11]. In an alternative on-shell scheme, the counterterm δβ can also be defined through the charged Higgs-boson-Goldstone-boson mixing two-point functions.
These schemes can be expected to produce small perturbative corrections to physical observables, because the finite counterterms of the mixing angles α and β cancel large terms originating from the Z-factors of the neutral, scalar sector and in the pseudoscalar/charged sector as mentioned in Ref. [12].
The p * scheme is derived from the same mixing two-point functions of Eq. (15), but computed at the external momentum such that the rotation matrices in Eq. (4) also remain scale-independent, see Ref. [25]. The additional terms in Eqs. (16) vanish and the two-point functions in Eq. (15) are equal to the 1PI two-point functions in the R ξ -gauge for ξ = 1. This results in the counterterms of the mixing angles Process-dependent renormalization scheme 1 A process-dependent renormalization of the parameters α and β can be obtained by imposing renormalization conditions on the physical decay processes H h → τ + τ − and H a → τ + τ − [21]. The two counterterms δα and δβ are fixed by the requirement that their NLO corrected partial decay width Γ NLO weak , which contains weak corrections only, is equal to the LO partial decay width Γ LO , i.e.
Since both vertices (τ −H h −τ and τ −H a −τ ) depend on t β , see Fig. 1 for illustration, the determination of both counterterms δα and δβ requires the solution of a linear system of equations. First, δβ can be determined from the τ −H a −τ vertex, which depends only on t β . Then, the result needs to be inserted into the τ − H h − τ vertex for the determination of δα. We refer to this scheme as proc1 in the following.
If the additional 2HDM particles are discovered, one can expect that these two processes can be measured precisely due to the leptonic final states. Therefore, they are in principle well suited for the determination of the parameters α and β. Due to the choice of the renormalization condition, no higher order electroweak corrections are present in the determination of α and β. Process-dependent renormalization may lead to large perturbative corrections, since the higher order corrections to the processes used for renormalization will appear in all other processes. Such a behaviour was, for example, observed for the decay process H ± → W ± H l in Ref. [21]. For light and heavy Higgs-boson production through gluon fusion, however, process-dependent renormalization leads to results similar to those obtained in the OS or p * scheme, as we will see in Section 5.

Process-dependent renormalization scheme 2
Unlike the previously discussed physical renormalization conditions, which rely on the partial decay widths of two leptonic Higgs-boson decays, we consider a process dependent renormalization which is based on the two vertices τ − H a − τ and Z − H h − Z in the following. Studying a potential decay of a heavy Higgs boson into two Z-bosons is experimentally less clean, since the unstable Z-bosons can further decay into a variety of other particles. In addition, its LO contribution is very small in physically relevant scenarios close to the alignment limit. Nevertheless, it turns out, as we will see in Section 5, that in the calculation of the processes gg → H l and gg → H h , this scheme leads to higher order corrections of moderate size.
As renormalization condition, we require that the purely weak corrected τ − H a − τ vertex is equal to its leading order value in order to fix the counterterm δβ. We impose a similar condition on the Z − H h − Z vertex in order to fix the counterterm δφ. Note that due to our choice of parameters, t β and c αβ , it is more convenient to renormalize φ = α − β rather than α. The Z − H h − Z vertex is proportional to g µν at LO, see Fig. 1. Computing higher order corrections to this vertex leads to a richer tensor structure, which also contains combinations of the four-momenta of the external Z-bosons. As renormalization condition for the counterterm δφ, we require that the coefficient in front of g µν of the Z − H h − Z amplitude is equal to its leading order value similar to the renormalization of the electric charge in QED. The coefficient in front of g µν can be extracted from the amplitude by using a suitable projector. This scheme will be called proc2 in the following. The use of the vertices τ − H a − τ and Z − H h − Z has the advantage that one has a separate condition for each of the two counterterms δβ and δφ, i.e. the value of δφ is not influenced by the renormalization condition for δβ, like it is the case for the previously discussed proc1 scheme.

Calculation
The partonic, leading order cross section for the Higgs-boson production process through gluon fusion reads where s is the square of the sum of the external gluon (g) momenta, α s is the strong coupling constant and G F is the Fermi coupling constant. The leading order amplitude where the index q runs over all quark flavours and The symbol H in Eq.
The production of the SM Higgs boson corresponds to the case where c 2HDM H is equal to one. In particular, c 2HDM H l of Eq. (23) becomes one, i.e. SM-like, in the alignment limit given in Eq. (8). One can also find a combination of c αβ and t β such that c 2HDM H becomes arbitrary small or even zero. In this case, the two-loop, electroweak corrected process becomes the true leading order. The symbol M q is the mass of the internal quark, which in general can be either the up (u), down (d), strange (s), charm (c), bottom (b) or top (t) quark. We consider only the top quark as massive and all other fermions as massless, i.e. only the case q = t contributes to Eq. (20).
There are no real corrections which have to be determined when computing the NLO electroweak corrections for this process. Consequently, one can straightforwardly apply the results for the electroweak percentage correction, which has been obtained for the cross section of the production process gg → H in Eq. (20), to the partial decay width of the process H → gg. The partial decay width and the production cross section are related by The Feynman rules for the 2HDM of type II, which has been defined in Section 2, have been produced in an automated way with the help of FeynRules [29]. The one-and twoloop diagrams have been calculated with QGS which is an extension of GraphShot [30]. It generates the Feynman diagrams with QGRAF [31] and performs algebraic manipulations of the amplitudes and loop integrals with the help of Form [32]. The program QGS performs the standard Dirac algebra operations and projects the expressions onto form factors. In addition, it removes reducible scalar products and uses integration-by-parts identities in order to simplify tadpole contributions. Finally, it reduces the number of loop integrals by means of symmetrization. Further details regarding these techniques can be found in Ref. [15]. After these manipulations, one remains with a two-loop amplitude which requires numerical evaluation, which is done with a Fortran code. In particular, the two-loop integrals are evaluated in Feynman-parametric space with the help of a Fortran library. For the numerical integral evaluation of the two-loop massive diagrams, the library uses the methods of Ref. [33] for self-energies and of Refs. [34,35,15] for vertex functions.
Compared to the production of a light, neutral, SM-like Higgs boson H l in gluon fusion in the alignment limit, which has been addressed already briefly in Ref. [11], the non-alignment limit case and the production of a heavy, neutral Higgs boson H h is computationally more involved. While the calculation of the integrals for the production of the SM-like, light, neutral Higgs boson H l in the alignment limit can be completely traced back to the calculation of the pure SM case, which was discussed in Refs. [14,15], new integral structures appear in the general case. The first diagram in Fig. 2 is non-planar and leads to a new rank 2 integral. In the SM and in the alignment limit, this diagram only exists with two Z-bosons. The same type of diagram exists for the W-bosons and the charged Higgs bosons (second diagram of Fig. 2). Since different types of fermions appear in this diagram, it leads to rank 3 integrals, and hence, the calculation becomes even more involved. In the SM, it was possible to cancel the rank 3 integrals since there were fewer different bosonic masses in the denominators. A new planar diagram, which appears in the general case of the 2HDM, is depicted on the right hand side of Fig. 2. After performing the obvious reduction as explained in Ref. [15], the second and third diagram of Fig. 2 lead to the collinear structures present in the first diagram of Fig. 3. The non-planar diagram also contains the more complicated collinear structure in the second diagram of Fig. 3. More details regarding the cancellation of these new collinear singularities are given in Section 4.1.
In Section 5, we provide numerical results for the NLO, electroweak corrections to the production processes gg → H l and gg → H h for various 2HDM scenarios. We define the NLO corrections in terms of the K-factor using the one-loop A H and two-loop A H contributions to the amplitude with the normalization defined in Eq. (20). Figure 2: Sample diagrams, which do not appear in light Higgs-boson production in the alignment limit, are shown. In the second and third diagram, the bottom quarks are considered to be massless. For light Higgs-boson production the contribution to the amplitude originating from these three diagrams is proportional to cos (α − β). This is the reason why they vanish in the alignment limit. H is very small too, and the two-loop amplitude A (2) H becomes the true leading order result. As discussed in the context of Higgs-boson production and decay with a fourth fermion generation [10], neglecting the term |A (2) H | 2 is no longer justified in this case, and we define a new K-factor Whenever K NLO EW differs significantly from K NLO EW , K NLO EW will be used in Section 5.

Analytical and numerical tests
We have performed several analytical as well as numerical tests to validate the new components of QGS, which are the 2HDM Feynman rules, automatic generation of Feynman diagrams and the appearances of new rank 2 and 3 integrals. The ultraviolet (UV) structure of the new integrals as well as the consistency of the 2HDM Feynman rules has been tested by extracting all UV poles in dimensional regularization, as explain in Ref. [15], and verifying their cancellation analytically. Secondly, the cancellation of collinear singularities can also be used to validate the implementation of the Feynman rules and the calculation of the new integral structures. Collinear singularities arise in some diagrams, if the external gluons couple to light fermions, but they have to cancel when summing up the contributions from all diagrams. These singularities have been regularized by fictitious small fermion masses and become manifest in terms of linear and quadratic logarithms. The two-loop electroweak amplitude is subdivided into the contribution which comes from the first and second generation of fermions and the contribution which originates from the third generation. The cancellation of the collinear divergent logarithms, which arise from the first and second generation of fermions, has been verified analytically. For the third generation of fermions, two new collinear structures, originating from the second and third diagram in Fig. 2, are presented in Fig. 3. Their analytic cancellation requires new integration-by-parts identities and more complicated partial-fraction decomposition due to new denominators. All quadratically collinear divergent logarithms can then be cancelled analytically. The linearly collinear divergent logarithms have been treated in a numerical approach. Here, the logarithm is extracted analytically, but its coefficient is evaluated numerically. We have verified that the sum of all these coefficients cancels numerically.
Another way to verify our implementation of the 2HDM is to test the behaviour of the NLO electroweak corrections in different limits. In the decoupling limit introduced in the end of Section 2, all new Higgs-boson masses are much larger than the electroweak scale, such that the new particles decouple from the SM. Therefore, the NLO corrections in the 2HDM should approach those of the SM. Fig. 4 shows a decoupling scenario for gg → H l . The 2HDM parameters have been chosen as where all masses are set equal to the scale of new physics M * , because perturbativity requires the mass splitting to be smaller than v 2 /M * [19]. At M * = 2400 GeV, the elec-   Table 1: 2HDM benchmark points (BP) in the alignment limit, i.e. c αβ = 0, taken from Ref. [7]. In the alignment limit, |c 2HDM  Table 2: 2HDM benchmark points (BP) outside the alignment limit taken from Ref. [36] (a-1, b-1) and Ref. [7]. account: since the LO of its defining process H h → ZZ vanishes in the alignment limit, it cannot be expected to show proper decoupling behaviour. Furthermore, several one-loop processes and all finite NLO counterterms have been cross-checked against the implementation in RECOLA2 [12] to verify the implementation of the Feynman rules and the generation of one-loop diagrams.

Results and discussion
In this section, we present the numerical results of phenomenologically interesting scenarios for light and heavy Higgs-boson production in gluon fusion. First, we consider the benchmark points (BP) in Tabs. 1 and 2 in different renormalization schemes. The BPs a-1 and b-1 are taken from Ref. [36]. All other BPs are from the LHC Higgs cross section working group report [7]. Then, in order to analyze the influence of the mass splitting on the perturbative behaviour of the Higgs boson production processes, we move to a benchmark scenario where we keep all new Higgs boson masses fixed at 700 GeV and vary the heavy Higgs boson mass between 600 GeV and 800 GeV. This scenario agrees with current experimental and theoretical exclusion limits [4,37,38].
For the numerical evaluation we use the following set of SM input parameters [2]: For the renormalization of the W -and Z-boson masses we use the complex mass scheme [39].
In contrast to Ref. [11], where the top-quark mass was renormalized on-shell, here we also use the complex mass scheme [39] for the top-quark mass.

Benchmark points
According to the analysis of the LHC Higgs cross section working group [7], the BPs in Tabs. 1 and 2 fulfill constraints like perturbativity and vacuum stability. In Tabs Tab. 3 contains the NLO electroweak K-factors for light Higgs-boson production in gluon fusion in different renormalization schemes in the alignment limit. For this process, we provide only the factor K NLO EW , because the difference between K NLO EW and K NLO EW is tiny -at the per mille level. The OS, p * and the two process-dependent renormalization schemes produce similar results and mostly, the K-factors are comparable in size with the SM value K SM EW = 1.051. For the MS-renormalized values, we have chosen µ 0 according to Eq. (13), different from the choice of the renormalization scale in Ref. [11]. The values for c αβ and t β of the benchmark points in Tabs. 1 and 2 are defined at the scale µ 0 . The uncertainty is obtained as described in Section 3 by taking one half and twice the central value µ 0 and running the parameters c αβ and t β to the corresponding scales.    For the remaining benchmark points with c αβ = 0, the K-factors of the electroweak corrections to the process gg → H l are shown likewise in Tab. 4. Again, the K-factors of the OS, p * and the two process-dependent renormalization schemes are close to the SM, while the MS scheme leads to different results. Nevertheless, the electroweak corrections are always moderate in size, i.e. they are mostly below 5% compared to the LO production cross section.
We now turn to heavy, neutral Higgs-boson production gg → H h . There are two aspects which we have to consider when looking at this process. On the one hand, the LO production cross sections in Tabs. 5 and 6 depend on the heavy Higgs-boson mass as well as on the coefficient |c 2HDM On the other hand, the coefficients |c 2HDM H h | 2 can be close to zero. This is different from the light Higgs-boson production case, where for all BPs the same light Higgs-boson mass enters, and where the coefficients |c 2HDM H l | 2 , as given in Tab. 2, are always close to 1 or even equal to 1, as it is the case in the alignment limit. For these reasons, the LO cross sections of the heavy Higgs-boson production can be very small and in particular comparable in size to the NLO contribution. Then, further terms of the perturbative expansion of the cross section are required where A H h is the three-loop amplitude. In Eq. (29), |A H h | 2 is the LO cross section up to a global factor. Likewise, A H h enter at NNLO. If |A (1) H h | 2 is very small, |A (2) H h | 2 can be of the same size as the NLO contribution, while A H h can be expected to be small again. Therefore, |A (2) H h | 2 should not be neglected. Furthermore, the imaginary part of the LO amplitude can be similar in magnitude or even larger than the real part when the heavy Higgs-boson mass becomes larger than twice the top-quark mass. This can lead to cancellations between real and imaginary parts of the NLO term.
In both cases, there can be a considerable difference between K NLO EW and K NLO EW as defined in Eqs. (25) and (26). Since K NLO EW contains the additional contribution, it seems to be more appropriate to quantify the electroweak corrections to heavy Higgs-boson production. In order to point out BPs where this difference occurs, the K NLO EW are shown in parentheses in Tabs. 5 and 6 even though they can take an unphysical, negative value for some benchmark points, e.g. for 2 1B and 3 B1 .
Different from the sister process gg → H l , many K-factors in Tabs. 5 and 6 are far from one. This is at least partly due to the different LO couplings between the t − H l − t and t − H h − t interactions. While the |c 2HDM H l | 2 in Tabs. 1 and 2 are generally one or close to one, the |c 2HDM H h | 2 are at most 0.57. In addition, the decoupling of the new Higgs sector from the SM can play a role for BP a-1 and 2 2A . While decoupling imposes the restriction that SM processes may not receive large corrections from the new sector, this is not the case for non-SM processes like heavy Higgs-boson production.
In the following, we have a closer look at the benchmark points 2 1B , 2 2A and 3 B1 , which have large electroweak NLO corrections and a large difference between K NLO EW and   K NLO EW for heavy Higgs-boson production. Benchmark point 2 1B does not only show very large corrections of more than -80%, but also the difference between K NLO EW and K NLO EW is of the same order of magnitude. Hence, |A H h | and |A (2) H h | must be similar in size, which may partially be due to the large value of |λ max i |/(4π) = 0.57. The same is true for benchmark point 2 2A . It has a large value of |λ max i |/(4π) = 0.64, and in addition, its LO cross section is quite small, such that the two-loop diagrams yield the true LO contribution. In order to verify whether the NNLO contribution is small in these scenarios, would, however, require the calculation of the three-loop contributions, which is not within reach in the near future.
The benchmark point 3 B1 has an even smaller LO cross section due to the almost vanishing coefficient |c 2HDM H h | 2 = 0.0003. Therefore, the extremely large K-factor of more than 27 can be understood, since the true LO contribution is given by the two-loop diagrams. However, since the coefficients in front of the counterterms of the mixing angles are large for this BP and the small LO does not factorize, we observe large differences between the results in different renormalization schemes. This leads to a large dependence on the choice of the renormalization scheme, and thus to large theoretical uncertainties on the production cross section. In particular, this shows that a small coupling of the heavy Higgs boson to top quarks does not automatically lead to a small cross section when higher order contributions are considered.
Finally, there are benchmark points which have moderate NLO corrections on the one hand, but a relatively large difference between K NLO EW and K NLO EW on the other hand. We have investigated this difference further and found that it is caused by cancellations between the real and imaginary contributions in the term A H h . Especially benchmark point a-1 exhibits this feature, but it can also be observed in BP 3 B2 and 4 4 .
As a measure of the scale dependence in the MS scheme, we can look at the size of the coefficients of the scale-dependent logarithms in the percentage correction presented in Tabs. 7 and 8 of Appendix A. In general, they can become very small or even vanish for certain choices of the parameters. In particular, in the alignment limit with M sb = 0, Eq. (32) becomes very simple  In the MS scheme the corrections are renormalization-scale dependent and we choose the central renormalization scale µ 0 according to Eq. (13). In order to estimate the impact of the scale dependence on the size of the NLO EW corrections, we vary the scale between µ 0 /2 and 2µ 0 . A more detailed analysis is discussed in Appendix C. In the alignment limit presented in the first plot of Fig. 6, the different renormalization schemes agree quite well for |λ i |/(4π) < 0.5. The curves for MS renormalization (blue) and the renormalization scheme proc1 (purple, dashed-dotted) show the biggest deviation from the SM, especially when entering the nonperturbative region. The corrections of the OS, p * and proc2 renormalization schemes are very close to each other for these scenarios, even for large λ i . For c αβ = 0.03 presented in the second plot of Fig. 6, the behaviour of the OS, p * , proc1 and proc2 scheme is very similar to the case of the alignment limit (upper plot), while the MS results show a much larger scale variation. Fig. 7 shows the K-factor of the NLO electroweak corrections in different renormalization schemes for the process gg → H h in the 2HDM as a function of the heavy, neutral    Higgs-boson mass for the same scenarios as in Fig. 6. Again, as for the benchmark points, we use K NLO EW as the K-factor for heavy Higgs-boson production. First of all, we can expect the NLO corrections for heavy Higgs-boson production to be larger: The LO contribution is suppressed since for t β = 2 and close to the alignment limit, the coefficient c 2HDM H h ≈ −0.5 is considerably smaller than c 2HDM H l ≈ 1. In addition, the LO depends on the heavy Higgs-boson mass M H h , which is not fixed in the scenario under consideration. In general, we can see that for |λ max i |/(4π) < 0.5, the K-factors do no longer lie roughly between 0.99 and 1.12 as for light Higgs-boson production, but now we observe a larger range of K-factors between 0.4 and values larger than 1.8. Just as for light Higgs-boson production, the proc1 scheme differs considerably from the other renormalization schemes for large M H h . The behaviour of the results for the MS renormalization seems to be more similar to the other schemes when compared with the light Higgs-boson production, even though there are sizable numerical differences in some regions. However, this scheme strongly depends on the definition of the running procedure requiring a more detailed analysis, which is performed in Appendix C.
As expected, due to cancellations among the finite counterterms, the OS and p * scheme lead to small perturbative corrections. The differences between these two schemes may be too small to give an estimate of missing higher-order uncertainties. As long as M H h does not become too large the two process-dependent schemes also perform very well for the analyzed scenarios, both for light as well as for heavy Higgs-boson production.

Summary and conclusion
We have computed the two-loop electroweak corrections to the production of a light and a heavy neutral, scalar Higgs boson through gluon fusion in the 2HDM. We have renormalized the new Higgs-boson masses in the on-shell scheme and provide the electroweak percentage correction in different renormalization schemes for the mixing angles α and β for these two processes. In particular, for the mixing angles we have employed the on-shell, p * , MS and two process dependent schemes. We can determine the next-toleading order, electroweak percentage correction for essentially any scenario of the new mass parameters and the mixing angles of the 2HDM. In particular, we have computed the two-loop electroweak corrections for benchmark points collected by the LHC Higgs cross section working group as well as for individual other example scenarios. For Higgsboson production through gluon fusion, the on-shell scheme performs well for all chosen scenarios. The MS scheme produces reasonable results only for certain choices of the new input parameters. For the production of the light Higgs-boson, the electroweak corrections are always moderate in size, i.e. they are mostly around 5% compared to the LO production cross section for the chosen benchmark points, while for the production of a heavy Higgs-boson, the electroweak corrections strongly vary depending on the details of the selected scenario. We have solved new technical challenges which was required to accomplish this calculation. Our results are also directly applicable to determine the electroweak percentage corrections for the partial decay widths of the light and heavy  Table 7: The coefficient in front of the scale dependent logarithm of the percentage correction is shown for the benchmark points which are in the alignment limit (c αβ = 0).  Table 8: The coefficient in front of the scale dependent logarithm of the percentage correction is shown for the benchmark points which are not in the alignment limit (c αβ = 0).
In particular for the BPs in the alignment limit, the coefficient can easily be obtained by Eqs. (32) and (33). The explicit values of the coefficients are shown in Tabs. 7 and 8. The size of the coefficients reflects directly the magnitude of the scale dependence in the MS scheme without taking into account the running of the mixing angles c αβ and t β . Coefficients whose absolute value is larger than one, like for example for the BPs 2 2A , 3 B1 and 3 B2 for the process gg → H h , exhibit a large scale dependence and give rise to large corrections. Small coefficients, like for example for the BPs 2 1C , 2 1D , 2 2A and 3 B2 for the process gg → H l , exhibit a very small scale dependence. Fig. 8 shows the range of the parameters |λ i |/(4π) of the potential for the alternative scenarios discussed in Section 5, where we have fixed the values of the mixing angles to t β = 2 and c αβ = 0 or c αβ = 0.03. All new heavy mass scales are set to the same value of 700 GeV, except for one Higgs-boson mass, which is varied between 600 and 800 GeV. In the two plots of Fig. 8 the parameters |λ i |/(4π) are shown as a function of M Ha and M Hc , while the plot as a function of M H h is given in Fig. 5 of Section 5. The size of the couplings is less sensitive to a variation of M Ha and M Hc than to a variation of M H h as shown in Fig. 5.
see also Eq. (31). In particular in Figs. 10 and 11 we compare the scale dependence of the EW corrections in three cases: • The dependence on the renormalization scale µ enters just in the loop corrections. The running of the mixing angles α and β is not taken into account (red bands).
• The dependence on the renormalization scale µ enters in the loop corrections as well as in the running of the mixing angles α and β, see Eq. (14); the default scale is µ d = µ 0 , where the soft-breaking scale takes the value M sb = 700 GeV, and where the mixing angles t β and c αβ take the benchmark values t β = 2, and c αβ = 0 or c αβ = 0.03 (blue bands).
• The dependence on the renormalization scale µ enters in the loop corrections as well as in the running of the mixing angles α and β, see Eq. (14); the default scale µ d at which the soft-breaking scale takes the value M sb = 700 GeV and at which the mixing angles t β and c αβ take the benchmark values t β = 2 and c αβ = 0 or c αβ = 0.03 is µ d = vev (yellow bands).
The bands in Figs. 10 and 11 are generated by considering the region of the parametric space between the three curves corresponding to the values µ 0 , µ 0 /2 and 2µ 0 for the renormalization scale, where µ 0 is the typical scale of the processes under consideration as calculated in Eq. (13). The curve of the on-shell renormalization scheme (dashed, black curve) is shown as a reference curve for the scale-independent schemes. For an MS renormalization of the mixing angles, the couplings λ i of Eq. (2) are scale dependent. The grey shaded band below each figure shows the region, where at least one of the MS couplings |λ i |/(4π) of Eq. (2) for at least one of the values of µ becomes larger than 0.5 in the considered renormalization scenarios. Due to the scale dependence, the grey shaded band for an MS renormalization is in general different from the case of considering a scale-independent scheme for the mixing angles, as shown in Figs. 6 and 7. Therefore, the regions where one enters in the non-perturbative regime can also be different.
The K-factor for the process gg → H l is shown in Fig. 10 for c αβ = 0 (upper plot) and c αβ = 0.03 (lower plot). In the alignment limit, c αβ = 0, the corrections do not differ much from the SM in all considered MS schemes when the couplings |λ i |/(4π) remain below 0.5 (white region); they are just slightly smaller (1.03 < K NLO EW < 1.055). The situation changes drastically in the non-perturbative grey region where scale variation becomes large and seems to reveal a certain correlation between perturbativity and scale dependence for this process in the alignment limit. This correlation, however, can be a peculiar property of the chosen scenarios: from Eq. (32) we see that in the case of a heavy new scalar sector, the largest contribution to the scale dependence in the MS scheme comes from the term proportional to (M 2 H h − M 2 sb ) when the running of the mixing angles α and β as well as of the soft-breaking scale M sb is not taken into account; the couplings λ i of Eqs. (35)-(39) on the contrary, show a more complicated dependence on the masses of the neutral sector. In the case of exactly equal heavy masses in the alignment limit, they receive the largest contribution from the (M 2 H h −M 2 sb ) difference in λ 1 . Moving slightly away from the alignment scenario (lower plot with c αβ = 0.03), we notice a larger difference from the SM value and a wider scale dependence, even in the white (perturbative) region. In particular, for M H h above 700 GeV, even remaining in a region where all |λ i |/(4π) are below 0.5, the MS bands become quickly wider. Due to this behaviour, MS scale dependence can not serve as a good estimator of the theoretical uncertainty in the light, neutral Higgs-boson production away from the alignment limit.
In Fig. 11, we compare the different MS schemes for the production of a heavy, neutral Higgs-boson in the alignment limit (upper plot) and slightly away from it (lower plot). In this case we do not have a comparison with the SM and we base our considerations on a comparison with the curve of the on-shell renormalization scheme. Restricting the analysis of the case c αβ = 0 to the white (perturbative) region we notice a strong dependence on M H h for all results. In addition they change a lot for the various MS schemes and for different values of the scale µ. The situation is similar but less dramatic for c αβ = 0.03. Owing to these large differences the MS scheme can not reliably predict the NLO contributions to heavy, neutral Higgs-boson production.