One-loop radiative corrections to $e^+ e^-\to Zh^0/H^0A^0$ in the Inert Higgs Doublet Model

We compute the full one-loop radiative corrections (including both weak and QED corrections) for two processes $e^{+}e^{-}\to Z h^0,H^0 A^{0}$ in the Inert Higgs Doublet model (IHDM). Up to $O(\alpha_{w})$ and $O(\alpha_{em})$ order, we use FeynArts/FormCalc to compute the one-loop virtual corrections and Feynman Diagram Calculation (FDC) to evaluate the real emission, respectively. Being equipped with these computing tools, we investigate radiative corrections of new physics for both the degenerate and non-degenerate scenarios with three typical collision energies of future electron-positron colliders: 250 GeV, 500 GeV, and 1000GeV. By scanning the parameter space of IHDM, we identify the allowed regions which are consistent with constraints and bounds, from both theoretical and experimental sides. We find that the radiative corrections of the IHDM to $e^+ e^- \to Z h^0$ can be sizeable and are within the detection potentials of future Higgs factories. We also find that the new physics of IHDM could also be directly detected by observing the process $e^{+}e^{-}\to H^0 A^{0} $ which could have large enough production rate. We propose five benchmark points and examine their salient features which can serve as physics targets for future electron-positron colliders, such as CEPC/CLIC/FCC-ee/ILC as well as for LHC.


Introduction
The first LHC run with 7⊕8 TeV and the second one with 13 TeV were successful operations which led to the discovery of a new scalar particle [1,2] including the recent observation of its production in association with tt [3,4] and its decay to bb [5,6] among other achievements. The LHC program has already performed several precise measurements in term of production cross sections and branching fractions. These measurements demonstrate that the Standard Model (SM) works well to explain these observed phenomena at the electroweak scale.
One of the main goals of the future run of the LHC with 14 TeV and its High Luminosity option (HL-LHC) is to improve the aforementioned measurements and pin down the uncertainties to few percent level [7][8][9][10][11]. On the other hand, it is also expected from the future LHC run to establish a new measurement such as the triple Higgs coupling [12], and Higgs decay into γZ and µ + µ − . Moreover, it is well known that a precise measurement program which already began at the LHC, is expected to be performed at the e + e − machines [13][14][15] such as the Circular Electron Positron Collider (CEPC) [16,17], the Compact Linear Collider (CLIC) [18][19][20][21][22], the Future Circular Collider (FCC-ee) [23,24], and the International Linear Collider (ILC) [13,25]. The e + e − machines which are expected to deliver rather high luminosity and possess a very clean environment, would be able to improve the Higgs couplings and production cross section measurements below the percent level [13][14][15]. Such a precision, if achieved, will be very useful to discover the evidence of new physics beyond the SM.
In Table 1.1, we provide the projected experimental precisions for δg ZZh 0 at future Higgs factories, which is taken from Table 4.2 of reference [24]. It should be declared that the numbers quoted here are just for illustration on the potential precision which could be reached, but they are not necessarily exact what those Higgs factories are planned to achieve since the actual precisions depend on more detailed run programs. However, there are several evidences both theoretical and experimental which indicate that the SM could not be the ultimate theory. Instead, the SM should be viewed as a low energy effective theory of some more complete and fundamental one yet to be discovered. It is believed that a precise measurement of Higgs boson productions and decays can be a promising probe both to test the prediction of the SM as well as to search for new physics beyond the SM. After the discovery of the new scalar particle, there have been many theoretical and phenomenological studies devoted to non-minimal Higgs sector models that can explain such discovery and address some of the weakness of the SM. One of the simplest non-minimal Higgs model is the popular Two Higgs Doublet Model (2HDM) where both Higgs doublets possess a vacuum expectation value (VEV) and participate to electroweak symmetry breaking. A subclass of the 2HDM, that can provide a dark matter (DM) candidate, is the Inert Higgs Doublet Model (IHDM) where one Higgs doublet does not develop a VEV and may act as a dark matter candidate while the other one plays the role of the SM Higgs doublet [26]. We notice that the IHDM possesses an remarks and brief discussions in section 6.
2 Review of IHDM, its theoretical and experimental constraints 2

.1 A brief introduction to IHDM
The IHDM is a simple extension of the SM which can also provide a viable dark matter candidate. It is a version of the 2HDM with an exact discrete Z 2 symmetry. The SM scalar sector parametrized by H 1 is extended by an inert scalar doublet H 2 which can provide a stable dark matter candidate. Under Z 2 symmetry all the SM particles are even while H 2 is odd. We shall use the following parameterization of the two doublets : , with G 0 and G ± are the Nambu-Goldstone bosons absorbed, after electroweak symmetry breaking, by the longitudinal component of W ± and Z 0 , respectively. v is the vacuum expectation value (VEV) of the SM Higgs H 1 . The most general renormalizable, gauge invariant and CP invariant potential is given by : In the above potential, because of Z 2 symmetry, there is no mixing terms like µ 2 12 (H † 1 H 2 + h.c). In addition, by hermicity of the potential, all λ i , i = 1, · · · , 4 are real valued. The phase of λ 5 can be absorbed by a suitable redefinition of the fields H 1 and H 2 , therefore the scalar sector of the IHDM is CP conserving. After spontaneous symmetry breaking of SU (2) L ⊗ U (1) Y down to electromagnetic U (1) em , the spectrum of the above potential will have five scalar particles: two CP even H 0 and h 0 which will be identified as the SM Higgs boson with 125 GeV mass, one CP odd A 0 and a pair of charged scalars H ± . Their masses are given by: where λ L,S are defined as: From above relations, one can easily express λ i as a function of physical masses: 1 The IHDM involves 8 independent parameters: five λ 1,...,5 , two µ 1,2 and v. One parameter is eliminated by the minimization condition and the VEV is fixed by the Z boson mass, finestructure constant and Fermi constant G F . Finally, we are left with six independent parameters which we choose as follow : For completeness we list here the triple and quartic Higgs couplings that are needed in our analysis: As we will see later, these triple Higgs couplings are either directly involved in the processes e + e − → Zh 0 and e + e − → H 0 A 0 under investigation or in the experimental constraints that have to be fulfilled, while the quartic Higgs coupling H 0 H 0 A 0 A 0 enters only in e + e − → H 0 A 0 production.

Theoretical constraints
The parameter space of the scalar potential of the IHDM should be consistent with theoretical requirements. The important theoretical requirements in our consideration include perturbativity of the scalar quartic couplings, vacuum stability and tree-level perturbative unitarity conditions for various scattering amplitudes of all scalar bosons.
• Perturbativity: To guarantee the perturbation expansion, it is required that each of the quartic couplings of the scalar potential in Eq. (2) should obey the following conditions: • Vacuum Stability: The vacuum stability requires the potential V should remain positive when the values of scalar fields become extremely large [26]. From this condition, we have the following constraints on the IHDM parameters (for a review see [81]): λ 1,2 > 0 and λ 3 + λ 4 − |λ 5 | + 2 λ 1 λ 2 > 0 and λ 3 + 2 λ 1 λ 2 > 0 (9) • Charge-breaking minima: Likewise, a neutral, charge-conserving vacuum can be guaranteed by demanding that [82] which is a sufficient but not necessary condition for the vacuum to be neutral. Because, a neutral vacuum can also be achieved for positive λ 4 − |λ 5 | with suitable µ 2 1 and µ 2 2 , but in this case the dark matter (DM) particle can not be neutral. Condition (10) avoids this scenario. 2 • Inert Vacuum: In order to insure that the CP-conserving minimum described earlier is the global one, we need to impose the following conditions [82]: • Unitarity: To constrain the scalar potential parameters of the IHDM, the tree-level perturbative unitarity is imposed to the various scattering amplitudes of scalar bosons at high energy. From the technique developed in [83], we get the following set of eigenvalues: We impose perturbative unitarity constraint on all e i 's. e i ≤ 8π , ∀ i = 1, ..., 12.

