Precise predictions for charged Higgs boson pair production in photon-photon collisions

The charged Higgs pair production via photon-photon collisions is investigated in the framework of two Higgs doublet model (2HDM), taking into account a full set of one-loop-level scattering amplitudes, i.e., including electroweak (EW) corrections together with soft and hard QED radiation. The numerical evaluation is carried out for three different scenarios, so-called non-alignment, low-$m_H$ and short-cascade, defined in the presence of the up-to-date experimental constraints and consistent with theoretical constraints as well. The total cross sections of $\gamma\gamma \rightarrow H^{-}H^{+}$ are scanned over the plane ($m_{\phi^0},\sqrt{s}$), where $\phi^0$ is $h^0$ for low-$m_{H^0}$ scenario and $H^0$ for other two scenarios. The regions of the parameter space in which the production rates are sufficiently large are highlighted for each scenario. The production rates in different polarization collision modes of initial beams are also discussed. It can be enhanced up to two-times by oppositely polarized photons at high energies and right-handed polarized photons at low energies. Furthermore, decay channels of the charged Higgs boson are examined for each scenario. It is observed that the one-loop EW corrections, i.e., the virtual plus real radiation corrections, reduce the tree-level cross sections and the relative correction is typically in the range of $-10\%$ to $-30\%$, depending on model parameters.


I. INTRODUCTION
The Standard Model (SM) is the most successful model of particle physics so far, because its theoretical predictions are magnificently compatible with experiments. The electroweak symmetry breaking in the SM is successfully accomplished through the Higgs mechanism, which yields their masses to the fundamental particles and suggests the existence of a Higgs boson. A huge step was taken for high energy physics with a 125 GeV of Higgs boson observed by both CMS and ATLAS experiments at the LHC [1 -3]. Although the properties of the Higgs boson discovered are suitable with the SM predictions, it is possible that it constitutes but one state of a richer Higgs sector. In association with that, the models with extended Higgs sector are among the best motivated beyond SM models (BSM), as they may provide solutions to many deficiencies in the SM, such as the gauge hierarchy problem, the origin of dark matter, the generation of a baryon asymmetry and the strong CP problem.
One of the simplest such extensions is the two-Higgs doublet model (2HDM) [4,5] that only adds one more Higgs doublet to the SM under the same (SM) gauge symmetry. Its scalar sector consists of five scalars, one CP-odd Higgs (or pseudoscalar) A 0 , two CP-even Higgs bosons (h 0 and H 0 ) and two charged Higgs H ± . Versions of the 2HDM are often used in several of well-founded new physics BSM scenarios, both with and without Supersymmetry [6], where the additional Higgs doublet is either an essential by-product or a necessary component in indicating issues. * mehmetdemirci@ktu.edu.tr Since it has an ultimate point as a high-precision experimental machine and hence its complementary with LHC, it is useful to investigate the phenomenology of the Higgs sector in detail in the context of an electronpositron collider. In this respect, the clean environment in these colliders would ensure that the Higgs sector is precisely identified. International Linear Collider (ILC) [7,8] is an efficient machine for precise experiments which is aimed to yield equipment for electron-positron collisions along with other possibilities such as electronelectron, electron-photon and photon-photon collisions. The Compact Linear Collider (CLIC) [9] is a TeVscale high-luminosity linear collider operating with the centerof-mass energy energy up to 3 TeV. The primary task of these linear colliders is to expand and complement the results obtained in the hadron colliders, and to discover new physics BSM. Besides, the photon-photon-collision is of the opinion as other collision mode which can yield an integrated luminosity of the order of one hundred fb −1 per year. This collision can produce unique new physics compared to other types of collisions. Upgrading the machine is expected to reach the high energies √ s = 1 TeV with up to three hundred fb −1 per year of total integrated luminosity [10]. In addition to the possible discovery of the Higgs boson relatives through examining the observed-Higgs boson properties, the ILC presents great opportunities to explore new lighter Higgs bosons or, more commonly, any weakly interacting light (pseudo)scalar boson via the direct production [11]. Discovering the charged Higgs boson would become a clear proof of physics BSM and an eminent signal for extended Higgs sector. An extensive review on charged Higgs phenomenology is available in Ref. [12]. The future e − e + and γγ colliders with high energy and high luminosity have significant potential in the discovery of the charged Higgs boson and in the study of its properties. Furthermore, high energy γγ-collisions, where photons directly coupled to charged particles, can give a better understanding of the SM and its extensions for several aspects. Correspondingly, exact predictions are needed for the physical observables related to charged Higgs bosons.
The main channels of the pair production for charged Higgs bosons in the linear colliders are e − e + → H − H + and γγ → H − H + . Since the cross-section of e − e + → H − H + is suppressed by s-channel contributions, the production rate of γγ → H − H + is larger than that of e − e +collision mode. The scattering process e − e + → H − H + has been widely studied by including the one-loop corrections in the framework of both 2HDM and MSSM [13][14][15][16][17][18]. On the other hand, the scattering process γγ → H − H + has been studied at the one-loop level by including the full squark corrections [19] and the real radiative corrections [20,21] in the framework of SUSY but only yukawa corrections in the 2HDM [22]. Full one-loop level corrections for this production mode should be considered in the framework of the 2HDM. In order to benefit from the high precision measurements, we also need a high precision predictions from the theory, which means that there is a need to go beyond the leading order calculations for most processes, hence the full one-loop contributions become significant for physics analyses at the future colliders, e.g. ILC or CLIC. The scattering process γγ → H − H + have rich physics results and needs a detailed study in the light of current constraints.
In the present work, the production of the charged Higgs bosons pairs through photon-photon collision is investigated in the framework of 2HDM for the first time, including the full one-loop contributions, i.e. electroweak (EW) corrections, as well as hard and soft QED radiation. For numerical evaluation, six different benchmark points, which have a CP -even scalar with mass of 125 GeV and couplings compatible with those of the observed Higgs boson, are chosen. They are constructed from the scenarios so-called as "non-alignment", "low-m H " and "short-cascade" [23]. These scenarios are defined in the presence of the up-to-date experimental constraints and consistent with theoretical constraints as well. The total cross sections of γγ → H − H + are scanned over the plane (m φ 0 , √ s), where φ 0 is h 0 for low-m H 0 scenario and H 0 for other two scenarios. The parameter space regions in which the improvement of the production rate is significant enough to be accessible at the future colliders are highlighted for each scenario. The production rates in different polarization collision modes of initial beams are also discussed.
The rest of this work is organized as follows. In Section II, a brief review of the 2HDM is presented. In Sec. III, we review the theoretical and experimental constraints imposed on the 2HDM. Then, the benchmark scenarios used in the calculation are given. In Sec. IV, the Feynman diagrams, the corresponding amplitudes, and some useful analytical expressions are given for the considered scattering process. The general forms of the virtual and the real photon radiation corrections are also discussed. In Sec. V, we present the numerical results related to the scattering process and decay channels, and discuss in detail the corresponding model parameter dependencies of the cross sections. Finally, the concluding remarks are presented in Sec. VI.

