Full one-loop radiative corrections to $e^+ e^-\to H^+H^-$ in the inert doublet model

We compute the full one-loop radiative corrections for charged scalar pair production $e^{+}e^{-}\to H^{+}H^{-}$ in the inert doublet model. The on-shell renormalization scheme has been used. We take into account both the weak contributions as well as the soft and hard QED corrections. We compute both the real emission and the one-loop virtual corrections using the Feynman diagrammatic method. The resummed cross section is introduced to cure the Coulomb singularity which occurs in the QED corrections. We have analyzed the parameter space of the inert doublet model in three scenarios after taking into account theoretical constraints, the collider experimental bounds, and dark matter search bounds as well. It is found that the weak interaction dominates the radiative corrections, and its size is determined by the triple Higgs coupling $\lambda_{h^0 H^+ H^-}$, which is further connected to the mass of the charged scalar. In the scenario where all the constraints are taken into account, we find that for $\sqrt{s}=250$ GeV and $\sqrt{s}=500$ GeV, the weak corrections are around $-6\% \sim-5\%$ and $-10\% \sim -3\%$, respectively. While for $\sqrt{s}=1000$ GeV, the weak corrections can reach $-15\% \sim +25\%$. The new feature is that the weak corrections can be positive near the threshold when the charged scalar is heavier than 470 GeV. Six benchmark points for future collider searches have been proposed.


Introduction
The Standard Model (SM) particle spectrum has been completed with the discovery of the Higgs boson on July 4, 2012, by the ATLAS and CMS experiments at CERN [1,2].Furthermore, this discovery has confirmed that the SM of particle physics is the underlying theoretical framework valid at least for energies up to the electroweak (EW) scale.The two collaborations also carried out several Higgs boson couplings measurements at the Large Hadron Collider (LHC) during run-1 and run-2 such as the couplings of the Higgs boson to top quarks [3,4], tau leptons [5,6], bottom quarks [7,8], and all the electroweak gauge bosons, including the decays to ZZ * [9,10], W W * [11][12][13][14], and γγ [15].In addition, recently upper limits have been set on the h 0 → γZ signal strength [16,17] and on the Higgs boson production cross section times branching fraction to muons [18,19].
The aforementioned measurements will be improved at future experiments such as the High-Luminosity LHC [20,21], scheduled to operate from 2029, where the Higgs boson couplings are projected to be improved to a precision level of 5 − 10%.In addition, the experimental uncertainties will be further reduced in the clean environment of the future lepton colliders, such as the International Linear Collider (ILC) [22,23], the Circular Electron Positron Collider (CEPC) [24], the Compact Linear Collider (CLIC) [25][26][27], and the Future Circular Collider [28].For example, at the ILC, with a c.m.energy of about 250 GeV and a luminosity of 2 ab −1 , some Higgs boson couplings will most likely be measured at a precision level of 1% for h 0 → b b and below 1% for h 0 → ZZ, W W [29].
Although new physics beyond the SM has not yet been established by the current LHC dataset, it is necessary in order to understand several puzzles of the SM and the observed Universe, such as what is the nature of dark matter, what is the origin of neutrino masses, how to stabilize the vacuum of the SM, and so on.From the viewpoint of effective theory, the success of the SM can only be understood as the low energy limit of a more fundamental theory.In this more fundamental theory, there must be extra sectors not included by the SM.For example, the Higgs sector of the SM is assumed to be minimal and composed of only one Higgs doublet.In order to solve the problem of CP violation, it is well motivated to extend the Higgs sector of the SM by introducing one more doublet [30] and the Higgs potential sector can be extended, such a theory model is called the two Higgs doublet model (2HDM).A comprehensive review of the 2HDM can be found in the reference [31].
In order to offer a dark matter candidate, a Z 2 symmetry is introduced into the 2HDM, which leads to the so-called inert doublet model (IDM).In the IDM, the second doublet does not develop a vacuum expectation value (VEV) nor does it have a direct coupling to the SM fermions.In such a Z 2 symmetry, the fermions, gauge bosons, and the SM Higgs doublet are invariant, and the second doublet is odd, i.e.H 2 → −H 2 .The Z 2 symmetry ensures that this extra doublet has no direct coupling to the SM fermions at both tree and loop levels, and its lightest stable neutral component may play the role of a dark matter candidate.The IDM was proposed by E. Ma et al. [32] and it was initially suggested for studies on electroweak symmetry breaking.This model is very appealing because it can generate small neutrino masses [33], can provide a dark matter candidate [34][35][36][37][38][39][40] and can solve the naturalness problem [41].The phenomenology of the IDM has been extensively studied in the literature in the context of the LHC and at future Higgs factories such as the ILC or CLIC  .
After the electroweak symmetry breaking, the IDM possesses five physical scalars: One Higgs boson h 0 which is identified as the 125 GeV SM Higgs, two new neutral physical scalars H 0 , A 0 and two charged physical scalars H ± .Both H 0 and A 0 could be dark matter candidates.
Obviously, the discovery of a charged scalar boson would be a clear sign of nonminimal Higgs sectors and precise knowledge of its production properties would be useful to reconstruct the full Higgs potential.Since all new IDM scalars are odd under the Z 2 symmetry, they must be produced in pairs at the colliders.Moreover, these new inert scalars only couple to the electroweak gauge bosons and the Higgs boson of the SM.Consequently, at the LHC they should be pair-produced via processes of the Drell-Yan type: q q′ → W * → A 0 H ± , q q′ → W * → H 0 H ± , and q q → Z * (γ) → H + H − , or from gg, q q → h 0 * → H + H − [53,72], where q′ represents a different quark flavor.Furthermore, the charged scalars in the IDM could also be produced in the same-sign pair production process pp → H ± H ± jj [52,81,82] or from vector boson fusion-like production pp → H ± H ∓ jj [53].At lepton colliders, charged scalars can be produced via the pair production process e + e − → H + H − [43,49,50].Future muon colliders [83,84] may have better potential to cover a wide range of the mass of charged scalars.
The discovery of a charged scalar would be clear evidence of beyond the SM physics.In the case such a discovery happens in e + e − collider, a subsequent precision measurement of its properties will be crucial to determine its nature.In order to obtain sufficient accuracy, oneloop corrections to the various charged scalar production modes and its couplings have to be considered.Recently, one loop radiative corrections to e + e − → Zh 0 [85][86][87] and e + e − → H 0 A 0 [85] in the IDM have been evaluated.In addition, there are many other works dealing with radiative corrections in the IDM either at one-loop order [85][86][87][88][89][90][91][92][93][94][95][96] or beyond the one-loop level [97][98][99].Full one-loop calculations to e + e − → H + H − in nonsupersymmetric models such as the 2HDM and some supersymmetric models like the real minimal supersymmetric standard model (MSSM) and complex MSSM have been presented in Refs [100][101][102][103][104][105][106][107].For example, it has been found in the complex MSSM [107] that radiative corrections to e + e − → H + H − can go up to 20% or higher, hence a full one-loop contributions are important for future linear colliders such as the ILC or CLIC.
In this work, we present for the first time a full one-loop analysis to e + e − → H + H − within the IDM.Besides the full weak corrections, we include the soft and hard QED radiation as well as the treatment of collinear divergences and the Coulomb singularity.In our scan and numerical analysis, we will take into account all the existing constraints on the IDM including the theoretical ones, the experimental constraints from the LHC, such as the Higgs data from the decay of the Higgs boson into two photons, the invisible Higgs decay, the electroweak precision tests, as well as the constraints from dark matter relic density in the Universe and the direct bounds from the monojet searches at the LHC.It should be mentioned that we have included the recast results from LEP II data [108] in our scan and analysis.Additional LHC recast results, such as dilepton and multilepton final states, as well as vector boson fusion final states, can be found in Refs.[62,63,109,110].These recasts might bring new limits to the IDM parameter space.However, in this work, they are not taken into consideration.
We perform a scan over the whole parameter space of the IDM in three scenarios: Scenario I assumes that the inert scalars are degenerate in mass and has taken into account all the constraints except the invisible Higgs decay and dark matter constraints, scenario II is defined by imposing all the constraints without dark matter ones and with a nondegenerate spectrum for the new inert scalars and scenario III is the same as scenario II but with dark matter constraints.We find that the dark matter constraints can greatly reduce the number of points in the parameter space.Furthermore, it is found that QED corrections are rather small and radiative corrections are dominated by the weak interactions.The size of weak corrections depends on the coupling λ h 0 H + H − , which is further related to the mass of charged scalar.In scenario III we find that for √ s = 250 and √ s = 500 GeV, the weak corrections are around −6% ∼ −5% and −10% ∼ −3%, respectively.While for √ s = 1000 GeV, the weak corrections can reach −15% ∼ +25%.The new feature is that the weak corrections can be positive near the threshold when the charged scalar mass is larger than 470 GeV.We have proposed six benchmark points for future lepton collider searches.The outline of this paper is as follows: In Sec. 2, we briefly describe the model, including its mass spectra, key trilinear and quartic scalar couplings, and list various theoretical and experimental constraints that we will take into account in this work.In Sec. 3, we provide the leading-order (LO) formula for the cross sections of the e + e − → H + H − process, introduce the on-shell renormalization scheme for the IDM, and set up basic notations and conventions.Moving on, we study the one-loop contributions to the e + e − → H + H − process and examine the importance of soft and hard photon emission in order to guarantee the cancellation of the infrared (IR) divergences at the next-to-leading order (NLO) calculation.Furthermore, we tackle the challenge posed by the Coulomb singularity using efficient resummation techniques.We present our numerical results in Sec. 4. In Sec. 5, we propose some benchmark points (BPs) and examine their radiative corrections for future e + e − colliders.We end this work with discussions in Sec. 6.