Experimental constraints
The parameter space of the scalar potential of the IHDM should also satisfy experimental search constraints. We will consider the following experimental constraints: 1) from Higgs data at the LHC, 2) the direct collider searches from the LEP, 3) the indirect searches from electroweak precision tests, and 4) the data from dark matter searches. Below we elaborate more of these constraints.

Constraints from Higgs data at the LHC (HD)
In the IHDM, because of the exact Z 2 symmetry, H 2 does not couple to SM fermions which lead to natural flavor conservation. Only SM Higgs doublet couples to fermions, therefore all SM Higgs couplings to fermions and gauge bosons W and Z are the same as in the SM. Therefore, Higgs production cross section through conventional channels at the LHC such as gluon fusion, vector boson fusion, Higgsstrahlung and tth 0 are exactly the same as in the SM. Similarly, all tree level Higgs decays h 0 → bb, τ + τ − , ZZ * , W W * are identical to SM. Nevertheless, Higgs data provides important constraints to the IHDM.
. So condition (10) implies that m H 0 ≤ m H ± . In the other case when A 0 is the DM candidate i.e (HD.1) The one-loop decay channels h 0 → γγ and h 0 → γZ can receive contributions from charged Higgs boson via loop processes which may modify the SM predictions [78]. Therefore, in our analysis we will take into account the existing constraints on h 0 → γγ (see Sec. 2.4 for details).
(HD.2) Besides, when a DM candidate is lighter than m h 0 /2, the Higgs decay into invisible channel can be open. To be precise, in the IHDM, h 0 can decay into a pair of dark Higgs χ = H 0 , A 0 : h 0 → χχ if kinematically allowed where χ is the LOP. Such invisible decay of the SM Higgs has been investigated experimentally by both ATLAS and CMS and set an upper limit on the branching ratio: Br(h 0 → invisible) to 26% (resp 19% ) at 95% confidence level (CL) for ATLAS (resp for CMS) [84][85][86].
Moreover, global fit studies performed on LHC data can in turn put limits on the invisible decay of the SM Higgs which is constrained to be less than about 8.4% [87].
Note that the upper limit on the invisible decay of the SM Higgs can be inverted into a limit on the coupling h 0 χχ [68]. A more recent global fit [88] found that a tighter constraint on BR(h 0 → invisible) 5% at 95% CL is achievable. In this work, we adopt the most recent result from ATLAS [89] for which we have BR(h 0 → invisible) < 11% at 95% CL.

Direct search from LEP (DS)
The direct search bounds from collider signatures of the inert Higgs bosons at hadron collider or at lepton collider are rather similar to charginos and neutralinos production of the Minimal Supersymetric Standard Model (MSSM) [90,91]. We will follow the strategy adopted in [78,79,90,91]. These constraints can be roughly summarized as follows: (DS.1) m H ± > 80 GeV (adapted from charginos search at LEP-II), (DS.2) max(m A 0 , m H 0 ) > 100 GeV (adapted from neutralinos search at LEP-II),

Electroweak Precision Tests (EWPT)
The EWPT is very sensitive to extra electroweak multiplets, which can contribute to the vacuum polarization processes. Therefore it is reasonable to constrain the Higgs spectrum of the IHDM by using the global electroweak fit through the oblique parameters S, T and U [92]. The contribution to S and T parameters of an extra weak Higgs doublet [32] can be written as where the function F (x, y) is defined by To study the correlation between S and T, we define the χ 2 test as follows: Where S and T are the computed quantities in the IHDM given above, while S and T are the experimental central measured values of ∆S and ∆T , σ S,T are their one-sigma errors and ρ is their correlation. Using the PDG values of S and T with U fixed to be zero, we allow S and T parameters, to be as follows [93] : with the correlation coefficient ρ S,T = 0.92.
4. DM relic density, direct, indirect and collider searches (DM) (DM.1) From Plank data, the dark matter relic density Ωh 2 has been determined as 0.1200 ± 0.0012 [94]. In the current work, we set this value as an upper bound since it is possible that some particles like the right-handed neutrinos or the axion [95][96][97] can contribute to the relic density of our universe.
(DM.2) Generally speaking, dark matter can be detected by direct search (via the dark matter particle scattering with nuclei), indirect search (via a pair of dark matter particles annihilating into cosmic rays), and collider search (via mono-jet, mono-photon, mono-W ± /Z 0 /h 0 boson processes to measure the missing energy). In micrOMEGAs 5.2 [98], experimental constraints from XENON1T [99], PICO-60 [100], CRESST-III [101] and DarkSide-50 [102] have been implemented. From 6 GeV to 1000 GeV, XENON1T provided the most stringent bounds. We will use these bounds in this work. These bounds are obtained by assuming equal proton and neutron spin-independent cross section and assuming a specific choice of astrophysical parameters (i.e the DM velocity distribution is Maxwellian). In the IHDM, dark matter candidates can be either A 0 or H 0 , they can scatter with nuclei via t-channel by exchanging the SM-like Higgs boson (so-called Higgs portal). At leading order, the scattering cross section is isospinsymmetric and spin-independent. Therefore, these bounds will be taken into account.
(DM.3) At the LHC, one can search for dark matter candidate directly by looking for events with high p T monojet balanced by a large missing transverse energy [103,104]. The mono-jet final states can arise from one of the following processes: gg → χχ + g, qg → χχ + q and qq → χχ + g where χ = H 0 , A 0 is the dark matter candidate. In the IHDM, most of these processes can proceed by producing SM Higgs in association with monojet: {gg, qq} → h 0 g or h 0 q followed by the decay h 0 → χχ. Therefore these monojet processes are proportional to λ L/S . An attempt to set a limit on IHDM parameters using the ATLAS and CMS data was done in [58] where a projection for future LHC run including the high luminosity option is also provided. Ref.
[60] did a recasting for ATLAS analysis and set a limit on λ L/S . It is found that for m χ < m h /2, λ L/S > 3 × 10 −2 is excluded [58,60] while for m χ > m h /2 the limit is weaker and λ L/S > 5 apply only to the case where 62.5 < m χ < 100 GeV [60]. For higher dark matter mass, the cross section is suppressed, therefore the constraints on λ L/S would be weaker.