II. REVIEW OF THE TWO HIGGS DOUBLET MODEL
In this section, a brief overview of the 2HDM with CP-conserving, concerning only relevant details for this study, is introduced. The 2HDM is constructed by augmenting the complex scalar doublet, Φ 1 , of the SM by another doublet, Φ 2 , which changes the dynamics of EW symmetry-breaking. The most general scalar potential being invariant under the SU(2) L ⊗ U(1) Y gauge group is defined by (2.1) where λ i (i=1,...,7) are quartic coupling parameters 1 and Φ 1,2 are two complex scalar Higgs doublets. To comply with some low energy observables, the discrete Z 2 symmetry is put forward by the Paschos-Glashow-Weinberg theorem [24,25]. In particular, this symmetry is applied to prevent flavor changing neutral currents at tree-level. The Z 2 symmetry requires to λ 6,7 and m 2 12 be zero. However, allowing m 2 12 be non-zero, the Z 2 symmetry are softly broken. The charges assigned under the Z 2 symmetry allow that each type of fermion interaction with only one Higgs doublet, i.e., Φ 1 or Φ 2 . According to the Z 2 assignment, there appear four 2HDM-types, which are mostly categorized as type-I, II, III (or "lepton-specific") and IV (or "flipped") of 2HDM [4,5]. Table I shows of how fermions couple to each Higgs doublet (Φ 1,2 ) in the allowed types where flavor conservation is naturally obeyed. This work is concentrated on the Type-I and Type-II of 2HDM. In Type-I, only the doublet Φ 2 interacts with both leptons and quarks like in the SM. In Type-II, the doublet Φ 1 couples to leptons and d-type quarks, while Φ 2 couples to u-type quarks.  I  Φ2  Φ2  Φ2  II  Φ2  Φ1  Φ1  III  Φ2  Φ2  Φ1  IV  Φ2  Φ1  Φ2 After the spontaneous symmetry breaking of the SU(2) L ⊗U(1) Y gauge symmetry associated with the electroweak force, the neutral components of scalar doublet acquire vacuum expectation values v j such that where ρ j and η j are real scalar fields. The combination v 2 = v 2 1 + v 2 2 ≃ (246 GeV) 2 is set by its relation to the mass of W and the Fermi constant as follows: v 2 = 1/( √ 2G F ) = 4m 2 W /g 2 . These Higgs doublets include initially eight degrees of freedom. The three of them, G ± , G 0 bosons, are eaten by the longitudinal components of the EW vector bosons Z and W ± . The remaining ones are five physical Higgs bosons: two CP-even h 0 and H 0 , a CP-odd Higgs A 0 , and two charged scalars H ± . They are related to the weak eigenstates via with the generic rotation matrix 2 For any given value of t β , the parameters m 2 1 and m 2 2 are determined by the minimization conditions of potential. The mass parameters m 2 1,2 and quartic couplings λ 1 -λ 5 can be defined in terms of the physical masses m h , m H , m A , m H ± , along with t β = v 2 /v 1 (the ratio of vacuum expectation values), and the mixing term of neutral sector s β−α . The soft Z 2 symmetry breaking parameter m 2 12 is defined by where the last equality is only for the tree level. Fixing λ 6 and λ 7 to zero to respect the Z 2 symmetry, m 2 12 , tan β, mixing angle α and four Higgs boson masses are enough to specify the model completely in the physical basis 3 . Consequently, there are seven independent free parameters encountered in the Higgs sector of the 2HDM. The phenomenology of the 2HDM depends significantly on the size of the mixing angle α together with angle β. There appears an alignment limit [26], where the CP-even Higgs boson h 0 looks SM-like Higgs boson if s β−α → 1 or c β−α → 0. The alignment limit is the most favoured condition by the experimentalists. In this limit, the couplings of the CP -even Higgs boson h 0 in 2HDM are similar to Higgs boson of SM and can be identified as the discovered 125 GeV Higgs boson. Furthermore, the CP -even Higgs boson H 0 acts as gauge-phobic such that its coupling to the vector bosons Z/W ± is greatly suppressed. However, in the limit c β−α → 1, H 0 looks SM-like Higgs boson.
Furthermore, a decoupling limit appears when c β−α = 0 and m H 0 ,A 0 ,H ± ≫ m Z [27]. At this limit, the Higgs boson h 0 coupling to SM particles completely appear like the couplings of the SM-Higgs boson that contain the coupling h 0 h 0 h 0 .

III. PARAMETER SETTING ON 2HDM
The parameter space of 2HDM is subjected to both the bounds coming from experimental searches and theoretical constraints that arise from the model itself. These have to be imposed to the free parameters of the model.