Review of the inert doublet model 2.1 A brief introduction to IDM
The IDM is one of the simplest extensions beyond the SM.This model has an extra doublet H 2 which is added to the scalar sector of the SM.This doublet does not generate any VEV and it does not have direct coupling to the fermions of the SM.An unbroken Z 2 symmetry is imposed such that fermions, gauge bosons, and the SM doublet are invariant while the additional scalar doublet is odd i.e.H 2 → −H 2 under this symmetry.The parametrization of the two doublets is given by where G 0 and G ± are the Goldstone bosons gauged out, after electroweak symmetry breaking, by the longitudinal components of W ± and Z, respectively.The v denotes the VEV of the SM Higgs doublet H 1 .
Then the renormalizable scalar potential can be given as: Note that because of the exact Z 2 symmetry, the above potential has no mixing terms like In addition, since the potential must be Hermitian, all λ i , i = 1, • • • , 4 are dimensionless and real whilst the phase of λ 5 can be absorbed by a suitable redefinition of the fields H 1 and H 2 .After spontaneous symmetry breaking of the group SU (2) L ⊗ U (1) Y down to U (1) em , we have five physical scalars: h 0 which is the SM 125 GeV Higgs boson and four inert scalars: H, A, H ± .Their masses are given by: where λ L,S are defined as From above relations, one can easily find the expressions of λ i as a function of physical masses: The IDM involves eight independent parameters: five λ 1,...,5 , µ 1 , µ 2 and the VEV.One parameter can be eliminated by using the minimization condition, while the VEV is fixed by the Z boson mass, fine-structure constant and Fermi constant G F .Finally, we are left with six independent parameters, which we choose as follows One alternative parametrization is to use λ L or λ S in place of µ 2 2 , as can be seen from Eq. ( 3).The advantage of such parametrization is that it is directly connected to the evaluation of the relic density.
For completeness, we list below the triple and quartic scalar couplings that are needed in our numerical analysis: It is clear that the triple coupling 2 ) as given in Eq. ( 5).

The IDM parameter space and constraints
In this work, we study in the same parameter space as in our previous work [85].The parameter space is obtained by scanning the whole space with certain theoretical and experimental constraints.The constraints used are summarized as follows:

Theoretical constraints
The parameters of the IDM are subject to several theoretical constraints that we shall impose throughout our numerical analysis.

• Perturbativity:
We require that each of the quartic couplings of the scalar potential in Eq. ( 2) is perturbative: • Vacuum stability: In order to ensure vacuum stability, the potential V should remain positive when the values of scalar fields become extremely large [32].From this condition, we have the following set of constraints on the IDM parameters (for a review, see [31]): • Charge-breaking minima: A neutral charge-conserving vacuum can be guaranteed by demanding that [111] • Inert vacuum: We impose the following conditions in order to insure that the inert vacuum is the global one [111]: • Unitarity: As in the SM, the couplings of the IDM have to satisfy certain relations in order to obey unitarity.The tree-level perturbative unitarity is imposed on the various scattering amplitudes of scalar bosons at high energy.From the technique developed in [112], we find the following set of eigenvalues: We impose perturbative unitarity constraint on all e i s: e i ≤ 8π , ∀ i = 1, ..., 12.
It is important to note that we did not take into account higher-order corrections to the bounded from below, global minimum, and perturbative unitarity constraints.We have only included the constraints at lowest order.However, the parameter space will change if we use those constraints at NLO.For example, it has been shown in [113] that the bounded from below conditions at one-loop level can change the stability of the electroweak vacuum, and the allowed parameter space for coexisting inert and inertlike minima at NLO is larger compared to the one observed at LO. Furthermore, for the unitarity constraints beyond the lowest order that have been analyzed in the 2HDM [114], it was found that perturbative NLO unitarity constraints are stronger than the LO ones.Therefore, the allowed ranges of the quartic λ i couplings, the inert scalar masses, and their mass differences will change if we use the NLO unitarity conditions.

Experimental constraints
The parameter space of the scalar potential of the IDM should also satisfy experimental search constraints.We will consider the following experimental constraints (for further details about these constraints see our previous published work [85]): • Higgs data at the LHC [115,116].
• The bound on the invisible Higgs decay [117].
3 Radiative corrections to e + e − → H + H − 3.1 e + e − → H + H − at tree level Due to the extremely small value of the electron mass and the corresponding Yukawa couplings, it is numerically justified to neglect the contributions proportional to the mass of electron, such as Feynman diagrams involving e + e − h 0 , e + e − G 0 , e − ν e G + , and e + ν e G − vertices.For this reason, the Feynman diagrams contributing to e + e − → H + H − at the tree level are mediated by the γ and Z s-channel exchange as shown in Figure 3.1.Using the covariant derivative of the Higgs doublet, one can derive the scalar coupling to gauge bosons.We list below a part of the Lagrangian needed for our study where e is the electric charge and c W ≡ cos θ W (s W ≡ sin θ W ) with θ W being the Weinberg mixing angle.It is clear from above equation that the interactions involving charged scalar depend only on the electric charge and the Weinberg mixing angle θ W .Let p 1,2 and k 1,2 be the momenta of incoming e ± and outgoing H ± , respectively, and define the Mandelstam variables as The matrix elements at tree level for the two contributions have the following form: with The total matrix element is then given by and the total cross section can be written as where β = 1 − 4m 2 H ± /s is the velocity of outgoing H ± in the c.m. frame.Note that the first term comes from the photon exchange, the second term comes from the Z boson exchange while the third one is the interference between the photon and Z boson exchange.
This cross section is only related to the mass of H ± and it is independent of the other four parameters µ 2 2 , λ 2 , m H 0 , m A 0 .It should be reminded that, near the threshold regions, it drops quickly due to the phase space suppression factor β 3 .
Lastly, as the collision energy is assumed to be 250 GeV or even higher, the intermediate Z boson is always far away from its mass shell.Hence, we have neglected the effects of decay width of Z boson in this work.