More about the constraints from h 0 → γγ
We now discuss the impact of the LHC experimental searches on the IHDM. Taking into account the latest measurement of the di-photon signal strength, we study the constraint on the charged Higgs mass and h 0 H + H − = −vλ 3 coupling that are involved in h 0 → γγ with h 0 being the SM Higgs. Since in the IHDM, the Higgs boson production cross section is identical to the SM one. Therefore, the di-photon signal strength reduces to the ratio of branching fractions of Br(h 0 → γγ) in the IHDM and in the SM: Moreover, in case where the decay h 0 → invisible is not open, the above ratio reduces to: where Γ(h 0 → γγ) is the partial decay width of h 0 decay into two photons. The measured signal strength relative to the SM expectation from ATLAS and CMS are given respectively by: [105] and µ CM S γγ = 1.18 +0.17 −0.14 [106]. The combined ATLAS and CMS result at the 2 σ level is given by: In the SM, it is well known that Γ(h 0 → γγ) is dominated by the W loops. In the IHDM, the charged Higgs loops can interfere constructively (respectively destructively) with the W loops for λ 3 < 0 (resp λ 3 > 0). In Fig. (2.1) (left), we present R γγ in the (m H ± , λ 3 ) plane. It is clear that for light charged Higgs boson, its contribution is rather important and could violate R γγ measurement. This is translated into a severe constraint on λ 3 . Namely, for m H ± = 200 GeV, λ 3 is forced to be in the range [−3, 2]. The allowed range for λ 3 becomes larger as far as the charged Higgs mass increase. When the charged Higgs is heavier than 400 GeV, the contribution of charged Higgs boson loops will be suppressed, i.e. in the decoupling limit, there is no limit at all on λ 3 .
It is remarkable that there exists a small parameter region where a light charged Higgs boson is allowed, i.e. the light charged Higgs boson can be around 100 GeV but with a large and positive λ 3 (say λ 3 > 10). For such a large λ 3 , the contribution of charged Higgs boson is around twice larger than that of W bosons with opposite signs. While in the SM case, the contribution of W bosons is the dominant one while the top contribution is subleading. In Fig. (2.1) (right) we show the allowed region in (m H ± , µ 2 2 ) plane. From this plot one can see that values of −40 × 10 3 GeV 2 ≤ µ 2 2 ≤ 0 are excluded for charged Higgs mass 100 GeV ≤ m H ± ≤ 420 GeV. While, values of 28 × 10 3 GeV 2 ≤ µ 2 2 ≤ 90 × 10 3 GeV 2 are not allowed in a mass range of 100 GeV ≤ m H ± ≤ 170 GeV.