A. Theoretical and experimental constraints
The parameter space of the scalar 2HDM potential is reduced by the theoretical constraints: potential stability, perturbativity and unitarity. The V 2HDM is bounded from below respect to the vacuum stability of the 2HDM. Namely, V 2HDM ≥ 0 must be satisfied for all directions of Φ 1 and Φ 2 . Consequently, the following conditions are placed on the parameters λ i [27][28][29]: Another set of constraints enforce that the perturbative unitarity (for details, see Refs. [30,31]) need to be fulfilled for scattering of longitudinally polarized gauge and Higgs bosons. Besides, the scalar potential needs to be perturbative by demanding that all quartic coefficients satisfy |λ 1,2,3,4,5 | ≤ 8π. Furthermore, note that the global fit to EW measurements dictates ∆ρ to be O(10 −3 ) [32,33]. This forbids the occurrence of large mass splitting between Higgs bosons of 2HDM, and imposes that m H ± ∼ m A or m H or m h .
Besides the above theoretical constraints, 2HDMs have been researched in the past and still ongoing experiments, such as direct observations at the LHC or indirect B physics observables. Consequently, many results have been accumulated since then, and the parameter space of the 2HDM is restricted by all conducted results. In the Type-I of 2HDM, the following pseudoscalar Higgs mass regions 310 < m A < 410 GeV for m H = 150 GeV, 335 < m A < 400 GeV for m H = 200 GeV, 350 < m A < 400 GeV for m H = 250 GeV with t β = 10 have been excluded by the LHC experiment [34]. Furthermore, the CP-odd Higgs mass is bounded as m A > 350 for t β < 5 [35] and the mass range 170 < m H < 360 GeV with t β < 1.5 is excluded for the Type-I [36].
The charged Higgs boson mass is subject to a number of limits from direct experimental researches at the LHC (and previous colliders) as well as from various B-physics observables. In type-II and IV of 2HDM, the BR(b → sγ) measurement places bounds on the charged Higgs mass as m H ± > 580 GeV for t β ≥ 1 [37,38]. On the contrary, in the type-I and III of 2HDM, this bound is much lower [39]. Considering t β ≥ 2, there is a possibility for the charged Higgs bosons in type-I and III of 2HDM to be as light as 100 GeV [39,40] while at the same time compatible with LEP and LHC limits as well as with constraints of all B physics [41][42][43][44][45][46]. Additionally, according to LHC, Tevatron and LEP results [47], there is no exclusion about s β−α = 1 for m A,H,H ± = 500 GeV in the Type-I 2HDM.

B. Benchmark points scenarios
Let us now provide details of the benchmark scenarios used in this study. The three different benchmark scenarios, which are named non-alignment, low-m H and short-cascade, proposed in Ref. [23], are used to investigate the phenomenology of the charged Higgs boson. All of these scenarios contain a CP -even scalar with 125 GeV mass and its couplings consistent with the successfully observed Higgs boson. In addition, a significant portion of their parameter space is allowed by the constraints from the extra Higgs bosons searches. The four benchmark points (BPs) are selected from these scenarios as shown in Table II, and they are agree with experimental and theoretical constraints. For each BP, the quantities of potential stability, perturbativity and unitarity have been verified by using 2HDMC 1.7.0 [48,49].
These benchmark scenarios are constructed on the hybrid basis 4 where the input parameter is designated as {m h , m H , c β−α , tan β, Z 4 , Z 5 , Z 7 }, for the case of softly broken Z 2 -invariant 2HDM. Here, the parameters Z 4 , Z 5 , Z 7 are quartic couplings in the Higgs basis of O(1). The masses of the pseudoscalar Higgs and charged Higgs boson in this basis are given by • In the non-alignment scenario, the discovered Higgs boson is interpreted as the lightest CP -even scalar h 0 , with SM-like properties. In the so-called alignment limit, the c hV V coupling become the same values as in the SM. In this case, the heavier CP-even Higgs boson H 0 could not decay into the gauge bosons W − W + and Z 0 Z 0 . However, this scenario is defined with a non-alignment (c β−α = 0) as permitted by the current constraints. This leads to some interesting phenomenology of the H 0 . The other two Higgs bosons H ± and A 0 are decoupled as m h 0 = 125 GeV < m H 0 < m H ± = m A 0 . In this scenario, to obtain values of m H ± which satisfy the b → sγ constraint, quartic coupling parameters are set as Z 4 = Z 5 = −2. Consequently, tan β and m H 0 remain as free parameters. One benchmark point, BP1, is selected from this scenario by fixing c β−α = 0.1 with type I.
• The short-cascade scenario is established with a SM-like h 0 by taking exact alignment, c β−α = 0. The mass hierarchy could be arranged such that it is possible for either one or both of the decay channels H 0 → W ± H ∓ or H 0 → Z 0 A 0 to be open, and resulting Higgs-to-Higgs decays in a small cascade. In this scenario, two of the Higgs masses m A 0 , m H 0 , m H ± are choosen to be equal for simplicity. These mass degeneracies could be arranged by choosing Z 4 and Z 5 properly. In this study, two different mass hierarchy are considered for Type-I 2HDM. They are as follows: m A 0 < m H ± = m H 0 for Z 4 = −Z 5 = −1, and m H ± < m A 0 = m H 0 for Z 4 = 2 and Z 5 = 0, while keeping the remaining parameters fixed such as tan β = 2, Z 7 = −1. Two benchmark points, BP2 and BP3, are taken from this scenario by fixing c β−α = 0 with type-I.
• Another scenario is the low-m H scenario where both CP -even Higgs bosons (h 0 , H 0 ) are light, however; the heavier one (H 0 ) is assumed as SMlike Higgs boson so that m H 0 = 125 GeV. The coupling of the heavier CP-even Higgs to gauge bosons is proportional to c β−α . Because m h < m H , the couplings of lighter CP-even scalars to vector bosons must have been heavily suppressed to agree with direct search bounds which forces s β−α → 0. Similarly to non-alignment scenario, to decouple the Higgs bosons A 0 and H ± , the quartic parameters Z 4 and Z 5 are taken as Z 4 = Z 5 = −5. The

Scenario
BPs The parameter space for 90 < m h < 120 GeV is restricted by the searches h → bb, τ τ at the LHC, leads to an upper bound on tan β [23]. Therefore, it is taken as tan β = 1.5. The BP4 is selected from this scenario by fixing c β−α = 1.0 with type II.