e + e − → H + H − at one loop
The 't Hooft-Feynman gauge has been used to evaluate both the weak corrections as well as the QED ones.The generic Feynman diagrams for e + e − → H + H − are drawn in Fig. 3.2, which can be put into six categories: 1. One-loop corrections to the initial state vertices e + e − V (V = γ, Z) which are purely SM, G 1 and G 2 .

One-loop corrections to the vertices
3. One-loop corrections to photon and Z boson propagators as well as γ-Z, G 12 to G 18 .
Note that the mixing Z-G 0 and γ-G 0 in the s-channel contribution was not included.Due to Lorentz invariance, such mixing would be proportional to (p 1 + p 2 ) µ , which after contracting with the initial state e + e − will be proportional to m e and thereafter will vanish.
5. The various counterterms for initial and final states and also the γ and Z propagators, γ − Z mixings are also depicted in G 22 to G 24 .
Moreover, we also add real photon emission e + e − → H + H − γ which is depicted in diagrams G 25 to G 29 .Note that diagram G 27 does not contain IR divergence but it is an important piece for electromagnetic gauge invariance.
Let M V be the matrix element of all these diagrams, the virtual corrections at NLO can be expressed as Since we have omitted the decay width of the intermediate Z boson in the leading order matrix element M 0 , at fixed order we can only take into account the real part of M V .This also means that we can neglect the imaginary parts of diagrams G 12 to G 18 , which contribute to the decay width of the Z boson at higher orders when both vector bosons V attached to the loop are Z bosons.
Calculation of the one-loop corrections will lead to ultraviolet (UV) as well as IR divergences.The UV singularities are treated in the on-shell renormalization scheme after being regularized using dimensional regularization while the IR divergences arising from the diagrams involving a photon are regularized with a small fictitious photon mass λ and have to cancel with the ones from real photon emissions.
The unbroken Z 2 symmetry prevents the mixing between the SM doublet H 1 and the inert doublet H 2 which significantly facilitates the renormalization of the IDM.The full renormalization of this model has been presented in [90].In our work, we will use the on-shell scheme developed first for the SM in [132][133][134], completed by the on-shell renormalization scheme for the inert scalar fields and their masses.V stands for generic vector bosons which could be γ, Z, or W ± ; and S could be either a Goldstone G 0 , G ± or a Higgs boson h 0 , H 0 , A 0 , or H ± .
Concerning the renormalization of the SM parameters and fields, we refer to [134,135].For the renormalization of charge, an α(m Z ) scheme is used, following the procedure in Refs.[134,135].The charge is first renormalized in the Thomson limit and then switched effectively to m Z by resumming large logarithms from light fermions, which gives with and α(0) is the coupling constant renormalized in the Thomson limit.More details about this renormalization scheme can be found in Appendix A, where the uncertainties from different schemes are also discussed.Following the notation in PDG [120], Eq. ( 22) is regarded as an "on shell" definition of the running coupling constant at the scale m Z .
For the renormalization of the inert scalar bosons, we use a similar approach as in [134,135].Because of the exact Z 2 symmetry, there is no mixing between H ± -W ∓ and H ± -G ∓ , this makes the renormalization of the scalar fields easier.The Higgs wave function renormalization constant and mass counterterm are fixed by the two following on-shell conditions : • The on-shell condition for the physical mass m H ± .
• The propagator of the charged scalar must have residue 1 at the pole mass.These requirements will lead, respectively, to where ΣH ± H ± is the one-loop renormalized charged scalar self-energy, i.e.
Here Σ H ± H ± is the one-loop unrenormalized charged scalar self energy.All the SM renormalization constants defined above can be found in [134,135].Using the condition in Eq. ( 24) and the relation in Eq. ( 25), one can prove that The explicit expression for δZ H ± is presented in Appendix B.
For the counterterms, only two are new here comparing with the SM.They are listed as follows: For the counterterms of the initial state vertices: e + e − γ and e + e − Z, counterterm of the Z boson, the photon propagators, and their mixing , they are exactly the same as in the SM and can be found in [134,135].
With such a "completely on-shell" renormalization scheme, our results are totally independent of the renormalization scale µ r .As an alternative, they now depend on the new scale introduced during the renormalization of charge (set to m Z in this work).
As mentioned before, Feynman gauge is used in our calculation.Strictly speaking, to ensure our results are gauge invariant, calculation with another gauge is needed.However this will lead to a totally different model file and/or modifications in the automatic package, which is beyond our abilities at present.Hence, we can only believe that this is automatically ensured by the theory, and only check those within our abilities (such as current conservation in real emissions).
On the other hand, as already known in the renormalization of 2HDM, improper treatment of tadpole contributions may lead to gauge-dependent results (see e.g.[136][137][138]).Fortunately, compared with the general 2HDM, the IDM is simple on this issue since there is no mixing between the two doublets.Hence, we can apply the same treatment as in the SM and avoid such an issue.
Let us now discuss the treatment of the IR divergences, which includes two parts, one is from virtual corrections, the other is from real emission.As mentioned before, IR divergences in this work are regularized by introducing a small fictitious photon mass, λ.The IR divergences in virtual corrections are present in four sources: (i) wave function renormalization of charged particles such as electrons and charged scalars; (ii) vertex corrections to e + e − γ and e + e − Z, G 1 in Fig. 3.2 with V = γ, where incoming electron and positron exchange a virtual photon with each other; (iii) in the case of e + e − → H + H − vertex corrections to γH + H − and ZH + H − : G 6 in Fig. 3.2 with V = γ where the outgoing charged scalar pair exchanges a virtual photon ; (iv) in case of e + e − → H + H − box corrections: G 20 and G 21 in Fig. 3.2 with (V, V ) = (γ, γ) or (V, V ) = (γ, Z) where incoming fermions exchange one or two virtual photons with outgoing charged scalars.
IR divergence in real emission comes from phase space integration.Besides λ, two cutoffs ∆E and ∆θ are introduced to deal with the IR singularities in real photon emission based on the two cutoff phase space slicing method [139].∆E = δ s √ s/2 defines the soft photon energy cutoff for the bremsstrahlung process.It can be viewed as the photon energy cut that separates the soft and hard regions.The angle ∆θ further separates the hard part into hard-collinear and hard-noncollinear regions.The NLO corrections are then decomposed into the virtual, soft, hard-collinear, and hard noncollinear parts as follows: Here S, HC, and HC denote the contributions of soft, hard-collinear, and hard-noncollinear parts from real photon emission, respectively.V denotes the virtual corrections including loop and counterterm diagrams.CT in the third term of rhs denotes the extra contribution arising from the structure function of incoming electron and positron.More details about the treatment of IR divergence are presented in Appendix C, as well as the checks of our results for the independence of cutoffs.
The total cross section at NLO σ NLO , is the sum of LO cross section σ 0 and NLO corrections σ 1 , namely where ∆ is the relative correction.
As described in Sec.3.1 of Ref. [140], the NLO electroweak corrections can be safely grouped into two gauge-invariant parts: 1.The "QED" part, which includes all the diagrams that contain an extra photon attached to the LO diagrams.These diagrams can be easily found by investigating the generic Feynman diagrams in Fig. 3.2 and we list them here: • G 1 when the vector boson connected with the incoming electron pair is a photon; • G 6 when the vector boson connected with the outgoing charged scalar pair is a photon; • G 8 and G 9 when the vector boson connected with the outgoing charged scalar pair is a photon; • G 19 , G 20 , and G 21 when either vector boson is a photon; • G 25 to G 29 : all real emission diagrams.
Meanwhile, the photonics contribution to the wave function renormalization of the electron and charged scalar is also grouped into this part.
2. The "weak" part, which contains all the remaining contributions.
It should be noted that the "QED" part here is only a gauge-invariant subgroup of the whole QED corrections.There are remaining QED corrections in the "weak" part.Actually, it is very hard, if not impossible in this process, to separate the whole QED part because of γ − Z mixing.The relative correction is then separated into two parts correspondingly, i.e. we have Similar to the pair production of W bosons or top quarks, the Coulomb singularity appears in one-loop corrections, which is proportional to απ/β [141,142].It gives an important enhancement to the cross section near the threshold.This effect is already well known and can be resummed to all orders.After the resummation, the LO cross section becomes where we have used the "resum."as the abbreviation of resummed in the subscript to label the quantities after the resummation.The factor |ψ(0)| 2 , which was originally obtained by Sommerfeld [143] and Sakharov [144], is given by where X = απ/β.
Beyond LO, we assume that the resummed cross section can still be written into the product of two parts, similar to Eq. (31).The first part is a factor that resums all the Coulomb singularities.Here we use |ψ(0)| 2 as this factor since we only resum the leading Coulomb singularities.The second component represents a hard part which contains no more Coulomb singularities and can be calculated order by order perturbatively.Namely, we suppose Here the label "H" denotes the hard part.Expanding the rhs of Eq. (33) to LO in α and then matching it with LO cross section σ 0 , it gives This means that the LO resummed cross section in Eq. ( 31) is correctly reproduced.Now expanding the rhs of Eq. ( 33) again to NLO in α and matching it with NLO cross section σ NLO , it gives Inserting Eqs. ( 29) and (34) into Eq.( 35) and solving σ 1 H from it, it gives Here, ∆ H represents the ratio of σ 1 H to σ 0 H .In our work, numerical analysis has confirmed that ∆ H does not exhibit any further Coulomb singularities.This implies that Eq. ( 33) works at least up to one-loop order.Then, the resummed cross section at NLO is derived as follows: Expanding σ NLO resum.to even higher in α, it gives The first three terms on the rhs are the contributions of O(α 2 ), O(α 3 ), and O(α 4 ), respectively.The third term stands for two different cases of Coulomb singularity at next-to-next-to-leading order.σ 0 H X 2 /12 corresponds to the exchange of two soft photons in the final states, while σ 1 H X/2 corresponds to the exchange of one soft photon in the final states.On the other hand, ∆ H can also be separated into its QED and weak parts, similar to what is done for ∆ in Eq. ( 30).This separation is expressed as Due to the fact that Coulomb singularity is caused by soft photon exchange in the final states, it should belong to the QED part, namely we have It has also been confirmed numerically that ∆ H,QED remains finite when the velocity of outgoing charged pair goes to zero.Up to this point, we have effectively addressed the treatment of Coulomb singularities and obtained resummed cross sections.However, in this work, we are more interested in relative corrections, not only the cross sections.Let ∆ resum.be the ratio of relative corrections after resummation (σ NLO resum./σ 0 resum.−1).Due to the fact that |ψ(0)| 2 acts as a global factor in both LO and NLO cross sections, ∆ resum. is always equal to the ratio of hard part, ∆ H . Correspondingly, we also have In order to show the effect of this resummation technique, we present in a certain case the cross section and relative corrections before and after the resummation in Table 3.1 .There are a few comments on the results in the table: • The relative weak corrections remain constant both before and after resummation, meaning that ∆ resum.,weakmaintains its equality to ∆ weak throughout.
• The resummation mostly affects the results near the threshold, i.e. the region where β is small and close to 0. When β is varying from 0.9165 to 0.0283, it is observed that the QED corrections can change from 1.628% to 42.143%.
• Before the resummation, the LO cross section can be changed drastically by the QED corrections, especially near the threshold region.After the resummation, the QED corrections are small and well controlled within a size around 1%.The LO cross section is enhanced directly by the factor |ψ(0)| 2 during the resummation, which varies from 1.0134 to 1.4918 as β decreases.While the NLO cross section is much more stable in this process.
• Throughout the entire region, the QED corrections decrease by at least one order of magnitude after resummation.This observation strongly suggests that the dominant factor contributing to the QED corrections is the Coulomb singularity term.
All of this demonstrates that the resummation is necessary and effective.In the following discussion, we will use σ or ∆ to refer to the resummed values, and for convenience, we will omit the subscript "resum."