Allowed parameter space and selected Scenarios
Before ending this section, we present the effect of the various aforementioned theoretical and experimental constraints on the IHDM parameter space. To be precise, the degenerate spectra are defined as the case where all Hidden Higgs bosons are degenerate, i.e.
According to Eq. (5), we have λ 4 = λ 5 = 0 while λ 3 could be either positive or negative depending on the splitting between m H ± and µ 2 2 . Then IHDM is fully described by three parameters which are: We perform a systematic scan over these three parameters taking into account all the above theoretical constraints. We notice that from vacuum stability constraints λ 2 must be positive. Using unitarity constraints e 9,10 one gets the strongest limit on λ 1,2 which gives λ 1,2 ≤ 4π 3 . For non-degenerate spectra, there will be no relation of Eq.  In Fig. (2.2), there are only three free parameters and we present the allowed parameter space in (λ 2 , m S ) plane with µ 2 2 colored in Fig. (2.2a). It should be pointed out that large and negative value of µ 2 2 is excluded both by unitarity constraints as well as by the inert vacuum constraints Eq. (11). Only small region with negative µ 2 2 survives. As one can see, the decoupling limit m S m Z is achieved for very large µ 2 2 , namely µ 2 2 > 10 6 GeV 2 . On the other hand, as discussed before, the size of h 0 H + H − triple coupling is directly related to the value of λ 3 . Therefore, it is interesting to know the allowed space for λ 3 and its correlation with other parameters. In Figs. (2.2b) and (2.2c), we show the (λ 3 , m S ) plane with µ 2 2 colored on the vertical axis and (m S , µ 2 2 ) plane where the values of λ 3 are color coded as indicated on the right of the plot, respectively. It can be seen from the figures that above decoupling limit could be reached for quite a wide range of λ 3 . It is also clear that λ 3 could be either positive or negative depending on the splitting between m S and µ 2 2 . However, vacuum stability constraints, request that λ 3 could not take large negative values. This is clearly seen in Fig. ( (2.2b) and (2.2c), but we can see that the upper bound of λ 3 is changed from 12 to 16, which indicates that larger h 0 H + H − coupling is allowed in this case. It is well known that the electroweak precision observables S and T put a strong constraint on the splitting between the masses that contribute to these parameters. This is clearly illustrated in Figs. (2.3d-2.3f). From Fig. (2.3d), one can see that the splitting between A 0 and H ± can not exceed 600 GeV. One can also see from this panel that the decoupling limit with large CP-odd and large charged Higgs can be reached for large µ 2 2 . It is also visible that in the decoupling limit, the splitting between A 0 and H ± becomes small and does not exceed 100-200 GeV. In Figs. (2.3e) and (2.3f), we illustrate the values of S and T parameters on the vertical axis as a function of m A 0 and m H ± . One can see that in the degenerate case, as expected, both S and T parameters vanish along the diagonal line m A 0 = m H ± . Away from the diagonal line, S and T get non vanishing values but remain within the allowed range.  In Fig. (2.4), the allowed parameter space in Scenario III is shown. It is observed that the constraints from dark matter conditions can kill 99% of points in Fig. (2.3). It is remarkable that from last two plots in Fig. (2.4), it can be seen that they are the same despite the permutation of m H 0 to m A 0 in th x-axis. As mentioned before, the independent parameters in this work are chosen as in Eq. (6). Actually it is found that not only all the allowed parameter space, but also all the results of the processes studied in this work, are totally symmetric under the exchange m H 0 ↔ m A 0 . In the following sections, we will consider 5 Scenarios, which are tabulated in  We will propose benchmark points from Scenario III and V and examine their radiative corrections, since both of them can pass all current experimental bounds.

Lowest order results
In our calculation, due to the tininess of electron mass and the corresponding Yukawa couplings, it is justified numerically to neglect the contributions of the Feynman diagrams which involve e + e − h 0 , e + e − G 0 , e − ν e G + and e + ν e G − vertices. For this reason, at the tree-level, the only one contributes in these two processes is the s-channel Z-exchange diagram, as shown in From the covariant derivative of the Higgs doublet, one can derive the Higgs coupling to gauge bosons. We list hereafter a part of the Lagrangian needed for our study where c W = cos θ W , s W = sin θ W . We stress first that all the above couplings are fixed just in terms of gauge coupling. For e + e − → Zh 0 , it is easy to compute the differential cross section which is given by [107]: where The total cross section is obtained after integration over the scattering angle. The analytical result can be found in [107]: The total cross section for the associate production e + e − → H 0 A 0 is given by: Because of the presence of two scalars in the final state, e + e − → H 0 A 0 process has the suppression factor κ 3 A 0 H 0 which reduces tremendously the cross section. In the Higgstrahlung process e + e − → Zh 0 , since there is no mixing between the two CP even neutral Higgs bosons, the cross-section is the same as in the SM. It is clear that σ 0 (Zh 0 ) scales like 1/s and is significant only at low energy just after the production threshold √ s ≈ m h 0 +m Z . We stress that the term κ 2 Zh 0 in σ 0 (Zh 0 ) originates from the longitudinal component of the Z, therefore one could conclude that at high energy the cross section is dominated by the longitudinally polarized Z boson. On the other hand, the production cross section for σ 0 (H 0 A 0 ) drops quickly due to the phase space suppression factor κ 3 H 0 A 0 (see Fig. (4.6)). Similarly, we draw in Fig. (3.3) the generic Feynman diagrams for e + e − → H 0 A 0 and the corresponding counter-terms. For this process, the self energies are not drawn, they are similar to the previous process e + e − → Zh 0 .

e
Evaluation of the one-loop corrections will lead to ultra-violet (UV) as well as IR divergences. The UV singularities are regularized with dimensional regularization and treated in the on-shell renormalization scheme while the IR ones are regularized with a small fictitious photon mass λ and cancelled with real photon emissions.
For both processes, owing to Lorentz invariance, the mixing Z µ -G 0 is proportional to (p µ 1 +p µ 2 ) where p 1 and p 2 are the momentum of the external electron and positron. After contracting Z µ -G 0 mixing with the initial state vertices e + e − Z and using Dirac equation the amplitude would be proportional to the electron mass which is neglected in the present calculation.
The presence of Z 2 symmetry forbids the mixing between the SM doublet H 1 and the inert doublet H 2 which tremendously eases the renormalization of the IHDM. The full renormalization of the IHDM has been presented recently in [34]. In our study, we will use the on-shell scheme developed first for the SM in [108][109][110] for all SM parameters supplemented by an on-shell renormalization for the extra-inert Higgs fields and their masses. Concerning the renormalization of the SM parameter and fields such as: the electric charge, the on-shell definition of the W and Z masses, γ-Z mixing and Weinberg angle, we refer to [110]. For the renormalization of the inert Higgses, we use similar approach as in [110]. Because of Z 2 symmetry, there is no mixing between h 0 -H 0 , G 0 -A 0 and Z-A 0 . This simplifies the renormalization of the Higgs fields. Let us redefine the new Higgs fields and masses as follows: Inserting these redefinitions into the above Lagrangian in Eq. (25), we find the following counter terms: For the counter terms of the initial state vertices e + e − γ and e + e − Z, counter terms of the Z boson, the photon propagators and their mixing, they are exactly the same as in the SM and can be found in [110]. We first stress that the counter-term for e + e − h 0 would be proportional to the electron mass and then vanishes for m e → 0. Second, in the case when both e − and e + are on shell, the one loop contributions coming from the e + e − h 0 vertex are vanishing in the limit of zero electron mass. However, from diagrams G15 and G16, one of the fermions (e − or e + ) coupled to the Higgs boson, is off shell. Consequently, for vanishing m e the correction to the vertex e + e − h 0 are UV-finite but non-vanishing (because, in this case there are contributions where the suppression factor m e /m W is absent). We have checked analytically this feature for this process. The remaining part, is of the same order as the other Feynman diagrams and should be included in the computation.
The Higgs wave function renormalization constants and mass counter-terms are fixed by the on-shell conditions for the masses and the Higgs fields and also by requiring the residue =1 for the Higgses. These requirements will lead to: However, all the counter terms for the gauge boson masses and wave function renormalization and their mixing as well as Weinberg angle are fixed following Ref. [110].
The electric charge renormalization constant δZ e is fixed from the e + e − γ vertex. We require that the renormalized three point functionΓ µ e + e − γ satisfies at the Thomson limit: Γ µ e + e − γ ( p 1 = p 2 = m e , q 2 = 0) = ieγ µ , and the renormalization constant for electric charge δZ e is obtained as [108][109][110] There is no reliable theoretical predictions available to extract Π AA hadron (0), but this quantity can be extracted from the experimental data. A non-perturbative parameter ∆α (5) hadron (m Z ) is used to absorb the hadronic contribution, namely δZ e is rewritten as Another popular scheme, α(m Z ) scheme, is more preferred, in which the large logarithm from leptons are also absorbed into the redefinition of running coupling constant [73,110,111]. The corresponding renormalization constant can be converted from α(0) scheme as: and the running coupling constant is replaced with Let us now discuss the treatment of the IR divergences. In fact, the IR divergences are present in two sources: i) wave function renormalization of charged particles such as electrons; ii) vertex corrections to e + e − γ and e + e − Z: Fig. (3.2)-G 17 and Fig. (3.3)-G 10 with V = γ, where incoming electron and positron exchange an virtual photon.
As mentioned before, in our calculation, to deal with the IR divergences, a small fictitious photon mass λ is introduced to regularize the soft and the virtual emission of the photon. Meanwhile, two cutoffs, ∆E and ∆θ, are introduced to deal with the IR singularities in real photon emission process. ∆E = δ s √ s/2 defines the soft photon energy cut-off for the bremsstrahlung process. It can be viewed as the photon energy cut that separates the soft from the hard radiation. The angle ∆θ cut defining between photon and the beam θ γ is used to separate hard radiation into hard collinear and hard noncollinear parts.
With λ, ∆E and ∆θ, the next-leading-order (NLO) corrections are decomposed into the virtual (V), soft (S), hard collinear (HC), and hard non-collinear (HC) parts as follows: Here V denotes the virtual correction including loop diagrams and counter terms from renormalization. CT denotes the "counter term" from electron structure function, originated from the 2nd term in Eq. (44).
Notice that the soft bremsstrahlung for e + e − → Zh 0 and e + e − → H 0 A 0 processes can be found in the literature. For completeness, we give the analytical expressions in the Appendix A, as well as many other details. The λ independence can be checked when combining the soft bremsstrahlung (S part) with the virtual one-loop QED contribution (V part) and this has been verified numerically with a good precision. Notice that we found a good agreement when we compare the result from FDC and FormCalc. Moreover, we have numerically checked that our results do not depend on ∆E, ∆θ and log(m e ), as shown in Fig. (A.1).
The total cross section at NLO, σ N LO , is the sum of LO cross section σ 0 , and NLO corrections σ 1 , namely where ∆ is the relative correction. Thus ∆ can be decomposed into two gauge-invariant parts, In the calculation of dσ V , computation of all the one-loop amplitudes and counter-terms is done with the help of FeynArts and FormCalc [112][113][114] packages. Numerical evaluations of the scalar integrals are done with LoopTools [115,116] and we have also tested the cancellation of UV divergences both analytically and numerically. The soft part dσ S is also done with FormCalc, while dσ HC+CT and dσ HC are obtained with the help of FDC [117].