IV. THE CROSS SECTION OF THE CHARGED HIGGS BOSON PAIR PRODUCTION IN PHOTON-PHOTON COLLISIONS: THEORETICAL SETUP
The future e − e + colliders could also supply another possibility to measure the production of H − H + in photon-photon collision mode. This can be carried through the Compton backscattered photons, produced in the intense laser photon scatterings on the initial electron-positron beams. Similar with the study on the scattering process e − e + → H − H + , the evaluation of the corresponding pair production in photon-photon collisions is also crucial for the precise experimental measurements of production of H − H + at an e − e + linear collider.
Here, some details about the calculation of the treelevel and one-loop corrections of cross section are respectively provided. The relevant production channel is expressed by where the quantities p 1 , p 2 , k 1 and k 2 in parentheses label the corresponding particle 4-momenta. λ 1 and λ 2 denote the helicities of initial photons and can take the value of either +1 or −1, corresponding to a right-handedly (R) and left-handedly (L) polarized photon beam, respectively. All these momenta obey the on-shell equations p 2 1 = p 2 2 = 0 and k 2 1 = k 2 2 = m 2 H ± . Also note for further use the Mandelstam variables: The polarization vectors of the photons are introduced by which ensure ε i · p j = 0 for i, j = 1, 2. The analytical and numerical evaluation have been obtained by adopting Mathematica packages 5 as follows: The Feynman diagrams and the relevant amplitudes have been created with the help of FeynArts [53]. Then, the squaring the corresponding amplitudes, the simplification of the fermion chains, and the numerical evaluation have been carried out by FormCalc [54]. The scalar loop integrals have been evaluated by LoopTools [55]. The phase-space integration is computed by the Monte-Carlo integration algorithm Vegas as implemented in the CUBA library [56].

A. Leading Order Calculation
At tree-level, the leading contribution to the process occurs via t and u-channel charged Higgs-exchange diagram and the quartic coupling diagram. These contributions are at an order of O(α ew ) and they are based on pure electroweak interactions. The tree-level Feynman diagrams contributing to γγ → H − H + in the 2HDM is given in Fig. 1. The matrix elements corresponding to each diagram are respectively expressed as where ε µ (p 1 ) and ε ν (p 2 ) are polarization vector of incoming photons. The total amplitude at the lowest order could be determined by summing up all the above matrix elements: The tree-level total amplitude includes only the couplings which are independent of the 2HDM angles. These couplings are universal in a sense and come from the kinetic energy term of the Higgs fields. However, the other couplings, After squaring the total amplitude and summing over the helicities of the final states, the total cross-section σ(γγ → H − H + ) is obtained by using following formulâ where the average over initial photons spins represented by the factor (1/4). The parameters entering the treelevel cross-section are all standard except for the mass of the charged Higgs. Moreover, the non-standard parameters of 2HDM appearing at the one-loop level could be consistently taken as bare in the calculations.

B. NLO Corrections
In the analysis of high-energy processes observed at the current and future colliders, higher-order corrections (at least, next-to-leading order contributions) should be included for precise theoretical predictions. The process (4.1) has one-loop level contributions as the nextto-leading order. They are based on pure EW interactions at the order of one-loop. In the one-loop level, the total amplitude can be expressed as a linear sum of triangle, box, and bubble one-loop integrals. According to the type of loop correction, the virtual contributions are coming from three different type diagrams: self-energy, boxtype, and vertex-type (triangle&bubble-type s-channel) diagrams.
A complete set of one-loop Feynman diagrams for γγ → H − H + in the 2HDM and the corresponding amplitudes have been provided by the FeynArts. These diagrams are drawn in Figs. 2 to 4. In the other set of diagrams, particles in each loop are running counterclockwise. Note that Feynman diagrams with electron-Higgs couplings have been neglected. The internal particles in diagrams are labeled as follows: φ 0 indicates to all neutral Higgs/Goldstone bosons (h 0 , H 0 , A 0 , G 0 ); u ± indicates the ghosts; u m /d m can be u/d-type quarks (the subscript m represents the generation of quark) and l stands for leptons e, µ, τ . In loop diagrams, dashedlines indicate to neutral and charged Higgs bosons, and wavy-lines represent γ and Z, W ± -bosons. The Mandelstam variables (4.2) are also valid for in there. The contribution can also be topologically divided intoŝ,t andû-channel diagrams with intermediate the neutral The self-energy diagrams consist of all possible self-energy loops of quarks, gauge bosons, and neutral/charged Higgs/Goldstone bosons on the propagator of charged Higgs boson as shown in  Fig. 4. Theŝ-channel diagrams may be make a significant contribution to the cross section, however they are nearly negligible away from the mass pole of the propagator.
The relevant total amplitude can be written by summation over all contributions from self-energy, triangle, and box diagrams as (4.11) For virtual one-loop corrections, the differential cross section, summing over the helicities of the final states, can be calculated by where |M virt | 2 is not included since it is very small. The virtual contributions are ultraviolet (UV) and infrared (IR) divergent. The UV divergences are handled by dimensional regularization [57] in the on-mass-shell renormalization scheme. The counter terms are included via diagrams in Fig. 5. For 2HDM, all Feynman rules including counter terms and the renormalization conditions are described in Ref. [58]. Here, the calculation has been  carried out by using the FeynArts model file that includes these counter terms (for details see Ref. [58]). The calculation of the amplitude has been performed in 't Hooft-Feynman gauge. After the renormalization procedure, the virtual part becomes UV-finite. Though, it still contains the soft IR singularities originated from the exchange of virtual photons in the loops. These singularities are regularized by introducing a photon mass parameter, m γ 6 . All these singularities in the limit m γ → 0 are cancelled by adding the real photon bremsstrahlung corrections, according to the Kinoshita-Lee-Nauenberg theorem [59,60]. The real photon radiation process is denoted by where k 3 is the four-momenta of radiated photon. The relevant diagrams are given in Fig. 6. According to the 6 This is automatically carried out by LoopTools.
energy of the radiation photon k 0 3 = | − → k 3 | 2 + m 2 γ , the bremsstrahlung phase space is divided into a soft region and a hard region. Hence, the real photon radiation correction is written as follows: where δ s is the soft cut-off energy parameter δ s = ∆E γ /( √ŝ /2). If the energy of radiation photon is k 0 the radiation photon is hard. The soft part is calculated by using the soft photon approximation formula [61,62] where dσ 0 is the tree-level differential cross section and the soft photon cut-off energy ∆E γ satisfies k 0 3 ≤ ∆E γ ≪ √ŝ .  Although both soft and hard terms depend on soft cutoff parameter δ s , the real correction does not depend on the soft cut-off parameter. Furthermore, summing the virtual and soft contributions drops out the dependence of the IR regulator m γ . The result now depends on the parameter δ s , i.e., ∆E γ , and the hard photon radiation contribution must be added as well for dropping this dependence out.
Consequently, the UV and IR finite total one-loop corrections are expressed as a sum of the virtual, the soft photon radiation, and the hard photon radiation: which is independent of the IR regulator m γ and soft cut-off parameter δ s . We have numerically checked that our results do not depend on m γ or on ∆E γ = δ s √ŝ /2. For the repre-sentatively non-alignment scenario with t β = 10 and m H 0 = 150 GeV, the virtual plus soft correction, the hard photon radiation correction and the total one-loop correction are plotted as a function of the soft cutoff parameter δ s at √ŝ = 1 TeV in Fig. 7. As one can see from Non-alignment so that where F γ/e (x) represents the photon structure function. The photon spectrum is qualitatively better for larger values of the x-fraction of the longitudinal momentum of the e − -beam. In the case of x > 2(1 + √ 2) ≈ 4.8, the high-energy photons could vanish via the pair production of e − e + in its collision with a subsequent laser-γ. The energy spectrum of the photon provided as a Compton backscattered photon off the e − -beam [63] has been utilized for F γ/e (x) in this study.