Numerical results
In the present work, the computation of all the one-loop matrix elements and counter-terms is performed with the help of FeynArts and FormCalc [145][146][147] packages.The scalar integrals are numerically evaluated using LoopTools [148,149].The other parts are obtained with the help of FDC [150] and BASES [151].The cancellation of UV divergences are obtained both analytically and numerically.It should be noted that the model file of IDM, including all needed renormalization counter terms and renormalization constants, is obtained by manually modifying the SM model file.We have checked this model file with the output from FeynRules and confirmed their agreement.Recently, an automatic tool NLOCT [152], which can automatically generate counter terms and calculate renormalization constants, became available.It provides a model for 2HDM that can be used in one-loop QCD and EW calculation.However, compared with the general 2HDM, IDM is simple because there is no mixing between the two doublets.Hence, we use our own model file mentioned above.
In this section, we present our numerical results for our process e + e − → H + H − .It has been pointed out that the triple Higgs couplings can cause a large (can go up to 30%) radiative correction to the process e + e − → H + H − [102,105] in the 2HDM.Below we examine the radiative corrections to our process in the IDM and we expect that the electroweak radiative corrections to e + e − → H + H − would have some similarities with e + e − → H 0 A 0 [85].For the QED corrections, we will study both soft photon emission as well as the hard one.

Numerical input
We adopt the following numerical values of the physical parameters from PDG [120]: 1.The fine structure constant: α(0) = 1/137.036,α(m Z ) = 1/128.943with ∆α  In the IDM, the CP even Higgs boson h 0 is the only SM-like Higgs boson observed by the LHC collaborations, and we use m h 0 = 125.18GeV.For the other IDM parameters, we perform a scan over the whole parameter space, which includes the physical masses m A 0 , m H 0 and m H± , µ 2  2 and λ 2 parameters.We take into account all the experimental constraints as well as the theoretical requirements given in the Sec.2.2.
It is found that our numerical results are almost independent of the λ 2 parameter.Consequently, we will fix λ 2 = 2 in the next part.
In the following, we will use the α(m Z ) scheme described before to present our results.

