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,2] with a mass of M h = 125.09 ± 0.21 ± 0.11 GeV [3], 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. Such an extension was introduced in ref. [4] to provide an additional source of CP-violation, which may contribute to explain the observed matter-anti-matter asymmetry of the Universe. Special choices of the parameters, for example in the Inert Model, can also provide a dark matter candidate [5,6]. Thus, the 2HDM can assist to address problems which are not solved within the SM. Phenomenological studies of the 2HDM have been performed by the ATLAS [7][8][9][10][11][12][13][14] and CMS [15][16][17][18][19][20][21] 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. [22].
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 [23]. The next-to-leading order (NLO) electroweak corrections in this model have been computed in refs. [24,25], which have helped to exclude this model at the LHC.

JHEP09(2018)017
In the CP-conserving 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 H a . A possible set of new free parameters of the extended Higgs sector are the masses of the new Higgs bosons, M H h , M H ± , M Ha as well as the soft Z 2 -breaking scale M sb and two mixing angles α and β. The mass of the light Higgs boson is here fixed to M H l ≡ M h .
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. [26,27]. The NLO electroweak corrections to the decay of the light Higgs boson of the 2HDM into four fermions have been computed in refs. [28,29]. 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 [3] amount to 5.1% [30,31]. 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. [26]. 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 also discuss 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

JHEP09(2018)017
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 [32], which does not induce the FCNC problem. With these constraints, the Higgs potential has the following form (2. 2) 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 that 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 bosons, 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 that 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. (B.1)-(B.5) of appendix B. Finally, we also introduce the Higgs basis [33] 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 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. (2.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. [26] 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 that 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 [34,35]. For the type II 2HDM, down-type quarks and charged leptons need to be odd under this Z 2 transformation, i.e. Φ 1 → −Φ 1 , d R → −d R , R → − R , 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 and the mixing angles α and β, which we express through

JHEP09(2018)017
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 [22,36].
3 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. [26]. 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. [37] and which was extended for the 2HDM, and also for a general Higgs sector, in ref. [38] and ref. [26]. In the SM, the FJ tadpole scheme corresponds to the β t scheme of ref. [39]. 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. [26]. For an MS renormalization of M sb see also ref. [40]. At next-to-leading order, the counterterms δα and δβ of the 2HDM can be obtained from the off-diagonal elements of the field-renormalization constants of the Higgs-and Goldstoneboson fields [41]. In particular, for a gauge-independent MS renormalization, this relation reads where the Z-factors are defined by 3)