Numerical results
In this section, we present our numerical results for the two processes introduced above. We adopt the following numerical values of the physical parameters from PDG [93]  In the IHDM, the CP even Higgs boson h 0 is identified as the Higgs like particle observed by the LHC collaborations and we use m h 0 = 125.18 GeV. For the other IHDM parameters, we perform a systematic scan which include the physical masses m H 0 , m A 0 and m H ± , λ 2 and µ 2 2 parameters. We take into account all theoretical requirements given in the subsection 2.2 as well as all experimental constraints given in the subsection 2.3. It is found that our numerical results are almost independent of λ 2 . Therefore, in the following part, we will fix λ 2 = 2.
In what follows, we will use the α(m Z ) scheme described before to present our numerical results.
Here we investigate them in the proposed scenarios given in Table (2.1).
In Fig. (4.1), total cross section and relative corrections in Scenario I are shown. Three typical collision energies of future electron-positron colliders, namely: √ s = 250 GeV, √ s = 500 GeV and √ s = 1000 GeV are chosen to present the results. Once the mass of charged Higgs boson is fixed, the triple Higgs boson couplings are simply determined by the parameter µ 2 2 , as given in Eq. (5). Accordingly, Scenario I also means that λ 3 can be non-vanishing while λ 4 = λ 5 = 0.
Therefore, in the upper panels of Fig. (4.1), the effect of the triple Higgs coupling λ 3 on the cross section is examined by varying the parameter µ 2 2 as shown in Table 4 Table 4.1. It is observed that weak corrections in the IHDM are typically negative and can reach 9% − 14% at √ s = 250 GeV and become 18% − 23% for √ s = 1000 GeV. When the mass of Higgs boson gradually increases to 500 GeV or so, the triple Higgs couplings proportional to λ 3 also increases, which leads to an increased new physics contribution to the total cross section. This seems to be against the decoupling limit, but it can be seen from Fig. (4.2) that the decoupling limit is reached at an even higher scale, around 1∼ 2 TeV. The enhancement is mainly due to the triple scalar couplings h 0 SS, S = A 0 , H 0 , H ± , which are proportional to λ 3 in Scenario I. Such triple couplings contribute into the corrections both linearly through the virtual corrections and also quadratically through the wave function renormalization of h 0 . A careful reader can also find that the starting points of the Higgs boson mass are also different for different values of µ 2 2 , which can be attributed to the theoretical and experimental constraints. IDM1 IDM2 IDM3 IDM4 IDM5 µ 2 2 (GeV 2 ) 40000 6000 0 -10000 -30000    Another comment is about the real emission, which can be clearly seen from the lower panels in Fig. (4.1). At the O(α(m Z )) order, it is found that the real emission contribution is independent of the new physics parameter and can be 0.5%, 1% and 1% for three collision energies, respectively. In order to illustrate the effect of radiative corrections in the IHDM and to avoid counting the pure SM effects, we define the following ratio given as: where σ IHDM Zh 0 and σ SM Zh 0 denote the one-loop total cross section in the IHDM and the SM, respectively. This ratio is useful for this process. We emphasize that there are only the contributions of new physics in the IHDM survived in numerator of the quantity δ while the full SM one loop effect has been subtracted. In terms of Feynman diagrams, the QED corrections as well as corrections to the initial state vertex e + e − Z Fig. (3.2)-G 17,18 , correction to e + e − h 0 vertex  In Scenario I, from upper panels of Fig. (4.2), one can read that the ratio of new physics with respect to the SM can change from 0.25% to −4.5% with √ s = 250 GeV, while they are similar for √ s = 500 GeV and √ s = 1000 GeV. Generally speaking, the contributions of new physics in the IHDM contributions reduce the production cross section, except from a tiny bump near the region with m S = 100 GeV.
In Scenario II, the contributions can change from −6% to 2% with √ s = 500 GeV and −5.5% to 3% with √ s = 1000 GeV. Depending upon model parameters, the contributions can either increase or decrease the production cross section.
The different behavior of Scenario I and II can be attributed to the different features of triple Higgs couplings. In Scenario I, according to Eq. (5), it is clear that both λ 4 and λ 5 should vanish. Then the triple Higgs couplings λ h 0 SS depends only on λ 3 which is severely constrained from di-photon signal strength limit (see discussion in section 2.4). In such a case, as shown in the first row of Fig. (4.2), there is clear dip near the region m S = 500 GeV where the ratio can reach from −1.5% to −4%. This occurs because the triple Higgs coupling λ h 0 SS = −λ 3 v which is driven solely by λ 3 becomes large in magnitude (say −8 ∼ −9). In fact, there are terms in δ which are proportional to the triple Higgs couplings, linear and quadratic. The linear terms can come from Feynman diagrams Fig.(3.2)-G 5 while the quadratic terms come from the wave function renormalization of h 0 which contributes to the counter-term of ZZh 0 vertex in Eq. (31). In contrast, in Scenario II, all λ 3,4,5 can contribute to the triple Higgs couplings λ h 0 SS , which can lead to more complicated interferences for each of the terms and can change the signs of new physics contributions.
It is interesting to notice that as shown in the upper panel of Fig. (4.2), there exist a large upside down peak near the region m S ∼ 500GeV and a tiny bump around m S = 150 GeV. When we look at these structures closely, we find that they are not caused by the mass of any particles nor by threshold effects. Instead, they are caused by the theoretical constraints such as unitarity, vacuum stability, no charged minima, etc. As shown in Fig. (4.3), for a given µ 2 2 (say the case µ 2 2 = 0 with √ s = 250 GeV), the allowed region of m S by all theses theoretical constraints produces a tiny bump near 150 GeV and ends at 500 GeV. Similarly with √ s = 500 GeV and 1 TeV cases, such structures appear. In other words, these accidental structures are produced by the cuts on m S .  It is necessary to point out that Fig. (4.2) also illustrates the decoupling behavior of IHDM in the process e + e − → Zh 0 . To demonstrate this, we include points with the Higgs masses m S up to 1 or 2 TeV and also take µ 2 2 in a wide range in order to satisfy theoretical constraints. As expected, the radiative corrections become smaller when m S increases, such a decoupling behavior can be clearly seen in Fig. (4.2) in the region with a large m S (say m S = 1 ∼ 2 TeV). In Scenario III, it is noteworthy that the parameter points with large radiative corrections are ruled out by the dark matter constraints, especially by the direct search bounds from XENON1T (Only 1% of points in Scenario II are still survived). Meanwhile, the radiative corrections can only lead to a change in cross section from −2.5% to 0.5% for √ s = 250 GeV case. Similar ranges hold for √ s = 500 GeV and √ s = 1 TeV cases. It should be emphasised that dark matter constraints indeed can significantly affect the allowed parameter space and the range of allowed radiative corrections of IHDM.
The results of Scenario IV and V are given in Figure (4.4). The band shapes in Figure  (4.4) are related to the fact that the mass region of dark matter particle is taken from 20 GeV to 62.5 GeV in Scenario IV and from 55 to 65 GeV in Scenario V (the mass range of 20-55 GeV is excluded due to too large relic density). It is observed that only 10% of points from the parameter space of Scenario IV can survive the dark matter searches constraints, which are displayed in Scenario V.  In Figure (4.5), we illustrate the angular distribution as a function of cos θ, where θ is angle between outgoing Higgs boson and electron beam in the center of mass energy frame. At high energy, it is well known that in the SM, the angular distribution behaves like sin 2 θ = 1 − cos 2 θ [118] and keeps the same shape at the NLO. This can be seen from Eq. (26) where the dominant term is κ 2 Zh sin 2 θ. In the figure, we illustrate both the LO and NLO distributions in the SM, as well as the angular distributions in the IHDM for six benchmark points, BP1−BP6, which are given in Table 5.1. From the plots, for √ s = 250 GeV, one can see that the results of BP1 and BP6 (BP2 and BP4) overlaps and for √ s = 500 GeV, the results of BP1 and BP4 (BP2 and BP6) overlaps. While for √ s = 1000 GeV, the results of BP1 and BP4 (BP2 and BP5) overlaps. In all cases, the IHDM contributions reduce the SM differential cross section. In the SM with 250 GeV, away from the forward and backward direction, the relative correction at NLO is about −10% to LO, while in the forward and backward direction, because of the box contributions, the relative correction could be slightly larger depending on the CM energy. As one can see, the angular distributions in the IHDM have the same shapes as those in the SM. For the √ s = 500 GeV and 1 TeV cases, the curves of IHDM almost overlap with the one of the NLO SM, i.e. the difference is very subtle and will be challenging to measure. While for the √ s = 250 GeV case, sensible deviations from the NLO SM can be observed which are detectable by experiments hopefully.
Before ending this section, we would like to stress that the case of 350 GeV CM energy is quite similar to that of 250 GeV case and is omitted here.
In the general 2HDM or in the MSSM, there exists a sum rule between the two vertices Zh 0 A 0 and ZH 0 A 0 , which itself simply reflects of the mixing between two CP even Higgs bosons h 0 and H 0 and the mixing between two CP odd Higgs boson G 0 and A 0 as well. The sum rule implies that the process e + e − → H 0 A 0 and the process e + e − → h 0 A 0 could always happen together for some specific choise of the mixing angles. Thus it is natural to expect that both processes could be detected at the future electron-positron colliders.
In contrast, in the IDHM, due to the fact that there is neither mixing between two CP even Higgs bosons h 0 and H 0 , nor mixing between two CP odd Higgs bosons G 0 and A 0 . A natural consequence from this fact is that the process e + e − → h 0 A 0 is forbidden. Therefore, to detect the signature of e + e − → H 0 A 0 and to prove that there is no e + e − → h 0 A 0 occurred at the same time can help to distinguish the IDHM from other general 2HDM like the MSSM.
It was pointed out in the Refs. [119,120] that triple Higgs couplings can greatly enhance the tree level cross section of e + e − → h 0 A 0 /H 0 A 0 in the general 2HDM. Such processes are supposed to help to probe the structure of Higgs potential of the 2HDM. Below, we examine the radiative correction to the cross section of e + e − → H 0 A 0 . In this process, there is no SM results, only the ratio ∆ defined in Eq. (41) is used.
We first give the numerical size for the tree level cross section. In Scenario II, Fig. (4.6) is to show the cross section of e + e − → H 0 A 0 in a scatter plot as a function of m H 0 and m A 0 . As shown in Eq. (29), the cross section depends only on the collision energy, the mass of H 0 and the mass of A 0 .
For Scenario I, we show the total cross section of e + e − → H 0 A 0 in the IHDM with three typical collision energies in Fig. (4.7). From the upper panels, it is observed that when the collision energy is fixed, the total cross sections decreases with the increase of the Higgs masses. The same values of µ 2 2 given in Table 4.1 are chosen to show the effect of triple Higgs couplings.   It can also be observed in the figure that the allowed ranges of m S are different for different values of µ 2 2 due to theoretical and experimental constraints. From the left plots of lower panel, it is found that for CM energy √ s = 250 GeV, the weak corrections can be −5.8%, −6% and −8% when µ 2 2 are chosen 0, 4 × 10 4 GeV 2 and 6 × 10 3 GeV 2 , respectively. For √ s = 500 and 1000 GeV cases, radiative corrections become larger and change dramatically near the threshold regions, where m S ∼ √ s/2. For instance, in the case when √ s = 1000 GeV and µ 2 2 = 0, the radiative corrections change from −12% to 30%. The change occurs in the mass range from 400 GeV to 500 GeV. Similar thing happens to the cases of other two µ 2 2 . A careful reader might notice that the behaviour of endpoints looks different for the three CM energies. As a matter of fact, the different behaviour near endpoints in Fig.  (4.7) are simply determined by the collision energy and the mass range of the new particles. Obviously, with a higher CM energy, the e + e − machine can cover a larger region of parameter space in IHDM.
Therefore, the CM energy √ s = 1000 GeV can cover the largest region of parameter space than those of √ s = 500 and √ s = 250 GeV cases. Although the CM energy changes the loop integral functions and causes its values to change, it is observed that when m S is near 120 GeV, ∆(e + e − → H 0 A 0 ) is similar for these three cases of CM energies.
Similarly when m S = 240 GeV, which is not reachable for √ s = 250 GeV case, the behaviour of ∆(e + e − → H 0 A 0 ) near the endpoint is similar for both √ s = 500 GeV and √ s = 1000 GeV. In the case √ s = 1000 GeV near the region m S = 400 GeV, it is noticed that there exists a fixed point where new particle's contribution is independent of other parameters of the IHDM. When m S becomes larger than 400 GeV, the triple Higgs couplings become larger and sizeable, and the contribution of the IHDM can even produce a positive value near the threshold. Such a behaviour can also be observed in Figure (4.8) for Scenario I with √ s = 500 GeV at the region m S = 180 GeV.
It is also worthy to mention that as shown in Fig. (4.7), the ratio of QED corrections is quite small and almost independent of the Higgs mass while the ratio of weak corrections depend significantly on the scalar masses and could be quite large.
For Scenario I, II and III, we show the ratio of the weak corrections in the whole allowed parameter space with corresponding triple Higgs couplings in the color bar in Fig. (4.8).
For Scenario I, in the case with √ s = 500 GeV, the corrections start from −9.5% or so when m S = 100 GeV. When m S increases from 100 GeV to 240 GeV, the ratio keeps increasing and can reach −1.5%. In the case √ s = 1000 GeV, the corrections start from −12% or so when m S = 100 GeV. When m S increases from 100 GeV to 490 GeV, the ratio keeps increasing and can reach 35%.
For Scenario III, less than 1% of points from Scenario II can survive when the constraints of dark matter are implemented, and the corrections can only be negative. It should be also emphasised that the range of the radiative corrections for those allowed points is shrunk by the dark matter constraints significantly. For example, for √ s = 1 TeV, before the dark matter constraints, the allowed range of radiative corrections can spread from −25% to 40%. But after the dark matter constraints, the allowed range can only change from −18% to −8%. All points with positive radiative corrections have been killed by the dark matter constraints. There are a couple of comments on Figure (4.8): • As shown in Scenario I, the ratio of weak corrections is within the range from −8.5% to −2.5% in the √ s = 250 GeV case, −10.5% to −2% in the √ s = 500 GeV case, and −15% to 35% in the √ s = 1000 GeV case. Whether such corrections can be detected at future electron-position colliders is determined by the production cross section. Meanwhile, the increase in the magnitude of the ratio when collision energy increases from 250 GeV to 1000 GeV does not mean the breakdown of the perturbation expansion. As a matter of fact, it more or less demonstrates the decrease of the LO cross section. Moreover, the larger collision energy also means a larger theoretical parameter space can be probed.
• As shown in Scenario II, in the √ s = 250 GeV case, the ratio is around −9% when triple Higgs coupling is around 1 ∼ 1.5; while it is around −1% when the coupling is around −2 ∼ −1.5. For √ s = 1000 GeV case, when m H 0 and m A 0 are larger than 400 GeV and when the triple Higgs couplings normalized to the vev become larger than −14, the ratio becomes larger than +40%.
• As shown in Scenario III, only 1% points from the parameter space can survive, which have been displayed. The radiative corrections can vary from −7% to −5.5% for √ s = 250 GeV, from −11% to −7% for √ s = 500 GeV, and from −18% to −8% for √ s = 1000 GeV. Points with positive corrections, like those shown in Scenario II for √ s = 1000 GeV case have been removed by the dark matter constraints.
For Scenario IV and V, we display the scatter plots of the allowed points from the parameter space in Figure (4.9). Similar to Figure (4.4), the band shapes are related to the mass range of the dark matter particles. In Scenario IV and V, the radiative corrections are always negative. It is noteworthy that there is a considerable amount of points that have been ruled out by the mono-jet constraints, which is different from Scenario III where mono-jet constraints has no perceivable affects to the parameter space.  In Table 5.1, we propose six benchmark points for the future e + e − collider search. Among them, BP1-BP3 belong to Scenario V and BP4-BP6 belong to Scenario III. BP1-BP3 can also be examined at the LHC or full higher energy pp colliders, via mono-jet measurement as shown in [58] or mono-W/γ signal as shown in [121]. BP1-BP3 are chosen such that the LOP is the CP-even Higgs boson H 0 . It has been checked that we can exchange the mass of H 0 with A 0 . For BP4-BP6, the invisible decay of the SM Higgs is not open.