Weak corrections in three scenarios
In our numerical analysis, we will consider three scenarios, which are tabulated in Table 4.1.They are categorized in terms of degenerate/nondegenerate between inert scalar boson masses, with/without invisible Higgs decay, and with/without DM constraints.It is thought that scenario I is the simplest one with only three free parameters and is the easiest one to demonstrate the one-loop corrections, while scenario III might be more realistic after taking into account collider experimental bounds and dark matter constraints and six benchmark points are chosen from this scenario for the LHC and future collider searches.
Our numerical analysis begins with scenario I, which has only three free parameters to consider.When λ 2 is fixed, there are only two free parameters (µ 2  2 and m S = m H 0 = m A 0 = m H± ) to vary.Three typical collision energies of future e + e − colliders, namely: √ s = 250 GeV, √ s = 500 GeV and √ s = 1000 GeV are chosen to expose the effects of new physics.It should be mentioned that when the charged scalar mass is fixed, the h 0 H + H − coupling is simply determined by the parameter µ 2 2 , as given in Eq. (5).Therefore, in order to examine how the theoretical parameters like µ 2  2 and m H ± can affect the cross section and how theoretical constraints and experimental bounds can affect allowed parameter space, we introduce five cases with different values of µ 2 2 given in Table 4.2.These values of µ 2 2 corresponding to five cases are chosen in the allowed parameter space which have passed all theoretical and experimental bounds and constraints, and we have labeled these results by the labels IDM1-5 in Fig. 4 In Fig. 4.1, total cross section and relative corrections to e + e − → H + H − as a function of the inert scalar masses m S are shown in scenario I.It is observed that the allowed ranges of m S are different when µ 2  2 takes different values.More explicitly, it is found that in the range of 100−500 GeV, the lower bounds of m S in the cases of IDM1, IDM4 and IDM5 are constrained by the signal strength of h 0 → γγ at the LHC, whereas the upper bound of m S in the case of IDM5 is constrained by the unitarity condition.For example, in the IDM5 case, m S shall not be smaller than 340 GeV, which explains why there is no IDM5 in the plot with √ s = 500 GeV and there are no IDM4 and IDM5 in the plot with √ s = 250 GeV.Obviously, since the process e + e − → H + H − is s-channel dominant, when the mass of charged scalar is fixed, the cross section becomes smaller and smaller when the c.m. energy increases from 250 GeV to 1 TeV.For a fixed c.m. energy, the cross section drops quickly when the mass of charged scalar increases to the kinematic end points.It is noteworthy that a pair of light charged scalar (say 100 GeV or so) can be copiously produced and the total cross sections could be of the order of 100 fb which would lead to more than 2 × 10 5 events at LC machines with √ s = 250 GeV which can deliver about 2 ab −1 luminosity.
[GeV] It can be read from the lower panel of Fig. 4.1 that the ratio of weak corrections in the IDM can reach −6% at √ s = 250 GeV, which is almost independent of the mass m S .In contrast, at the c.m. energy √ s = 1 TeV, the ratio starts from −10% and can go up to 30% or higher.Such a behavior can be attributed to the large corrections of the h 0 H + H − coupling.For the degenerate scenario with µ 2 2 fixed, when the parameter m S = m H ± increases from 100 GeV to 500 GeV or so, the h 0 H + H − coupling proportional to λ 3 also increases almost monotonically, especially near the end point region (say 350 GeV < m S < 500 GeV).Such a h 0 H + H − coupling can contribute to the EW corrections both linearly via the virtual corrections to (γ, Z)H + H − vertices and quadratically via the wave function renormalization of the charged scalar.When m S increases, the total cross sections decrease due to the suppression of phase space, while the ratio of EW corrections becomes larger and larger, as shown in the last plot of the lower panel.
With these crucial lessons on weak corrections learned from scenario I, we are ready to further explore the whole parameter spaces of the IDM and illustrate the behavior of weak corrections for all three scenarios, as shown in Fig. 4.2.
The upper panel of Fig. 4.2 is devoted to scenario I, where the ratio of EW corrections to e + e − → H + H − as a function of the inert scalar masses m S are shown, while the λ h 0 SS /v = −λ 3 coupling is labeled as a color bar.For the c.m. energy √ s = 250 GeV, the parameter m S can be measured in the range [80,125] GeV, and the ratio ∆σ can change in the range [−5.9%, −3.8%]; for the c.m. energy √ s = 500 GeV, the parameter m S is kinematically accessible up to 250 GeV and ∆σ can change in the range [−7.8%, 1%].For these two cases, the ratio of EW corrections is not very significant, which is due to the fact that the range of λ 3 is limited to [−3.5, 1.5].In contrast, for the c.m. energy √ s = 1 TeV, the parameter m S is kinematically accessible up to 500 GeV, which in turn implies a larger range for λ 3 ∈ [−9, 2].Then the ratio of EW corrections can be in the range [−13%, 38%].In this case, the relative corrections become substantial for m S close to the threshold production, where the absolute values of λ 3 can be quite large (say |λ 3 | ∼ 9).Consequently, the relative size of loop corrections increases and can reach 40%.
For scenario II, the theoretical parameters λ 4,5 are nonvanishing and could contribute to the vacuum stability and unitarity constraints.Therefore, when compared with scenario I, two more dimensions are added to the theoretical parameter space.Similar to scenario I, the parameter λ 3 should be constrained by the Higgs data, like the branching fraction measurement of h 0 → γγ, since a light charged scalar can have a significant contribution to this branching fraction.In the case with √ s = 250 GeV, as shown by the plots of the middle row of Fig. 4.2, due to kinematics, the range of charged scalar mass that the machine can measure is [80,125] GeV.It is found that the diphoton signal strength puts a severe constraint on the parameter λ 3 , which can only be within the range [−1, 1].Thus, the effect of the h 0 H + H − coupling is as small as the scenario I where the allowed ratio of radiative corrections can spread from −6.5% to −3%.At √ s = 500 GeV, the allowed parameter λ 3 can be in the range [−1.5, 3.5] and the radiative corrections start from −10% and can go up to 2% near the end point region of m H ± .In contrast, in the case with the c.m. energy √ s = 1 TeV, the allowed range of λ 3 becomes wider, it can change from −10 to 2, which leads to a significant enhancement for the ratio of EW corrections as shown in the plot.The h 0 H + H − coupling can contribute through the charged scalar wave function renormalization and the triangle diagrams.The ratio of EW corrections is large and can even go up to 50% near the end point (for a charged scalar mass around m H + = 500 GeV).
For scenario III, we show the ratio of EW corrections in the lower panel of Fig. 4.2, where we have taken into account all existing constraints on the parameter space of the IDM, especially the dark matter relic density, direct search constraints, and collider search of dark matter candidates.We can see that at √ s = 250 GeV the dark matter constraints shrink the allowed range of λ 3 and the ratio of EW corrections can only vary between −6% and −4.8%.At √ s = 500 GeV, the ratio of EW corrections can only be negative, and it can only change from −9% to −3%.At √ s = 1 TeV, a larger theoretical parameter space can be probed.Accordingly, a wider range of λ 3 leads to a larger ratio of EW corrections, which can spread from −15% to 30% near the end point region.
When the allowed points in the parameter space are compared for scenarios II and III, it is noteworthy that the dark matter constraints can kill almost 99% points, as demonstrated in the lowest panels of Fig. 4.2, only sporadic points are allowed for scenario III.
It should be pointed out that the ratio of EW corrections could be large and can reach 20% or higher near the end point region, like in the √ s = 1 TeV case, due to the small cross section, the production rate might also be small.For example, in the IDM5 case, when m H ± = 465 GeV, the cross section is 2 fb or so.Nonetheless, when the integrated luminosity is large enough (in the case of 2/ab, we can have 4000 signal events), there is a chance to examine this loop induced effect if we could know the model parameters precisely from other measurements.Finally, we would like to point out that the EW corrections are rather significant for high c.m. energy √ s = 1 TeV than they are for low energy cases √ s = 250 TeV or 500 GeV.This is because the contributions from the boxes are rather important for high energy.