V. NUMERICAL RESULTS AND DISCUSSIONS
The numerical results of the production of the charged Higgs boson pairs through photon-photon collisions are discussed in detail, considering full one-loop corrections in the 2HDM, including soft and hard QED radiation. For each benchmark scenario, the tree-level and the NLO corrections of the cross sections are numerically evaluated as a function of the center-of-mass energy and the mass of Higgs boson selected as a non-fixed free parameter. The regions of the parameter space where the production rates are large enough to be detectable are highlighted. The longitudinal polarizations of the initial beams are significant to improve the production rate; therefore, some polarization distributions are presented, as well. Decay channels of the charged Higgs boson are also investigated for the scenarios interested.
The following notations are used here: i. σ LO := the tree-level total cross section.
ii. σ NLO := the NLO corrections of the cross-section i.e., the virtual plus real contributions.
iii. σ LO+NLO := the full cross section including all oneloop corrections.
In the following subsections, the numerical evaluations for each BP is separately presented. The numerical results are of course dependent on the choice of the 2HDM parameters. Nonetheless, they provide an opinion of the relevance of the full one-loop contributions. As a general comment, it can be pointed out that the tree-level cross sections depend solely on the parameters of SM (and m H ± ). As a result, any dependence on the free parameters of 2HDM can appear firstly at the one-loop level (except for m H ± ).