Benchmark points
A few more explanations on these BPs are provided below: • In these BPs, BP1 represents a favourable case where the loop correction to the process e + e − → Zh 0 can be detectable and the new particles, rather light, can be directly produced at the future Higgs factories via the process e + e − → H 0 A 0 at all energy cases.
• BP4 provides a case where for all the CM energies, the effects of new physics are small and difficult to be detected in the process e + e − → Zh 0 . Although there are light particles like H 0 , A 0 and H ± for BP4, the contributions of new particles yield a small correction to the cross section of e + e − → Zh 0 due to the smallness triple Higgs couplings. But when we consider the process e + e − → H 0 A 0 , it is possible to produce these new particles directly at the future Higgs factories.
• BP6 provides a case where the mass of new particles is too heavy to be produced directly at future e + e − colliders with CM energy less than 1000 GeV. But the effects of new physics can be detected via the loop effects in e + e − → Zh 0 . The reason for such a sizeable correction can be attributed to the large triple Higgs couplings for this case.
• BP3 and BP5 represent more complicated cases where new physics contribution induced via loop to the process e + e − → Zh 0 can be sizeable, but to confirm the case we need future colliders with high CM energies like 500 GeV or higher.
• BP2 represents the case that a 250 GeV Higgs factory might be difficult to detect the effects of NP, but Higgs factories with a collision energy larger than 350 GeV might be better. The production of new particles needs a CM energy higher than 500 GeV.
Except the collider searches, these BPs can also be searched by the future dark matter searches, which can be left for future study.
In Table 5.2, the total cross section for e + e − → Zh 0 in the SM with various center of mass energies are presented, as well as the total cross section for the benchmark points in the IHDM. For the SM results, we present the LO cross section, one-loop QED corrections, one-loop weak corrections and full NLO cross section. For the IHDM, as the LO cross section and one-loop QED corrections are exactly same as the ones in the SM, only one-loop weak corrections and full NLO cross section are presented, with the addition of ∆ and δ, where ∆ is the relative one-loop corrections to LO and δ is the relative correction of IHDM to the full NLO SM result defined in Eqs. (40) and (42), respectively.
There are a few comments on the results in Table 5.2: • As shown in the column δ, when the CM energy √ s = 250 GeV, the ratio of the contribution of new physics in BP1 and BP6 can increase the cross section by a factor 0.414% and 0.525%, respectively. In contrast, the contribution of the new physics in BP3 and BP5 decreases the cross section sizeably by a factor −1.221% and −2.294%, respectively. The contribution of the NP in BP2 and BP4 is −0.093% and −0.098%, respectively, which is less than the projected precision 0.1% of the future Higgs factories.
• When the CM energy increases to 350 GeV, the ratio in BP1 and BP6 keeps to be positive and increase the cross section slightly. The corrections in BP3 and BP5 decreases the cross section by a factor −0.994% and −2.140%. While the contribution in BP2 increases the cross section and the effect can reach 0.1%. However, the contribution in BP4 is still small and the effect is −0.091%.
• When the CM energy increases to 1000 GeV, the correction in BP1 becomes negative and is given by −0.071%, while in BP6 it remains positive and is 0.85%. The corrections of BP3 and BP5 are negative: −1.824% and −0.715%. While the one of BP2 and BP4 becomes −0.456% and 0.098%, respectively.
In Table 5.3, we present the LO and NLO results for e + e − → H 0 A 0 with various center of mass energies. We give the weak contributions and the QED ones as well as the total NLO cross section. We also show the relative corrections with respect to the LO results, as demonstrated by the results given in ∆ column. Roughly speaking, the contributions of NP always reduce the cross sections from −5.5% to −17% for these BPs. When the precisions of future Higgs factories are considered, obviously such large corrections must be taken into account in any experimental analysis. For example, in term of the cross section of BP4 and the projected precisions given in Table 1.1, the CEPC with √ s = 250 GeV can find 3.075 × 10 5 raw events of e + e − → H 0 A 0 , which corresponds to a precision of 0.18% on the cross section when only statistic errors are taken into account. Instead, the LO result predicts that a total number of events is 3.284 × 10 5 , which will lead to a deviation of 60 σ. Meanwhile, such a huge number of signal events can also lead to a precise measurement on the masses of H 0 and A 0 , and impose a strong constraint on the parameters in the loop.