2% 100%
Table 5.1: Benchmark points consistent with collider experiments and dark matter constraints on the relic density are proposed.Decay information of H 0 , A 0 , and H ± are also given.
In Table 5.1 we propose six benchmark points from scenario III that are consistent with current collider experiments and dark matter searches.It is important to mention that several constraints from current long-lived particles searches on the IDM exist in the literature [153,154].However, we did not take them into account in the present analysis.But as it has been discussed in [53], limits from quasistable charged particle searches can be evaded if we set an upper bound on the charged scalar lifetime of τ ≤ 10 −7 s, this implies a lower bound on the total decay width of the charged scalar of Γ H + ≥ 6.582 × 10 −18 GeV, which is respected by our BPs.
It should be mentioned that we assume that IDM alone is not sufficient to accommodate the whole dark matter content of the Universe, therefore, we only demand that the relic density of IDM is smaller than the one required by experiments.As a matter of fact, there might be other extra sectors that can contribute to the relic density of the Universe but have not been taken into account in this work.One example is the right-handed neutrino sector, another possible candidate for dark matter are the axion or axionlike particles, and there could be others as well.Coupling constants and decay information are also presented in the table.In order to examine the radiative effects of weak interactions, we have deliberately proposed three BPs, i.e.BP2, BP4, and BP6, which are close to the threshold region of certain c.m. energies of 250 GeV, 500 GeV, and 1 TeV, respectively.
There are a few comments on the features of these benchmark points.(1) According to our previous findings given in [85], it is found that both A 0 and H 0 can be the dark matter candidates.For the sake of simplicity, in these BPs, we assume that the dark matter candidate is H 0 .(2) For BP2, the charged scalar has a long life time, which can trigger the signature of a displaced vertex caused by a massive charged particle.For BP3 and BP5, the A 0 has a long life time, which can trigger the signature of a displaced vertex caused by a neutral heavy particle.While for the rest of BPs, both A 0 and H ± can promptly decay when produced.(3) For BP1 and BP4, the extra invisible decay mode h 0 → H 0 H 0 can be open, which will lead to a larger invisible decay branching for the SM-like Higgs boson.According to the CDR of CEPC [155] and a recent Monte Carlo study [156], the future Higgs factory of CEPC with √ s = 250 GeV could have the potential to determine the branching fraction of the invisible decay of h 0 down to 0.3%.Therefore, it is possible to probe BP1 at the CEPC with √ s = 250 GeV not only through e − e + → H + H − but also via the process e − e + → Zh 0 .
Another comment is on the decay widths of charged scalars.When decay widths are tiny, the narrow width approximation is appropriate, as benchmark points 1−4 and 5 are narrow enough since decay widths are only 1% of masses.For the case where decay widths are not tiny, as benchmark point 6 demonstrated, where decay width is around 4% of its mass, then the narrow width approximation might not be precise enough to describe the production H + H − and the consequent decay products.A better method to describe e + e − → H + H − → W + W − H 0 H 0 is to include the decay widths of H ± in the matrix elements, instead of producing e + e − → H + H − and then decaying H ± to W ± H 0 .Nonetheless, when its decay width is 10% smaller than its mass, our numerical calculation demonstrate, the difference is acceptable.Only when its decay width is too big, say larger than 20%, the cross section in this case should be computed carefully.Furthermore, we have explored the dependence of the decay width on the charged scalar mass, and we have found that the total decay width of the charged scalar over its mass in scenario III cannot exceed 5%.
In Table 5.2, weak corrections, QED corrections, the LO and full one-loop cross sections of these BPs are provided.Generally speaking, at √ s = 250 GeV and √ s = 500 GeV, the weak corrections are negative since the size (absolute value) of the λ h 0 H + H − coupling is small (say from -1.5 to 1.5), as shown in Fig. 4.2.Only when the charged scalar is heavy enough (say 450 to 500 GeV) and the absolute value of the λ h 0 H + H − coupling can be large (say from -10 to 2 for scenario II and from -8 to 1 for scenario III), the weak corrections ∆ weak can be positive when the masses of charged scalar are close to the threshold region when √ s = 1000 GeV, as demonstrated by BP6.It should be mentioned that the positive contribution comes from diagrams with charged scalars in the loop.From the above table, one can see that at the e + e − colliders, in BP1, BP2, BP4 and BP6, e + e − → H + H − would lead to W + W − H 0 H 0 final state, which in turn could lead to a final state with dileptons and missing energy.While for BP3 and BP5 one could have the following final states: which could give dilepton events as well as multileptons.
At the LHC, the situation is slightly different because we have the following production mechanisms for charged scalar: pp → H + H − , pp → H + H 0 , and pp → H + A 0 .For the charged scalar pair production, similar final states as for e + e − → H + H − can be obtained.While the process pp → H + H 0 (respectively, pp → H + A 0 ) would give W + H 0 H 0 (respectively, , which leads to final states with single/multileptons and missing energy.At the LHC, it should be pointed out that in order to pinpoint the electroweak corrections, the QCD corrections should first be reliably calculated, which are typically quite large. Moreover, there are two more comments on our BPs: (1) We have checked our BPs in Table 5.1 and confirmed that they can survive even in light of the new direct detection results from the LUX-ZEPLIN Experiment [157].(2) For multicomponent dark matter situations, a scaling factor is included in Ref. [56]; however, this factor is not taken into account in our work.Such a factor will ease the restrictions from direct search because we require that the relic density from IDM is not greater than the one from P lanck data.Based on the fact that our BPs have already passed all the constraints in our analysis, it is easy to conclude that they will all survive if the factor is included.

Conclusions and discussions
IDM is a simple model that can solve the problem of dark matter.Such a model predicts a charged scalar in its spectrum.The observation of such charged scalar either at LHC or at the future e + e − colliders would be a conclusive evidence of physics beyond SM.
Future e + e − colliders would provide an opportunity to measure the charged scalar cross section and its properties precisely.Following the standard renormalization scheme of the SM, in this work, we study the radiative corrections of the new physics process e − e + → H + H − in the IDM.The dimensional regularization is used to evaluate the one-loop Feynman amplitudes in the Feynman−'t Hooft gauge.We employed comprehensive on-shell scheme renormalization, which means that not only the particle masses and fields, but also the coupling constant, were renormalized using on-shell conditions.As a consequence, our predictions are totally independent of the renormalization scale µ r .Nevertheless, a new scale m Z is introduced via Eq.( 22) where the large logarithms from light fermions are absorbed into the redefinition of the running coupling constant.The running coupling constant is changed from 1/137.036 to 1/128.943correspondingly.This yields a 13% difference in the LO predictions.But this new scale dependence is greatly reduced at the NLO.
We have considered the QED corrections and checked that the IR divergences cancel when we add the virtual and real photon emissions.We have also used the resummed cross section to cure the well known Coulomb singularity.In addition, collinear divergences in our calculation appear as terms proportional to log(m e ).After including the counter term from the structure function of an electron, such divergent terms should vanish in the final result.In order to check this, we vary the mass of the electron by a factor of k from 2 −6 to 2 6 , namely m e is taken k × 0.511 MeV.It is found that the result remains unchanged when k varies.
We have examined the size of the weak corrections for three representative c.m. energies: √ s = 250, √ s = 500, and √ s = 1000 GeV.After taking into account the theoretical constraints and experimental bounds for the scenario III, we have found that the weak corrections are still sizable; i.e. the weak corrections are around −6% ∼ −5% for √ s = 250 GeV, −10% ∼ −3% for √ s = 500 GeV, and −15% ∼ +25% for √ s = 1000 GeV, as shown in Fig. 4.2.The origin of those sizable corrections is the h 0 H + H − coupling which could become large for some configuration of the parameters.We have shown that the size of the radiative corrections, are typically of the order 5−25%, which makes their proper inclusion in any phenomenological studies and analyses for e + e − colliders indispensable.
From those allowed parameter points, we have proposed six benchmark points as given in Table 5.1 for future lepton collider searches.We provide the parameters, the total decay width, as well as the branching fractions of the neutral and charged scalars.According to those BPs, we have also discussed the signature of e + e − → H + H − followed by the charged scalar decays and show that it would lead to a final state with multileptons and missing energy.At the LHC, one can get similar signatures through pp → H ± H ∓ ; H ± A 0 ; H ± H 0 production.
It was noteworthy that initial state radiation can reduce the production rate of the process e + e − → Zh 0 for √ s = 240/250 GeV by a factor of 10% or so [158], according to the simulation done with the public code WHIZARD [159].For the process e + e − → H + H − , the NLO cross section σ 1 satisfies σ ), hence the initial state and final state radiations are simply linearly summed in the matrix element of M 1 .Therefore, at the O(α) order, the initial state radiation of the process e + e − → H + H − is just the same as that of the process e + e − → Zh 0 , i.e., such a reduction in the production rate also holds for the process e + e − → H + H − for √ s = 240/250 GeV.The reduction can be around −10% in the M Z scheme, as shown in Ref. [160] for √ s = 240 GeV.While for the cases √ s = 500 GeV and √ s = 1 TeV, due to the effects of radiative return [161][162][163], it is expected that the initial state radiation can enhance the cross section of e + e − → H + H − when m H ± is small, as in the cases of BP1 and BP2.
Meanwhile, at the O(α) order, the initial state radiation and final state radiation can be treated individually.In order to cure the Coulomb singularity arising from the final state interaction, the final state radiations can be partially resumed, and the cross section can be modified as given in Eq. (37).In terms of our numerical results given in Table 3.1, the Coulomb singularity can be reliably removed.And the resummation effect is represented by the factor |ϕ(0)| 2 .Thus the final state radiation can be correctly evaluated, which can increase the total cross section by a factor from 1% to 49%.It reaches its maximum when the invariant mass of the pair of charged scalar approaches the total collision energy, as shown in Table 3.1 for √ s = 500 GeV.For higher order O(α) corrections, there exist interference terms between initial state radiation and final state radiation, appropriate treatment of these two radiations is beyond the scope of current work.