A. Non-alignment scenario
For BP1 defined in the non-alignment scenario, the tree-level and the NLO corrections of cross sections of , respectively. Additionally, to describe the full EW corrections to the tree-level cross section quantitatively, the corresponding relative correction is plotted as functions of m H 0 -√ s in Fig. 8(c). The scan parameters are varied as follows: 150 ≤ m H 0 ≤ 600 GeV in steps of 10 GeV, and 850 ≤ √ s ≤ 3000 GeV in steps of 50 GeV. Here, it is assumed that the lightest CP-even scalar is the SM-like Higgs boson. For this scenario, the vacuum stability and perturbativity do not allow a large split between m A 0 and m H ± ; for this reason it is set m A 0 = m H ± . Also, any value of m A 0 is allowed via the EW precision tests. Consequently, this scenario have a mass hierarchy as follow: We can see from these figures that the total cross sections are sensitive to the mass of Higgs boson m H 0 since m H 0 is directly related with m H ± so the phase space of the final state particles. In the case of m H ± > √ s/2, the charged Higgs pair production is kinematically inaccessible as shown by white regions in parameter space. When m H ± ≪ √ s/2, the cross section nearly scales as 1/s and reaches its larger values. When the mass of charged Higgs boson m H ± with along m H 0 becomes larger, the tree-level cross section decreases as expected while NLO corrections slowly increases. Particularly, the cross section reaches its larger values for m H 0 < 400 GeV in the scan region. The NLO corrections make negative contributions to total cross section except for extreme points in the parameter regions. The NLO correction of total cross section σ NLO (γγ → H − H + ) ranges from -7 to -2 fb at most of the parameter space. The relative correction is always negative in the whole region and mostly ranges from −15% to −30% as seen from the contour lines in Fig. 8(c). Its magnitude increases with increasing of √ s γγ and reaches about −31% at √ s γγ = 3 TeV. However, the relative correction becomes larger near the production threshold ( √ s γγ ≈ 2m H ± ) because the cross section is very small in this region, so this enhancement is phenomenologically insignificant. At √ s γγ = 1 TeV, the relative correction is about −13.06% for m H ± = 401.25 GeV and −8.38% for m H ± = 492.63 GeV. For example, at √ s γγ = 900 GeV for m H 0 = 150 GeV (m H ± = 379.05 GeV), the full one-loop corrected cross section is σ LO+NLO (γγ → H − H + ) = 51.97 fb with δ r = −11.28%. Overall, the full one-loop corrected cross section is at a visible level of O(10 1 fb) in the range of 4 to 55 fb for considered parameter regions of non-alignmentscenario.
In Fig. 9, the initial beam polarisation dependence of the integrated tree-level and full one-loop EW-corrected cross sections are plotted as a function of √ s γγ for BP1, where we take m H 0 = 200 GeV and t β = 10. The √ s γγ varies from the value little larger than the threshold 2m H ± to 3 TeV. The curves correspond to the integrated cross section with oppositely-polarized photons (+−), right-handed polarized photons (++) and unpolarized photons (UU), respectively. Note that the integrated cross sections with the (+−) and (−+) photon polarization are equal: σ +− = σ −+ . As indicated in this figure, all curves reach to their maximum cross section values as the center-of-mass energy goes from the threshold value to the corresponding position of peak. This is also the expected behavior. These peaks appear at √ s γγ = 890 GeV for σ UU , √ s γγ = 880 GeV for σ ++ and √ s γγ = 1600 GeV for σ +− . At high energies, the integrated cross sections with oppositely polarized photons (−+) or (+−) are enhanced by a factor of 1.9 as compared to the unpolarized case. In the case of both photons with right-handed or left-handed polarized (++) or (−−), although at high energies, the integrated cross sections are highly suppressed, at low energies they are amplified up to about 2 times. These results imply that having both photons polarized can turn out to be significant to ensure a measurable production rate. The σ UU When m H ± ≪ √ s/2, the cross section approximately scales as 1/s, the t-channel contributions to production rate become important. When m H 0 runs from 250 to 500 GeV, m H ± varies from 250 to 500 GeV for BP2, while it varies from 48.75 to 436 GeV for BP3. Due to presence of exact alignment c β−α = 0 in the short-cascade scenario, the coupling C h 0 H − H + would reach its largest value affecting the cross section. However, the effect of C H 0 H − H + on the total cross-section will be reduced.
For both benchmark points, the NLO corrections make negative contributions to total cross section except for threshold points in the parameter space. The tree-level cross section decreases with increment of m H ± as expected, while the NLO corrections increases quantitatively. The size of NLO-corrected cross sections reach up to a level of 10 2 and 10 3 fb for BP2 and BP3, respectively.
As shown in Fig. 10, the σ NLO (γγ → H − H + ) for BP2 ranges from -12 to -2 fb at the considered parameter space, and its maximum value of -11.93 fb reaches at m H ± ∼ m W ± + m A 0 . The corresponding relative correction mostly ranges from −10% to −35% as seen from the contour lines in Fig. 10(c). Its magnitude increases with increasing of √ s γγ and reaches around −35% at √ s γγ = 3 TeV in the region of 250 GeV ≤ m H ± ≤ 320 GeV. Note that the relative correction becomes positive in a small region due to the production threshold ( √ s γγ ≈ 2m H ± ). For example, at √ s γγ = 1 TeV, the relative correction increases from −16.59% to −7.13% when m H ± running from 250 GeV to 490 GeV. The NLO corrected cross section σ LO+NLO (γγ → H − H + ) reaches a maximum value of 130.46 fb with δ r = −7.23% at √ s γγ = 550 GeV for m H ± = 250 GeV. At √ s γγ = 1.5 TeV, when m H ± running from 250 to 500 GeV, the σ LO+NLO (γγ → H − H + ) decreases from 29.99 to 19.16 fb while δ r changes from −22.63% to −18.27%.
In Fig. 11, the initial beam polarisation dependence of the integrated tree-level and full one-loop EW-corrected cross sections are plotted as a function of √ s γγ for BP2, where we take m H 0 = 250 GeV. The curves here are like those in Fig. 9. All curves increase firstly, reach their maximal values, and then decrease with the increment of √ s γγ . The σ UU LO and σ UU LO+NLO have a peak around  decrease from −5% to −35% while δ ++ r become largest where the production cross section σ ++ goes to zero. When √ s γγ running from 520 to 3000 GeV, the polarization improvement varies from 0.007 to 1.98 and 2.00 to 0.02 with oppositely polarized photons (+−) and right-handed polarized photons (++), respectively, compared to the unpolarized case.
For BP3, the σ NLO (γγ → H − H + ), as seen from Fig. 12, ranges from -34 to -3 fb at the considered parameter space, and its maximum value of -34.34 fb reaches at √ s γγ = 150 GeV and m H 0 = 250 GeV. The corresponding relative correction mostly ranges from −5% to −30%. Its magnitude increases with increasing of √ s γγ and reaches about −30% at √ s γγ = 3 TeV in the region of 173 GeV ≤ m H ± ≤ 317 GeV. Note that the relative correction becomes positive in a small re-gion due to the production threshold ( √ s γγ ≈ 2m H ± ). For example, at √ s γγ = 1 TeV, the relative correction decreases from −8.8% to −14.11% when m H 0 running from 250 GeV to 400 GeV. At √ s γγ = 0.5 TeV, the σ LO and σ LO+NLO decrease from 436.8 to 111.9 fb and 419.4 to 107.9 fb, respectively, as m H ± running from 48.75 GeV to 243 GeV. The NLO corrected cross section σ LO+NLO (γγ → H − H + ) reaches a maximum value of 2.34 pb with δ r = −1.45% at √ s γγ = 150 GeV for m H ± = 48.75 GeV. At √ s γγ = 1 TeV, when m H 0 running from 250 to 500 GeV, the σ LO+NLO (γγ → H − H + ) decreases from 111.89 to 40.28 fb while δ r varies from −8.8% to −12.13%.
In Fig. 13, the initial beam polarisation distributions on σ LO and σ LO+NLO are given as a function of √ s γγ for BP3, where we take m H 0 = 300 GeV. The curves here are labeled the same as in Fig. 9. All curves increase firstly, reach their maximal values, and then decrease with the increment of √ s γγ . The σ UU LO and σ UU LO+NLO peak around √ s γγ = 380 GeV with a value of 300.14 fb and 290.77 fb, respectively, with δ UU with δ +− r = −9.37%. The σ ++ LO and σ ++ LO+NLO peak at √ s γγ = 380 GeV with a value of 585.02 fb and 565.51 fb, respectively, with δ ++ r = −3.09%. For BP3, both δ UU r and δ ++ r decrease from −2% to −29% while δ ++ r become largest where the production cross section σ ++ goes to zero. At high energies, the integrated cross sections with oppositely polarized photons are enhanced by a factor of 1.99 as compared to the unpolarized case. In the case of both photons with righthanded polarized, the integrated cross sections are highly suppressed at high energies; but at low energies they are amplified up to 1.99 times. The polarization considerably improves the production rate, as expected. This improvement is almost independent of scenarios considered in this study. At other scenarios, similar improvement appears.