Conclusions and Discussions
We have studied the one loop radiative corrections to the neutral processes: the Higgsstrahlung e + e − → Zh 0 and the neutral Higgs-pair production e + e − → H 0 A 0 in the IHDM. We have evaluated both the QED corrections, the soft and hard photon emissions and the full weak corrections. The Feynman diagrams are evaluated using dimensional regularization in the Feynman gauge. The full one loop analysis is done using the on-shell renormalization scheme. In addition, in the numerical analysis, we first performed a systematic scan over the IHDM parameter space taking into account theoretical as well as experimental constraints and localized allowed parameter space.
For e + e − → Zh 0 process, we first evaluated the one-loop radiative corrections in the SM and checked that they do agree with the existing results in the literature. Next we have evaluated the one-loop corrections in the IHDM and the relative corrections with respect to the one-loop SM result. We have shown that the pure IHDM effect could reach about −4.5% percent in Scenario I and could be slightly larger and reach −6% in Scenario II, and only at most −2% in Scenario III. It is remarkable that for Scenario III (V) after imposing the dark matter constraints, the allowed points in the parameter space significantly are reduced when compared with Scenario II (IV). Meanwhile, the range of the radiative corrections for the allowed points is also greatly shrunk by the dark matter constraints, as shown in the lower panel of Results are shown for 250 GeV, 500 GeV and 1 TeV center of mass energy. For 350 GeV, the situation is similar to 250 GeV. Such effect is large enough to be measured in the precise future linear collider program. We have also presented one-loop angular distributions for six BPs. In addition, we have demonstrated that for the heavy internal IHDM spectrum, the oneloop corrections decouple for large µ 2 2 . It has been demonstrated that the triple Higgs couplings h 0 SS with S = H 0 , A 0 , H ± which contribute into the one-loop virtual corrections as well as to the wave function renormalization of h 0 could contribute significantly to the relative corrections.
In the case of the pairwise production of neutral Higgs bosons e + e − → H 0 A 0 , since two scalars are involved in the final state, large effect is found mainly coming from triple Higgs couplings: either from the wave function renormalization of H 0 and A 0 or from the virtual correction Fig. (3.3)-G 1 . In the case of 250 GeV center of mass energy, the mass of H 0 and A 0 are restricted by m H 0 + m A 0 < 250 GeV, the λ i involved in the triple couplings could not be significant, the effect is rather small and could not exceed −9% both in Scenario I and II. While in the case of 500 GeV CM energy, the effect is slightly larger and could reach −15%. In the case of 1 TeV center of mass energy, with this energy one can cover a large range for m H 0 and m A 0 . Therefore, the quartic couplings λ i become large than in the previous cases which would make the triple scalar coupling large and the relative corrections to the tree level result could be of the order −30 → +40%.
In Scenario III, after taking into account the dark matter constraints, the range of weak corrections ∆ weak of allowed points is confined in a narrower range when compared with those of Scenario I and II. For example, for the case √ s = 250 GeV, ∆ weak of allowed points can only be in the range −6.8% to −5.5%. It should be emphasised that in both processes e + e − → Zh 0 /H 0 A 0 , as demonstrated by the six BPs, the radiative corrections in the allowed parameter space can be rather large. When the projected precisions of future e + e − colliders given in Table 1.1 are considered, it is mandatory to take them into account in any realistic experimental measurements and data analysis.
Here it would be interesting to explore whether the future e + e − colliders have the potential to distinguish different new physics models, like the IHDM, the general 2HDM, and the MSSM. For obvious reasons (like huge parameter space in MSSM), an exhaustive and thorough comparison is an impossible mission. Even a fair comparison is difficult. Instead, we confine to compare the results of a few scenarios considered in literature at their face values, and compile them in Table 6.1.
• In the scenarios of the general 2HDM, the one loop radiative corrections to e + e − → Zh 0 and e + e − → H 0 A 0 have been computed in Refs. [111,119,120]. We compile the results of scenarios in the recent Reference [111] for the purpose of comparison.
• In the context of supersymmetric models, the one loop corrections to e + e − → Zh 0 /H 0 A 0 have been presented within the MSSM and complex MSSM (CMSSM) [122][123][124]. We compare our results with those of scenarios presented in the Reference [122] the same ratio δ is defined and used there.  According to the projected precision presented in Table 1.1, the precision in measuring δg ZZh 0 /g ZZh 0 can reach 0.2% for future FCC-ee with √ s = 250 GeV, which roughly corresponds to a precision in measuring δ = ±0.4%. A recent global analysis given in [125] demonstrated an even more aggressive precision ±0.09% might be achievable at CLIC with a combo runs with three CM energy √ s = 380/1500/3000 GeV and polarized beams, which, roughly speaking, means an error in δ can reach ±0.18%. In order to address the issue of model discrimination, below we deliberately take an moderate optimistic assumption that a precision δ = ±0.2% can be achievable.
Obviously, given a possible precision, whether two models can be distinguished is determined by the central value of δ EXP . Suppose that the future experiment could determine δ EXP = 0.0%±0.2%, it could be able to rule out the parameter region of the MSSM analysed in [122] more than 5σ, while the parameter space of the CMSSM analysed in [122], the SM, 2HDM and IHDM might still be consistent with experimental bounds. In contrast, suppose the future experiment could determine δ EXP = −2.5%±0.2%, then the parameter region of the MSSM analysed in [122] and 2HDM and IHDM could interpret the experimental data while the parameter region of the CMSSM analysed in [122] and the SM could be ruled out in terms of 8σ and 12σ, respectively. While, suppose the future experiments could determine δ EXP = −5.5%±0.2%, then only 2HDM can interpret the experimental data while the other models could be ruled out with more than 10σ. In the case with δ EXP = 0.55% ± 0.2%, only IHDM could interpret with the experimental data comfortably and the deviation from the prediction of the SM can reach to 2.5σ.
For future searches at e + e − colliders, we have presented six benchmark points that satisfy dark matter and LHC constraints. For these BPs, we have given the weak and the QED corrections for various center of mass energies. These BPs can be explored at the future e + e − colliders, the LHC and future proton-proton colliders, and future dark matter experiments.
For example, the discovery channel of e + e − → H 0 A 0 can lead to some interesting signatures, as explored in the reference [42,47]. Moreover, for BP1-3, H 0 is the LOP and the final state of e + e − → H 0 A 0 would lead to ZH 0 H 0 final state since A 0 → ZH 0 dominantly, we expect a signature with dilepton plus missing energy. Since the Z boson is on-shell, then we expect to observe two energetic leptons and a large missing energy in the final state. For BP4, the Z boson is off-shell and two leptons from off-shell Z boson decay in the final state are soft. Then two soft leptons and large missing energy would be the characteristic signature of BP4. In the parameter space, there are some points where A 0 can have a tiny decay width and its lifetime can be larger than 10 −15 s, we expect that the signature could be displaced vertex and large missing energy, as shown in [126].
At the LHC, it is believed that radiative corrections to pp → W h 0 , Zh 0 , H 0 A 0 would be, to some extent, similar to our finding for e + e − colliders. Therefore, our BPs can also be explored at the pp colliders via the processes pp → H 0 H ± → H 0 H 0 W ±( * ) for BP1-6. The signatures of these BPs can be a large missing energy plus a W boson, while the W boson can be either on-shell (for BP1-3 and BP5) and off-shell (for BP4 and BP6). For both BP4 and BP6, the charged Higgs boson is almost degenerate with the LOP H 0 and has a large lifetime, which might lead to a signature of charged displaced vertex at the LHC. financial support where part of this work has been done. A More details about IR divergences As mentioned in section 3.2, IR divergences in this paper are regularized with a small fictitious photon mass λ. Meanwhile, two cutoffs, ∆E and ∆θ are used to separate the phase space of real photon emission. Thus full NLO corrections are separated into four parts, as given in Eq. (39). The hard non-collinear part, dσ HC , is obtained using traditional Monte-Carlo integration techniques. Here we present results for the soft and hard non-collinear parts, as well as the check on the independence of the cutoffs.
For the Higgsstrahlung e + e − → Zh 0 and the associate production e + e − → H 0 A 0 , the real photon is only emitted from the initial state electron and positron. The analytical expression for the soft bremsstrahlung is given by: where ∆E is the cut on the photon energy and λ is a small fictitious mass for the photon. In the above formula, the IR term log 4∆E 2 λ 2 (resp log is also canceled by the virtual QED diagram. One-loop radiation correction includes collinear singularities when m e goes to zero. In our calculation electron has nonzero mass, but the singularities will become terms proportional to log(m e ). Some of them are cancelled when summing up virtual and real corrections, and some of them are absorbed into the redefinition of running coupling constant as mentioned above, but some are remained. To deal with this, we used following fixed order electron structure function which can be derived [127] f ee (x, s) = δ(1 − x) + α 2π log s 4m 2 e P + ee (x, 0) with P + ee (z, 0) = being the regularized Altarelli-Parisi splitting function. The 2nd term in Eq. (44) gives an additional "counter term" which can be combined with hard collinear part. The HC + CT part is obtained as dσ HC+CT ≡ dσ * HC+CT + dσ SC dσ * HC+CT = α 2π 1 + z 2 1 − z log ∆θ 2 − 2z 1 − z × dσ 0 (zk 1 ) + dσ 0 (zk 2 ) dz dσ SC = − α π log s 4m 2 e 3 2 + 2 log δ s dσ 0 , where the approximation ∆θ m e / √ s has been taken and ∆E is replaced with a dimensionless parameter δ s = 2∆E/ √ s. ∆E and ∆θ are unphysical cutoffs we introduced to deal with IR singularities. Our final results should not depend on them. In Fig. (A.1), we show our check for this independence, taking one of the processes as an example. From first subfigure, it can be seen that the independence on ∆E is found in a wide range and we choose δ s = 10 −3 as our default choice. In second subfigure, we can see that the result becomes cut dependent when ∆θ is smaller than 10 −4 . It is because the approximation used in Eq. (46) demands ∆θ m e / √ s ∼ 2 × 10 −6 . Thus we choose ∆θ = 10 −3 as our choice.
Collinear divergences in our calculation appear as terms proportional to log(m e ). After including the counter term from electron structure function, such divergent terms should vanish in the final result. In order to check this, we vary the mass of electron with a factor of k from 2 −4 to 2 8 , namely m e is taken k × 0.511 MeV. The cancellation is shown in last subfigure of Fig. (A.1), from which we can see that the result remains unchanged when k varies . Also, we can see that singular terms only appear in σ V +S and σ SC parts.