A Details about the renormalization of charge
As mentioned in the main text, we have used an "on-shell" renormalization scheme not only for the masses and fields of particles but also for the charge.By adopting this approach, our results are independent of the renormalization scale µ r .However, in order to resum the large logarithms arising from light fermions vacuum polarization, a new scale Q = m Z is introduced.Further details are presented below.
The electric charge renormalization is carried out following Refs.[134,135].The renormalization constant is obtained in the Thomson limit with the condition by imposing the condition which gives where the renormalization constants are defined as follows (bare quantities are denoted by a subscript "0"): Here ε ≡ (4 − D)/2 with D being the space-time dimension, µ is the arbitrary mass parameter introduced for bare charge, µ r is the renormalization scale of charge, e 0 is the bare dimensionless charge and e is the renormalized one.For the remaining quantities in Eq. ( 43), "Π" is defined as and denotes the transverse part of the AA(AZ) self-energy.This is just the "zeromomentum" scheme in Ref. [152].
It should be noted that in the calculation of Π(0), nonperturbative strong interaction effects cannot be neglected.In order to treat them properly, we first introduce a universal quantity: to denote the difference of Πs when the external photon is on shell and off shell, respectively. Here where k is the momentum of an off-shell external photon.At the one-loop level, it is easy to separate this ∆α into several parts according to the particles inside the photon's self energy, i.e. ∆α =(∆α) f + (∆α) b (47) where the labels "f " and "b" are used to denote the contributions from fermions and bosons.
For the convenience of discussion below, the (∆α) f part can further be expressed as (∆α) f =(∆α) ℓ + (∆α) (5)  q + (∆α) top (48) and the labels "ℓ", "q", and "top" are used to denote the contributions from leptons, light quarks and top quarks, respectively.Here "5" in (∆α) q refers to the number of quark flavors.Π can also be separated in the same way, so we will also apply the labels to Π.
Although Π (5) q (0) can be calculated perturbatively, it is unreliable due to its nonperturbative nature at low energy regions.Hence we try to exact it from data.At the scale Q = m Z , the lhs of Eq. ( 49) can be obtained from experiment measurements using the dispersion relation, i.e.
As pointed out in Sec.8.2.1 of Ref. [134] , the above renormalization of charge [labeled as α(0) scheme for convenience] is carried out at zero momentum transfer, where the relevant scale is set by the masses of fermions.They are much smaller than the relevant scales in high energy experiments.Hence, when working on processes in high energy experiments, such as this work, the large ratio of these different scales will cause large logarithms.These large logarithms can be summarized in the universal quantity introduced above, ∆α(Q).
The leading logarithms in ∆α arise from fermionic contributions and can be correctly resummed to all orders in perturbative theory by the replacement: where a label "LL" is used to denotes the leading logarithms.Meanwhile, since not only the leading logarithms but also the whole fermionic contributions are gauge invariant, we can resum the latter, i.e.
On the other hand, the large logarithms in ∆α originate from the renormalization constant δZ e (0) given in Eq. ( 43); hence they will appear wherever α appears in lowest order and can be taken into account by replacing the lowest order α by a running α(Q) as Here a subscript "lowest" is used to denote the lowest order α.This replacement can be taken as an effective renormalization of α at momentum transfer Q 2 .More details about this part can be found in Refs.[134,135].In our work, we set the relevant scale Q to m Z and resum over the contribution from light fermions, namely the running coupling constant is defined as Here the label "f ̸ = top" denotes the contributions from light fermions, which are just the sum of (∆α) lepton and (∆α) q in Eq. ( 48).This is how our α(m Z ) scheme is introduced.And it is regarded as an "on-shell" definition of the running coupling constant at the scale m Z following the notation in PDG [120].Corresponding renormalization constants can be derived easily with Eq. ( 56), which is The above discussion is based on our choice of physical parametrization, {α, m W , m Z }.Sometimes people use the Fermi constant G F as an input.Here we try to discuss the uncertainties arising from different input schemes.The relationship between G F and above parameters is already well known, which is ∆r can be calculated perturbatively.It vanishes at LO and its expression at NLO can be found in Ref. [134].Here we take G F as an observable, which is used to extract the value of certain renormalized physical parameter(s).The renormalization scheme in this case is exactly the same as the α(0) scheme.
If the input scheme is {α, G F , m Z }, then this means the mass of W boson has to be extracted with Eq. (58).However, there are some difficulties in this extraction: 1.There are new physics contributions in ∆r, which make it vary in the IDM parameter space.
2. ∆r itself also depends on m W .
Hence, we use its SM value ∆r = 0.03652 [120] as an approximation to estimate m W , which leads to It is very close to the value of m W we use in the {α, m W , m Z } scheme, which means the four inputs are precisely consistent.This indicates that the difference between these two schemes is negligible.
If the input scheme is {G F , m W , m Z }, then the situation is similar, since we have verified the consistency of the four inputs.However, in this case, similar to what we have done in Eq. ( 55), we can introduce another effective coupling constant, which is The value of α G F can be obtained directly from Eq. ( 58), with no more approximation.We replace all the α with this α G F in our calculation.This can also be taken as an effective renormalization of α, with the renormalzation constant Since ∆r also contains Π(0), it will cancel with the one in δZ e (0).Hence δZ e | G F no longer contains Π(0), and the large logarithms discussed above are also resummed in this α G F scheme. with and Here the two-point functions B 0 and B 1 are defined as with µ r being the renormalization scale and A 0 being the one-point function given by In our work, as shown in Eq. ( 44), the mass parameter µ in bare charge is replaced by µ r after the renormalization.Consequently, µ r appears instead of µ in above one-and two-point functions.