C. Low-mH scenario
In Fig. 14, the tree-level and the NLO corrections of cross sections, and the corresponding relative correction of process γγ → H − H + are scanned over the regions of m h 0 -√ s in the low-m H scenario, where the heaviest CP-even Higgs H 0 behaves like the SM Higgs boson and its mass is taken as m H 0 = 125 GeV. The scan parameters are varied as follows: 65 ≤ m h 0 ≤ 120 GeV in steps of 1 GeV, and 1150 ≤ √ s ≤ 3000 GeV in steps of 50 GeV. Similar to the non-alignment scenario, the two Higgs bosons, A 0 and H ± , decouple sufficiently such that they do not affect the phenomenology. The corresponding mass hierarchy is m h 0 < m H 0 = 125 GeV < m H ± = m A 0 . The values of the mass of charged Higgs boson m H ± , which are calculated in terms of the mass of neutral Higgs boson m h 0 , are shown with contour lines, and the m H ± changes very slowly with increasing of the m h 0 . When m h 0 runs from 65 to 120 GeV, m H ± varies from 554 to 563 GeV for BP4. Due to being c β−α = 1 in the low-m H scenario, the coupling C H 0 H − H + would reach its largest value affecting the cross section. However, the effect of C h 0 H − H + on the total cross-section will be reduced.
The integrated cross sections σ LO (γγ → H − H + ) and σ NLO (γγ → H − H + ) change very slowly with m h 0 , due to the small range of mass m H ± . They decrease with increasing of √ s γγ when m H ± ≪ √ s/2. Particularly, they reaches its larger values for √ s γγ < 1500 GeV in the scan region. The NLO corrections make negative contributions to total cross section as in the previously discussed scenarios. The σ NLO (γγ → H − H + ) ranges from -3.5 to -3 fb at most of the parameter space. On the other hand, the relative correction mostly ranges from −10% to −32% as seen from the contour lines in Fig. 14(c). Its magnitude increases with increasing of √ s γγ . For example, at √ s γγ = 1.2 TeV, the relative correction increases from −9.85% to −9.53% when m h 0 running from 65 GeV to 120 GeV. The NLO corrected cross section σ LO+NLO (γγ → H − H + ) reaches a maximum value of 25.47 fb with δ r = −11.1% at √ s γγ = 1250 GeV for m H ± = 554 GeV. Overall, the σ LO+NLO (γγ → H − H + ) appears usually in the range of 6 to 25 fb for considered parameter regions of the low-m H scenario. Figure 15 presents the initial beam polarisation dependence of the integrated tree-level and full one-loop EW-corrected cross sections versus √ s γγ for BP4, where we take m h 0 = 65 GeV. The √ s γγ varies from the value little larger than the threshold 2m H ± to 3 TeV. The curves here are labeled the same as in Fig. 9. It is seen that all curves increase firstly, reach their maximal values, and then decrease with the increment of √ s γγ . The are enhanced by a factor of 1.7 as compared to the unpolarized case. The σ ++ LO and σ ++ NLO are significantly suppressed for the region of √ s γγ > 2.5 TeV ; but for the region of 1.1 < √ s γγ < 1.5 TeV they are amplified by between 1.5 and 2 times. It can be seen that the longitudinal polarization of initial photons increases the production rate of H − H + signal in the photon-photon colliders.

D. Angular Distribution of the Differential Cross Section
The tree-level and the virtual plus soft photon corrected of differential cross sections of γγ → H − H + are presented in Fig. 16 as a function of the angle between the initial photon and the charged Higgs boson at √ s γγ = 1.5 TeV. Also, in the same figure, the corresponding relative corrections as a function of cos θ are shown on the bottom panel for each BPs. It can be seen in this figure that all curves are rather symmetric according to cos(θ) = 0. The tree-level differential cross sections relatively flat, particularly in the region of −0.6 < cos θ < 0.6, but the virtual plus soft photon corrections significantly depend on the angle θ. The corrections reach their maximum values when cos θ have values of +1 and -1. Namely, the charged Higgs pairs are dominantly produced in the backward and forward directions and it will be much more possible to detect them in that region of the collision. On the other hand, the differential cross-section is smaller in large values of the charged Higgs mass. When cos θ goes from 0 to +1 or -1, the relative correction δ r varies from −25.58% to −17.01% for BP1, −32.48% to −17.21% for BP2, −29.67% to center-of-mass energy. Also, in the same figure, the corresponding relative correction as a function of √ s is shown on the bottom panel. From this figure we find the expected behavior: a rapid increase near the production threshold, followed by a decrease with the increment of the colliding e − e + center-of-mass energy. It is obvious that the production rate in BP3 is always larger than those in the other BPs.