JHEP09(2018)017
where the subscript B denotes a bare field. For the renormalization conditions of the Z-factors we refer to ref. [26]. The choice of an MS renormalization leads to scale-dependent renormalized amplitudes. In order to account for the different scales in the new Higgs sector, we use the average as a central renormalization scale, see refs. [28,29] 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. (3.4) 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) (3.5). 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. (2.4), see refs. [42,43]. However, using finite parts of the mixing two-point functions to define the renormalization constants of α and β leads to gauge-dependent renormalization constants in addition to the gauge-dependent tadpole terms discussed in ref. [26]. The problem has been overcome in ref. [38] using the pinch technique [42][43][44][45][46][47][48][49][50][51][52][53][54][55]. Recently, it has been treated in a more general framework in ref. [27] within the background field method (BFM), leading to the same result. In the JHEP09(2018)017 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 [38] Σ 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 [56,57].
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. [26]. 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. [27].
p * renormalization scheme. The p * scheme is derived from the same mixing two-point functions of eq. (3.6), but computed at the external momentum such that the rotation matrices in eq. (2.4) also remain scale-independent, see refs. [42,43]. The additional terms in eqs.  of the mixing angles (3.10) 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 → τ + τ − [38]. 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 figure 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 JHEP09(2018)017 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. [38]. 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 g + g → H l and g + g → 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 figure 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 δβ, contrary to 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 A LO reads where the index q runs over all quark flavours and in eq. (4.2) originates from the Higgs-boson-quark coupling and is in the 2HDM, type II different for up-and down-type quarks and also different for the production of a light or a heavy, neutral Higgs boson: The production of the SM Higgs boson corresponds to the case where c 2HDM H,q is equal to one. We consider only the top quark as massive and all other fermions as massless, i.e. only the case q = t contributes to eq. (4.1). In this case only c 2HDM of eq. (4.4) contributes to eq. (4.2). For simplicity we will drop the label t in the following, i.e. c 2HDM H l of eq. (4.4) becomes one, i.e. SM-like, in the alignment limit given in eq. (2.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 process.
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 g + g → H in eq. (4.1), to the partial decay width of the process H → g +g. 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 [58]. The one-and two-loop diagrams have been calculated with QGS, which is an extension of GraphShot [59]. It generates the Feynman diagrams with QGRAF [60] and performs algebraic manipulations of the amplitudes and loop integrals with the help of Form [61,62]. 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. [31]. After these manipulations, one remains with a two-loop amplitude that JHEP09(2018)017 Figure 2. Sample diagrams that 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. requires numerical evaluation, which is done with a Fortran code. In particular, the twoloop 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 refs. [63,64] for self-energies and of refs. [31,[65][66][67] 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. [26], 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 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. [30,31], new integral structures appear in the general case. They are generated by new diagrams, a sample of which is shown in figure 2. The first diagram 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 figure 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 at the rightmost side of figure 2. After performing the obvious reduction as explained in ref. [31], some diagrams lead to collinear divergent structures; the new ones, not present in the alignment limit, are shown in figure 3. In particular, the second non-planar diagram of figure 2 leads to the more complicated collinear structure depicted in the second diagram of figure 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 g + g → H l and g + g → 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. (4.1). 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 [25], neglecting the term |A (2) H | 2 is no longer justified in this case, and we define a new K-factor

JHEP09(2018)017
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 explained in ref. [31], and verifying their cancellation analytically.
Furthermore, 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 that comes from the first and second generation of fermions and the contribution that 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 figure 2, are presented in figure 3. Their analytic cancellation requires new integration-by-parts identities and a 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. Figure 4 shows a decoupling scenario for g + g → H l . The 2HDM parameters have been chosen as

JHEP09(2018)017
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 * [36]. At M * = 2400 GeV, the electroweak corrections in the 2HDM have almost approached those of the SM. In addition, we observe that the scale dependence of the MS-renormalized corrections shrinks, cf. appendix A for the analytic formula. The proc2 renormalization has not been taken into 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 [27] 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 for the CP-conserving 2HDM. First, we consider the benchmark points (BP) in tables 1 and 2 in different renormalization schemes. The BPs a-1 and b-1 are taken from ref. [68] and they correspond to best-fit constraints on the triple Higgs couplings from the Higgs signal strengths. All other BPs are from the LHC Higgs cross section working group report [22] and are sample points of phenomenologically interesting allowed scenarios compatible with the 2HDM type II.  Table 1. 2HDM benchmark points (BP) in the alignment limit, i.e. c αβ = 0, taken from ref. [22]. In the alignment limit, |c 2HDM  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 except for the heavy Higgs boson mass which is varied between 600 GeV and 800 GeV. This scenario agrees with current experimental and theoretical exclusion limits [7,69,70].
For the numerical evaluation we use the following set of SM input parameters [3]: For the renormalization of the W -and Z-boson masses we use the complex mass scheme [71][72][73]. In contrast to ref. [26], where the top-quark mass was renormalized onshell, here we also use the complex mass scheme [71][72][73] for the top-quark mass.
Benchmark points. According to the analysis of the LHC Higgs cross section working group [22], the BPs in tables 1 and 2 fulfill constraints like perturbativity and vacuum stability. In tables 1 and 2 we also provide the maximum of the couplings |λ max i |/(4π) of the Higgs potential of eq. (2.2). Compared to the SM value λ SM /(4π) = 0.02, all |λ max i |/(4π) values are large. This can potentially lead to large NLO corrections.  Table 3. The NLO electroweak K-factors for the process g + g → H l for the 2HDM benchmark points (BP) in the alignment limit are shown. They are defined in table 1. The renormalization scale has been set to µ 0 as defined in eq. (3.4) for the MS scheme. The lower (upper) value of the MS result corresponds to the change when taking µ 0 /2 (2µ 0 ) as renormalization scale. Table 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. (3.4), different from the choice of the renormalization scale in ref. [26]. The values for c αβ and t β of the benchmark points in tables 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. The results in the MS scheme are different from the other results and the benchmark points 2 1A and 2 1B show a large scale dependence as already analyzed in ref. [26], where no running of the parameters was considered yet. The running results in an enhanced scale dependence for these two benchmark points.

JHEP09(2018)017
For the remaining benchmark points with c αβ = 0, the K-factors of the electroweak corrections to the process g + g → H l are shown likewise in table 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 g + g → H h . There are two aspects that we have to consider when looking at this process. On the one hand, the LO production cross sections in tables 5 and 6 depend on the heavy Higgs-boson mass as well as on the coefficient |c 2HDM   Table 4. The NLO electroweak K-factors for the process g + g → H l for the 2HDM benchmark points (BP) that are not in the alignment limit are presented. They are defined in table 2. The renormalization scale has been set to µ 0 as defined in eq.
H h is the three-loop amplitude. In eq. (5.2), |A H h | 2 is the LO cross section up to a global factor. Likewise, A H h corresponds to the NLO contribution, while |A 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 contributions in the NLO term.
In both cases, there can be a considerable difference between K NLO EW and K NLO EW as defined in eqs. (4.6) and (4.7). 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 tables 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 g + g → H l , many K-factors in tables 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 tables 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 JHEP09(2018)017  Table 6. The NLO electroweak K-factors for the process g + g → H h for the 2HDM benchmark points (BP) that are not in the alignment limit are shown. They are defined in table 2. The renormalization scale has been set to µ 0 as defined in eq. 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 that 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 JHEP09(2018)017 in tables 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, the scale dependent logarithmic term for g + g → H l can be obtained from eq. (A.1) of appendix A and becomes very simple which even vanishes for the special case M 2 H h = 2m 2 t /(t 2 β − 1). This scenario is almost realized for BP 2 1C and 2 1D , and hence, these benchmark points exhibit only a small scale dependence in table 3.
In addition, we observe that the MS corrections for light Higgs-boson production are usually moderate, while for heavy Higgs-boson production, this scheme often leads to Kfactors far from one in addition to a very strong scale dependence. This can, again, partly be traced back to small LO couplings between the heavy Higgs boson and the top quarks. The first scenario uses the alignment limit and the second scenario sets c αβ = 0.03. As we see in figure 5, both scenarios fulfill the perturbativity restriction |λ i |/(4π) < 1 (i = 1, . . . , 5). In addition, they respect the experimental constraints considered in refs. [7,69].
In both scenarios we vary the heavy, scalar Higgs-boson mass between 600 and 800 GeV. The size of the couplings λ i is more sensitive to these variations compared to variations in the other Higgs-boson masses M Ha , M Hc (compare figure 5 to figure 8 in appendix B). For that reason we will present the NLO corrections to Higgs-boson production and also the scale dependence in the MS renormalization scheme as a function of the heavy Higgs-boson mass in the following. Figure 6 shows the K-factor of the NLO electroweak corrections for the process g +g → H l in the 2HDM as a function of the heavy, neutral Higgs-boson mass for the different renormalization schemes of the mixing angles α and β. We consider the OS, p * , two processdependent and the MS scheme. The grey shaded band shows the region where at least one of the couplings |λ i |/(4π) of eq. (2.2) becomes larger than 0.5. The region of variation of the corresponding |λ i |/(4π) values is displayed in figure 5. These values are not valid for the MS scheme, where they adopt different values depending on the choice of the renormalization scale. As a result of this, the region of perturbativity generally looks different for the MS renormalized results. In the MS scheme the corrections are renormalization-scale dependent and we choose the central renormalization scale µ 0 according to eq. (3.4). 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 figure 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 non-perturbative 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 figure 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. Figure 7 shows the K-factor of the NLO electroweak corrections in different renormalization schemes for the process g + g → H h in the 2HDM as a function of the heavy, neutral Higgs-boson mass for the same scenarios as in figure 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 Higgsboson 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 JHEP09(2018)017

JHEP09(2018)017
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-to-leading order, electroweak percentage correction for essentially any scenario of the new mass parameters and the mixing angles of the CP-conserving 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 can suffer from a large scale dependence and can in general not provide reliable predictions for all scenarios. 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 neutral, scalar Higgs-boson decay into two gluons within the 2HDM.  Table 7. The coefficient in front of the scale dependent logarithm of the percentage correction is shown for the benchmark points that are in the alignment limit (c αβ = 0). • 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).

JHEP09(2018)017
• The dependence on the renormalization scale µ enters in the loop corrections as well as in the running of the mixing angles α and β, see eq. (3.5); 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. (3.5); 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 figures 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. (3.4). 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 JHEP09(2018)017 Figure 10. The K-factors of the process g + g → H l in various MS schemes are compared to the on-shell scheme and to the SM result. Solid lines are for µ = µ 0 , while dotted and dash-dotted lines refer to µ = µ 0 /2 and µ = 2µ 0 , respectively. The grey shaded band denotes the region where at least one of the couplings |λ i |/(4π) (i = 1, . . . 5) of the potential in eq. (2.2) becomes larger than 0.5 for at least one of the scales µ considered. The curve with µ = µ 0 for the case of the MS scheme without running coincides with the µ = µ 0 curve for the MS scheme with running, where the parameters are defined at µ d = µ 0 , since in both cases the parameters have been defined at µ 0 . renormalization of the mixing angles, the couplings λ i of eq. (2.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.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 JHEP09(2018)017 Figure 11. The K-factors of the process g + g → H h in various MS schemes are compared to the on-shell scheme. Solid lines are for µ = µ 0 , while dotted and dash-dotted lines refer to µ = µ 0 /2 and µ = 2µ 0 , respectively. The grey shaded band denotes the region where at least one of the couplings |λ i |/(4π) (i = 1, . . . 5) of the potential in eq. (2.2) becomes larger than 0.5 for at least one of the scales µ considered. The curve with µ = µ 0 for the case of the MS scheme without running coincides with the µ = µ 0 curve for the MS scheme with running, where the parameters are defined at µ d = µ 0 , since in both cases the parameters have been defined at µ 0 .

JHEP09(2018)017
band for an MS renormalization is in general different from the case of considering a scaleindependent scheme for the mixing angles, as shown in figures 6 and 7. Therefore, the regions where one enters in the non-perturbative regime can also be different.
The K-factor for the process g + g → H l is shown in figure 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. (A.1) 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. (B.1)-(B.5) 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 figure 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.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.