C More details about treatment on IR divergences
As mentioned in the main text, IR divergences in this work are regularized with a small fictitious photon mass λ.Meanwhile, two cutoffs, ∆E and ∆θ, are introduced to deal with the IR singularities in real correction processes based on the two cutoff phase space slicing method [139].The three-body phase space of the real correction process e + e − → H + H − γ is then divided into three parts: • Soft (S) part: Where the energy of photon E γ is smaller than ∆E.
• Hard collinear (HC) part: Where E γ ≥ ∆E and the angle between photon and the beam θ γ is smaller than ∆θ.
The full NLO corrections are then separated into four parts, as given in Eq. ( 67): Here dσ V denotes the virtual correction, including loop diagrams and counter terms from renormalization.And CT in the third term of rhs denotes the extra contribution arising from the structure function of the incoming electron and positron.
In our work, the first two parts, dσ V and dσ S , are obtained using FeynArts and FormCalc [145-147] packages, in which numerical evaluations of the scalar integrals are done with LoopTools [148,149].The other two parts, dσ HC+CT and dσ HC , are obtained with the help of FDC [150] and BASES [151].
The soft part is calculated in the soft limit.After the phase space integration of the soft photon, it can be expressed as the product of a factor and the LO result.The result of this part can be found in the literature [102,106].For completeness, we give the analytical expressions here: where Li 2 represents the known dilogarithm function.FormCalc has the option to include soft bremsstrahlung automatically when calculating the virtual corrections.We have confirmed that the results from FormCalc are consistent with the above expression.Meanwhile, by varying the value of λ inside FormCalc, we confirm the λ independence of our results.The hard noncollinear part, dσ HC , needs no more special treatment.It is obtained using traditional Monte Carlo integration techniques with the FDC package, in which BASES [151] is used for multidimensional numerical integration.
It is found that in the limit ∆θ 2 ≫ m 2 e /s, the hard collinear part given in Eq. ( 69) can be expressed in a simpler form as given below We will use the form given in Eq. ( 70) as our default formula to compute the full cross section given in Eq. (67).It is interesting to explore the difference between the results of Eq. ( 69) and Eq. ( 70), and we will explore this issue near the end of this section.One-loop radiation correction includes collinear singularities, which lead to terms proportional to log(m e ).Some of them are canceled when summing up virtual and real corrections.Some of them are absorbed into the redefinition of the running coupling constant, as shown in Eq. (56).But there are still some remaining.To deal with this, the structure function approach [164] is applied.According to the approach, the cross section of e + e − annihilation can be expressed as Here f ia (x, Q 2 ) is the so-called structure function, which denotes the possibility to find parton i with energy faction x from particle a at scale Q.In Eq. ( 71), all remaining collinear singularities are absorbed into the structure functions and resummed there.The parton-level cross section dσ ij is then free of those large logarithms.The particle-level cross section dσ e + e − is obtained by convoluting the partonic cross section with the structure functions.This factorization theorem is extended from QCD with i, a ∈ {e + , e − , γ}.Strictly speaking, what we are studying in this work is the parton-level cross section of H + H − production only.Since we are interested in the NLO corrections in the whole parameter space, this is sufficient.But if we want to study the exact cross section of production, convolution with the structure functions is needed.
In QED, all the structure functions f ia (x, Q 2 ) are perturbatively calculable.They can be expanded into a series of α, namely Meanwhile it is obvious that the first term in above expansion always takes the form Based on this, in this work up to one-loop level, Eq. ( 71) can be rewritten into (using f ee ≡ f e − e − = f e + e + ) The latter term on the rhs is just the counterterm mentioned above.
On the other hand, f ee satisfies following equation [164]: with being the regularized Altarelli-Parisi splitting function.This equation can be solved recursively according to the power of α, which leads to In this work, the scale Q is taken as the energy of beam in the center-of-mass frame, Q = √ s/2.The CT parts are finally obtained via Eqs.( 76) and ( 79) as [after same notations as in Eq. ( 69 Now all the parts in Eq. ( 67) are available.Table C.2: Data used in the checks on the independence of ∆θ in unit of fb are shown.The seventh column and the eighth column with the superscript "origin."denote the results by using Eq.(69).In contrast, the third column and the sixth column denote the results by using Eq.(70).Numbers in the brackets are integration errors to the last digits of the data.
It should be stressed that ∆E (δ s ) and ∆θ are unphysical cutoffs we introduced to deal with IR and collinear singularities.Our final results should not depend on them.In order to check this, we choose a certain point (same as the one used in Table A.2) in the IDM parameter space and compare the NLO corrections σ 1 obtained with different cutoffs.The results are shown in From the first subfigure, it can be seen that the independence on δ s is found in a wide range and δ s = 10 −3 is used as our default choice.
Collinear divergences in our calculation appear as terms proportional to log(m e ).After including the counterterm from the structure function of an electron, such divergent terms should vanish in the final result.In order to check this, we vary the mass of the electron by a factor of k from 2 −6 to 2 6 , namely m e is taken k × 0.511 MeV.The cancellation is shown in the second subfigure of Fig. C.1, from which we can see that the result remains unchanged as k varies.Furthermore, singular terms appear only in the σ V +S and σ SC parts.
In the third subfigure, we can see that the result becomes cut dependent when ∆θ is smaller than 10 −4 .It can be attributed to the fact that Eq. ( 70) can only hold when ∆θ ≫ m e / √ s ∼ 2 × 10 −6 .
It is interesting to compare the results computed by using Eqs.( 70) and ( 69), which are provided in Table C.2.It is clear that the result which is labeled by a superscript "origin." is almost independent of the values of ∆θ (when ∆θ is small enough), and the difference between the results calculated by using Eqs.( 70) and ( 69) is tiny when ∆θ is chosen appropriately (say ∆θ ∼ 10 −3 or larger).In our practical computation, we choose ∆θ = 10 −3 in this work.
In addition to the above checks for independence on λ, δ s , ∆θ, and log(m e ), our results also passed many other self-checks: • The IDM model used in FDC is generated with its own code.We have confirmed that both FDC and FormCalc give the same results at LO.
• The soft part dσ S has been calculated individually and checked with analytic results, as mentioned before.
• The other parts of real emission, dσ HC+CT and dσ HC , obtained with FDC are also treelevel calculations.This part of FDC has been used in many other calculations (see e.g.[165,166]).
• FormCalc is a public package that has been used in many one-loop EW calculations (see e.g.[167,168]) as well as in our previous works [85,160].
• We have also checked our model file using the output from FeynRules.
However, it should be noted further cross-checks are still important and welcome.

Figure 3 . 2 :
Figure 3.2: Generic one-loop Feynman diagrams for e + e − → H + H − are shown, where F stands for SM fermions;V stands for generic vector bosons which could be γ, Z, or W ± ; and S could be either a Goldstone G 0 , G ± or a Higgs boson h 0 , H 0 , A 0 , or H ± .

2 .
The gauge boson masses: m W = 80.379 GeV and m Z = 91.188GeV. 3. The fermion masses: m e = 0.511 MeV m u = 0.134 GeV m t = 173.2GeV, m µ = 0.106 GeV m d = 0.134 GeV m b = 4.660 GeV, m τ = 1.777GeV m s = 0.095 GeV m c = 1.275GeV.The masses of u and d are taken as effective parameters to reproduce ∆α

Figure 4 . 1 :
Figure 4.1: The total cross section and the ratio of EW corrections to e + e − −→ H + H − as a function of the inert scalar masses mS with three collision energies √ s = 250 GeV, 500 GeV, and 1000 GeV are shown in scenario I in the upper and lower panels, respectively.The corresponding values of µ 2 2 are given in Table4.2.

Figure 4 . 2 :
Figure 4.2: Electroweak corrections to e + e − → H + H − for √ s = 250 GeV, 500 GeV, and 1000 GeV as a function of mS and the λ hSS coupling normalized to the VEV of the SM Higgs are shown.Plots in the upper panel show scenario I, and the middle and lower panels show the nondegenerate scenario before (scenario II) and after applying dark matter constraints (scenario III), respectively.

Fig. C. 1 ,
with corresponding data in Tables C.1 and C.2.

Table 3 .
1: Cross section and relative corrections before (second block) and after (third block) the resummation of the Coulomb singularity with √ s = 500 GeV are demonstrated.The remaining IDM parameters are chosen as m H 0

Table 4 .
.1.We have deliberately introduced two with positive values and two with negative values of µ 2 2 .1: Scenarios and their conditions are tabulated.Scenario I where all the inert scalars are degenerate is described by three parameters which are (mS = m H 0 = m A 0 = mH±, µ 2 2 , λ2), scenarios II and III with a nondegenerate spectrum for the inert scalars are described by five parameters, which are (m H ± , m H 0 , m A 0 , µ 2 2 , λ2).

Table 4 .
2: In scenario I, five cases with typical values of µ 2 2 labeled as IDM1-5 are given.

Table 5 .
2: Weak corrections, QED corrections, the LO and full one-loop cross sections of BPs are provided.