F. Decay channels of the charged Higgs boson
The final decay products of the produced charged Higgs bosons will be analyzed for all scenarios in this section. The decay channels are calculated by using 2HDMC 1.7.0. To explore the process in a collider, we must firstly identify all the possible charged Higgs products. The total decay widths of charged Higgs boson H ± versus the mass of Higgs boson, m φ 0 where φ 0 is h 0 for low-m H 0 scenario and H 0 for other two scenarios, are plotted in    In Fig. 19, we show the branching ratios of charged Higgs boson for the dominant decay modes as a function of m h 0 ,H 0 for all scenarios discussed in the text. The mass of charged Higgs boson as a function of m h 0 ,H 0 is also shown by star symbols on the right axis. Note that the decay modes for which BR value is less than 10 −4 are omitted here for clarity. When considering only decays to SM particles, the dominant decay modes are those involving the heaviest lepton or quark pairs accessible such as follows: τ ν τ , tb, ts and cs. Among these, the decay mode tb of the charged Higgs boson is available in all BPs significantly with varying branching ratios. In scenarios with Type I, due to the cot β dependence of coupling, the Br(H ± → τ ± ν τ ) is suppressed by m 2 τ /m 2 t over Br(H ± → tb).
In the non-alignment scenario, the dominant decay mode of charged Higgs is in the W ± H 0 channel, following sub-dominant channels tb and W ± h 0 for BP1, followed by other suppressed modes such as, H ± → ts and cs, particularly in the range m H 0 < 500 GeV. Then, Br(H ± → W ± H 0 ) gradually decreases at larger values of m H 0 , and the decay mode W ± h 0 becomes dominant. Br(H ± → W ± h 0 ) increases from 1.2 to 66.3% in the mass window 150 GeV < m H 0 < 600 GeV. For m H 0 < 500 GeV, the process becomes γγ → H + H − → W + H 0 W − H 0 . On the other hand, the dominant decay mode of H 0 is W + W − with branching ratio of 88.7 − 50.2%. When the hadronic decays of W ± -boson are considered, there appear 12 jets at the final state. Consequently, it is difficult to reconstruct the W ± , so the charged Higgs bosons.
In the short-cascade scenario, once the decay H ± → tb opens up (when m H ± > m t + m b ), it quickly becomes dominant, leading almost 100% Br. On the other hand, the decay modes W ± A 0 and ts (τ ± ν τ and cs ) get suppressed for m H 0 > 300 GeV in the BP2 (BP3). In the short-cascade scenario, the process becomes γγ → H + H − → tbtb, and the decays of t may be an ideal option for reconstructing the process at m H 0 > 300 GeV. The subsequent decays of t → W b, W → qq(lν l ) will form the signature of H ± at a detector. Consequently, it can be tagged with 8-jets plus 2-b-tagged jets for the short-cascade scenario.
In the low-m H scenario, the decay channel H ± → W ± h 0 is clearly dominant over the full mass range (65 GeV < m h 0 < 120 GeV), because its decay width is proportional to c β−α leading it dominant (∼ 100%) for the choice of s β−α = 0. The sub-dominant channel for BP4 is H ± → tb with around 18% Br. The Br(h 0 → bb) changes between 90.0 − 87.7%. Therefore, the process can be tagged with 4-jets plus 4-b-tagged jets.

VI. SUMMARY AND CONCLUSIONS
The 2HDM is the simplest extension of the SM which contains the charged Higgs boson. In the case of a discovery of charged Higgs boson, a subsequent exact measurement of its properties will be important for determining its nature and the corresponding model parameters. In order to provide enough precision, full one-loop contributions need to be included in the production channels of charged Higgs boson. The pair production of charged Higgs is one of the main channels that would provide an observable signal in a wide range of the parameter space in 2HDM. In this study, the charged Higgs pair production has been studied via γγ collisions, considering a complete set of one-loop EW corrections in the framework of 2HDM. In the one-loop diagrams, the UV divergences have been regularized by dimensional regularization in the on-mass-shell renormalization scheme, and IR divergences have been canceled by the inclusion of soft and hard QED radiation. The numerical evaluation was carried out for three different scenarios, so-called non-alignment, low-m H and short-cascade, defined in the framework of 2HDM, in the presence of the up-to-date experimental constraints. The tree-level and full one-loop EW corrections of total cross sections have been scanned over the plane (m φ 0 , √ s), where φ 0 is h 0 for low-m H 0 scenario and H 0 for other two scenarios (non-alignment and short-cascade scenarios). The regions of the parameter space in which the production rates including the relative one-loop corrections are sufficiently large have been highlighted for each scenario.
The results show that the one-loop EW corrections mostly reduce the tree-level cross section and the relative correction is typically few tens of percent for both the γγ → H − H + and the e − e + → γγ → H − H + depending on chosen parameter space. The virtual plus the real corrections are mostly negative for selected BPs. The overall effect range between −10% and −30% in a wide range of the model parameters. Since the tree-level cross-section is mostly from QED, the model-dependent parameters appear firstly at one-loop level. The cross section in the short-cascade scenario is always larger than those in the other two scenarios and the cross section of the low-m H 0 scenario is the smallest one among all of the three scenarios. The production rate with light charged Higgs bosons is larger than that with heavy charged Higgs bosons owing to the larger final state phase space volume for either an electron-positron or photon-photon collider. The full one-loop corrected cross sections of γγ → H − H + reach their maximum values as follows: σ BP1 LO+NLO = 51.97 fb with δ r = −11.28%, σ BP2 LO+NLO = 130.46 fb with δ r = −7.23%, σ BP3 LO+NLO = 2.34 pb with δ r = −1.45% and σ BP4 LO+NLO = 25.47 fb with δ r = −11.1%. The absolute relative corrections increase with the increment of √ s for all cases. The production rates of γγ → H − H + in different polarization collision modes of initial beams have been also discussed. The production rate of γγ → H − H + is enhanced up to around 2 times with oppositely polarized photons at high energies and right-handed polarized photons at low energies, as independent of the scenarios interested. Consequently, having both photons polarized can turn out to be significant to ensure a measurable production rate.
The reconstruction of the charged Higgs boson has been presented for each scenarios, studying its dominant decay modes. In the non-alignment and low-m H scenarios, the bosonic decay channels H ± → W ± H 0 and W ± h 0 are dominant, respectively, while in the short-cascade scenario, the dominant decay channel is H ± → tb. The bosonic decay channels are highly suppressed due to alignment limit and limited phase space.
In summary, the first phenomenological results in the context of 2HDM for the one-loop EW corrections to the charged Higgs pair production via photon-photon collisions have been produced, and in the light of this, the main distinctive features between the selected scenarios have been highlighted. The precise measurements for the associated production would be possible at the future colliders, and our results will be helpful for determining new physics signals based on the 2HDM and putting more precise limits on the model parameters.