Charged Higgs-boson pair production in association with the $Z^0$ boson at electron-positron colliders

In this paper, the production of the charged Higgs pair associated with the $Z^0$ boson is analyzed in the minimal extension of the standard model the so-called two-Higgs-doublet model (2HDM). The process $e^+e^- \rightarrow H^+H^-Z^0$ is calculated at the tree level including all the possible diagrams in 2HDM. The numerical analysis is performed in consideration of the current experimental constraints and various scenarios for the free parameters of the model. The results are presented as a function of center-of-mass energy, the charged Higgs mass ($m_{H^\pm}$), and the ratio of the vacuum expectation values ($t_\beta$). The unpolarized cross section, taking into account the results in the flavor physics, gets up to $0.278\text{ fb}$ for $m_{H^\pm}=175\text{ GeV}$, and it declines with decreasing $m_{H^\pm}$ in Type-I. However, it gets down to $0.073\text{ fb}$ for $m_{H^\pm}=500\text{ GeV}$ in Type-II. Further, the calculation is also carried out in the non-alignment scenario and low-$m_{h^0}$ scenario. The effect of the polarized incoming $e^+$ and $e^-$ beams shows that the cross section is enhanced by a factor up to 2.5 at P(+0.60,-0.80) polarization configuration. Decay channels of the charged Higgs, possible final states of the process, and some differential distributions belong to the charged Higgs and $Z^0$ boson are examined for each scenario. The analysis shows that some channels have higher branching ratio such as $H^+\rightarrow t \bar{b}$, $H^+\rightarrow W^+h^0$, and $H^+\rightarrow W^+H^0$. These decay channels are essential for the charged Higgs searches in the lepton colliders regarding the scenarios interested. The detection of the charged Higgs is a powerful sign for the extended scalar sectors, and the results show the potential of a future lepton collider.


I. INTRODUCTION
Over the last couple of decades, many extensions to cure the quadratic divergence at the scalar sector of the Standard Model (SM) has been proposed, and the implications of the new physics have been studied intensively. One possible extension of the SM is to add a second Higgs doublet to the scalar sector. These two Higgs doublets defined to have the same quantum number so that they together could give mass to leptons, quarks, and electroweak bosons. Addition of an extra doublet gives new couplings and interactions, in a result rich phenomenology to the 2HDM. In a general 2HDM, there are two charged Higgs bosons (H ± ) and three neutral Higgs bosons (h 0 , A 0 , H 0 ) [1,2] playing with the free parameters of the model h 0 could be set to resemble the discovered Higgs boson.
Nowadays, there is an ongoing effort for another project named Linear Collider Collaboration (LCC), it is an organization that brings the two most likely candidates for the next collider program, the Compact Linear Collider Study (CLIC) and the International Linear Collider (ILC), together under one roof. When this project is constructed, there e + e − , e − e − and γe collisions will be studied. One of the primary task at the future lepton colliders is to complement the LHC results, and also searching for clues in beyond the SM such as supersymmetry, an extension of the scalar sector or exotic models. Both of the collider projects are designed to study the properties of the new particles and the interactions they might undergo according to the vast amount of theories. As it is known, the lepton colliders compared to the LHC have a cleaner background and it is possible to extract the new physics signals from the background more easily.
There has been a long time effort to observe a hint associated with a charged Higgs boson in the past and current experiments. However, it was not discovered at the LEP, Tevatron, and yet the search is still going on at the LHC [3]. The main discovery channel of the charged Higgs boson in 2HDM is through e + e − → H + H − or e + e − → H ± f f channels.
Such signals in 2HDM are studied in Ref. [4,5]. On the other hand, the pair production of the charged Higgs boson associated with gauge boson is also possible at the lepton colliders where this process could be studied with more precision. The process e + e − → H + H − Z 0 is complementary to the discovery of the charged Higgs pair production. The process would also give a way to measure the couplings c H + H − h 0 and c H + H − H 0 . These couplings along with the other trilinear Higgs couplings will help us to reconstruct the Higgs potential [? ]. The process is investigated before in Ref. [6] for a case study where the triple Higgs couplings are studied, and the diagrams which are sensitive to the triple Higgs couplings are included in the calculation. The same process is investigated in left-right twin Higgs model where doubly charged Higgs pair production with Z 0 boson is analyzed [7], and a similar process in the Higgs triplet model is examined in Ref. [8]. In this work, the production of the H + H − Z 0 including all the possible diagrams at the born level is calculated. Numerical analysis of the total cross section is performed as a function of the center-of-mass (COM) energy, the charged Higgs mass, and the free parameters of the 2HDM for the benchmark scenarios with various motivations. In addition to these, the results with different beam polarizations are presented. A discussion is carried on the couplings essential in the production of H + H − Z 0 and decay channels of the charged Higgs boson for each scenario. A preliminary calculation of the differential cross section as a function of the kinematical properties of the Z and the H ± boson is obtained. Each scenario with subsequent decay channels of H ± is explored for later Monte Carlo studies.
The content of this paper is organized as follows. In Sec. II, the scalar sector, free parameters of the model and the masses of the Higgs bosons in 2HDM are reviewed. The machinery, the workflow of the analysis, and the kinematics of the scattering process are explained in Sec. III. The constraints coming from the experimental results are underlined in Sec IV.
Numerical analysis is performed on three scenarios, which are named as non-alignment scenario, low-m H mass scenario, and the favored region in light of the flavor physics results.
The numerical results for each of the scenario, including the beam polarization are presented thoroughly in Sec V. Decay channels of the charged Higgs and identification of the process are discussed along with the differential distributions in Sec VI. The conclusion is drawn in Sec VII.

II. SCALAR SECTOR AND THEORETICAL FRAMEWORK FOR 2HDM
In this section, the phenomenology of the 2HDM, the scalar sector, and the free parameters in the model are presented. The model itself and the detailed introduction of the framework are studied before by many authors, and it is given in Ref. [2,[9][10][11]. Therefore, only a short review of the 2HDM is presented which is relevant to the analysis. 2HDM is constructed by adding a second SU (2) L Higgs doublet with the same a hypercharge (Y = 1) to the scalar sector. If we denote the Higgs doublets as with the scalar potential of the 2HDM given in Eq. 2, there will be in total 14 free parameters.
Both of the doublets have the same charge assignment, and they could couple to leptons and quarks as in the SM. However, to suppress the CP violation and the flavor changing neutral currents (FCNC) at the tree level, the construction of the model needs to be constrained in some way. Traditionally, a discrete symmetry (Z 2 ) was introduced which puts restrictions on the most general form of the Higgs scalar potential and the Higgs-fermion interactions [12][13][14]. Discrete Z 2 symmetry is simply defined as the invariance of the Lagrangian under the interchange of Φ 1 → Φ 1 and Φ 2 → −Φ 2 (in a generic basis). If the discrete Z 2 symmetry is extended to the Yukawa sector, Higgs-fermion (Yukawa) interactions could be written in a couple of four different and independent ways. Since we do not care the CP violation, we set λ 6 = λ 7 = 0. Then the complex parameters m 12 and λ 5 are taken as real. As a result, the free parameter number reduces to 8 under these assumptions. If this symmetry is allowed to violate softly, then FCNCs are naturally suppressed at the tree level. The m 2 12 term in Eq. 2 ensures the breaking of the discrete symmetry softly.
The mass spectrum of the Higgs bosons is computed first by imposing the constraints obtained from the potential minimum condition (∂V /∂Φ i = 0, for i=1,2) and eliminating m 11 and m 22 . Second, the Higgs doublets are defined as given in detail in Ref. [1,9,15,16]. After decomposing the scalar potential into a quadratic term plus cubic and quartic interactions, the mass terms are extracted. Finally, diagonalizing the quadratic terms, we easily obtained the physical Higgs states and their masses. The masses of the charged and the CP-odd Higgs states are defined as follows: where β is the ratio of the vev of the Higgs doublets (tan β = v 1 /v 2 ). The mass of the CP-even states becomes where R is a unitary rotation matrix which diagonalizes the CP-neutral Higgs mass matrix [1,10] as a function of the angle (β − α), GeV. Then, the physical neutral CP-even scalar states are obtained by orthogonal combinations of ρ 1 and ρ 2 given in Eq. 1. A lighter h 0 and a heavier H 0 bosons are defined as h 0 = ρ 1 sin α − ρ 2 cos α and H 0 = ρ 1 cos α − ρ 2 sin α. Accordingly, the SM Higgs boson would be where the angle which rotates the CP-even Higgs states are defined as s βα = sin(β − α), c βα = cos(β − α). If s βα = 1 is assumed, which is called the SM-alignment limit [17][18][19][20], an important feature shows up; the ratio of the couplings between h 0 (H 0 ) and the SM gauge bosons (V = W ± /Z 0 ) to the corresponding SM Higgs one will be , respectively. Therefore, the lighter CP-even Higgs boson (h 0 ) becomes indistinguishable from the Standard Model Higgs boson, and H 0 acts as gaugephobic (c βα = 0) in this limit. In the literature, to explore the phenomenology of heavier CP-even boson (H 0 ), it is custom not to set s βα to unity. That is also explored, and small deviations from unity are considered in the numerical calculation. The rest of the degree of freedom makes the prominent property of the model; two charged Higgs bosons and three neutral Higgs bosons [2]. In conclusion, the free parameters of the model are the masses of the neutral Higgs bosons (m h/H 0 /A 0 ) and the charged Higgs bosons (m H ± ), the ratio of the vacuum expectation values (t β = v 2 /v 1 ), the mixing angle between the CP-even neutral Higgs states (α), and the soft breaking scale of the discrete symmetry m 12 [17].

III. CALCULATING THE CROSS SECTION
In this section, the analytical expressions, the vertices and the Feynman diagrams relevant to the scattering process e + e − → H + H − Z 0 are presented. Throughout this paper, the process is denoted as where FeynRules [21,22]. If there is a scalar particle in the model, there is a possibility of a scalar mediator between the incoming and outgoing states as it is seen in Fig. 1. Then, these schannel diagrams make a significant contribution, but they are almost negligible away from the mass pole of the mediator. In any case, the narrow-width-approximation for the scalars Another way to suppress the FCNCs in the Yukawa sector is to impose a natural condition which is to take the two Yukawa matrices to be aligned (Eq. (2.7) in [24]), and the Z 2 symmetric types are assumed as the particular cases of this aligned model. There, the ζ L is the factor defines the Yukawa coupling structure namely the Type-I through -IV. As follows, ζ L = 1/t β in Type-I and Type-IV, and it is defined as ζ L = −t β in Type-II [24] and -III. For that reason, the numerical results presented in this work hold for Type-III and -IV as well.
Having defined all the couplings, the amplitude for each of the diagrams are constructed using FeynArts [25,26]. Next, the simplification of the fermion chains, squaring the corresponding amplitudes, and the numerical analysis is accomplished using FormCalc [27] routines.
After squaring the amplitude, a summation over the polarization vectors of the final states and averaging over the helicities of the initial states is performed. Finally, the total  angle is defined as s w = sin θ w , c w = cos θ w . s αβ = sin(α+β), and c αβ = cos(α+β) are abbreviated.
unpolarized cross section is defined as an integral over a 4-fold differential cross section [28] σ = 1 4 µ,ν where the energy of the k 3 and k 5 are defined below and a = √ s − k 0 5 .
One other option at the LC is to polarize the incoming beams which could maximize the physics potential, both in the performance of precision tests and in revealing the properties of the new physics beyond the SM. In this study, we also explored the dependence of the cross section on the polarization of the incoming electron and positron beams. The polarization of the incoming beam is significant especially to enhance some of the helicity channels.
In annihilation diagrams, typically in s-channel, the helicities of the incoming beams are coupled to each other. For that reason, the helicities of the incoming particles in the SM need to be opposite from one another to recombine into the vector boson mediator, the Z 0 boson or the photon. In exchange diagrams, the helicities of the incoming beams are directly coupled to the helicities of the final states. In this case, all helicity configurations for the beams are in principle possible.
The expression for the cross section for an arbitrary degree of longitudinal beam polarization is defined as: where σ LR denotes the cross section for 100% left-handed positron and 100% right-handed electron polarization. P e − and P e + denote the percentage of the electron and positron beam polarization, respectively. The σ LL and σ RR configurations are omitted due to the negligible contributions in Eq. 11. In this study, we examined various polarization configurations in an LC and presented in the next chapters. Two of these configurations are inspired from the ILC which in pol-1 (pol-2) is defined as right-handed positron with 30% (60%) polarization [29,30] and an electron with 80% left-handed polarization on both cases, respectively.

POINTS
There are two sets of constraints in 2HDM; one is the theoretical constraints which come from the theory itself, the other one is the results coming from the measurements carried out in the past and current experiments. These constraints need to be applied to the free parameters defined in the previous section.

A. Theoretical Constraints
• Perturbative Unitarity: This requirement comes from the fact that the scattering amplitudes need to be flat at the asymptotically high energies. Due to the additional Higgs states in 2HDM, we need to make sure that Higgs-Higgs and Higgs-Vector boson scattering cross sections are bounded by 16π [31].
• Perturbativity: The theory needs to be inside the perturbative region. The perturbative region is defined as where all the quartic couplings in the theory are small, and we take them to be |λ i | ≤ 4π.
• Vacuum Stability: The scalar potential defined in Eq. 2 needs to be positive in any direction of the field space, even at the asymptotically large values [32][33][34][35][36]. Stability constraint translated into the following conditions:

Experimental Constraints
• ATLAS and CMS experiments at the LHC have reported the so long hunted resonance with a mass of 125.4 ± 0.4 GeV [37][38][39][40][41][42]. We know that the announced peak must have a nature of a CP-even, and for that reason 125.4 GeV peak needs to correspond to the one of the CP-even h 0 or H 0 states in 2HDM. In the numerical calculation, we considered both of these possibilities and investigated the implications on the cross section of e + e − → H + H − Z 0 . If we assume that h 0 is the discovered particle at the LHC, that puts a restraint on the couplings. Therefore, s β−α is pushed to unity. In the opposite case, where the H 0 is assumed that it is the discovered one then c βα is set to unity. However, we let that factor to deviate from unity just for the sake of phenomenological curiosity.
• 2HDM needs to be compatible with the full set of electroweak precision observables which are measured in the previous experiments [43]. There are parameters S, T, and U which are called oblique parameters [44,45], and they merely represent the radiative corrections to the two-point correlation functions of the electroweak gauge bosons.
These parameters are sensitive to any new physics contribution, and they take the value of the top-quark and Higgs masses so that they are set to vanish for a reference point in the SM (S = T = U = 0). According to that, a sizable deviation from zero would be an indicator of the existence of new physics. These oblique parameters have been calculated by [46][47][48][49] for the scenarios presented in the next section with the help of 2HDMC, and in all cases the oblique parameters are much less than 10 −2 .
• The LEP experiment excluded the charged Higgs boson with a mass below 80 GeV (Type II scenario) or 72.5 GeV [50] (Type I scenario, for pseudo-scalar masses above 12 GeV) at the 95% CL. If it is assumed that B(H + → τ + ν) = 1, then charged Higgs mass bound increased to 94 GeV for all tan β values. [51]. The Tevatron experiments D0 [52][53][54] and CDF [55] excluded the charged Higgs mass in the range of 80 GeV < m H ± < 155 GeV at the 95% CL assuming B(H + → cs). The search on charged Higgs is also carried out at the LHC in the decay of top quark [56,57], and upper limits are set on the B(t → H + b) and B(H + → τ ν). More recent results are given in [58] and the references therein.
• It is known that charged scalar states in 2HDM affects the flavor physics, particularly B s → X s γ or B s → µµ. In general, the flavor observables in these models are sensitive to the m H ± and tan β. According to Refs. [59,60],B −B mixing disfavors tan β < 0.5, and also lets the couplings of the Higgs to heavy quarks to be in the perturbative region t β < 60 [61,62]. More discussion is carried out in Ref. [50] and the references therein.

C. Benchmark Points
Recognizing all these experimental and theoretical constraints, the following scenarios and benchmark points are chosen. These points are preferred by aiming at a broader survey of the region of phenomenologically interested. As it is stated in the experimental constraints, we employed the 2HDMC [23] to check whether the theoretical constraints for each benchmark point are fulfilled.
• Non-alignment scenario : This benchmark point is taken from Ref. [63] where the particle with the mass of 125 GeV is interpreted as the lightest CP-even Higgs boson (h 0 ) with SM-like couplings. In the so-called alignment limit where |c βα | → 0 the c hV V coupling approaches the corresponding SM value, and in which case the heavier CP-even Higgs boson H 0 could not decay into W + W − /ZZ. The alignment limit is also endorsed by the additional Higgs boson searches at the Large Hadron Collider.
However, to allow some interesting phenomenology for the heavier CP even state (H 0 ), this benchmark is defined with a non-alignment (c βα = 0) as allowed by the present constraints. With this motivation, the mass of the other Higgs states A 0 and H ± are taken as degenerate and they are let to decouple Detailed analysis is carried out in [63] and the charged Higgs mass is given in the Hybrid base as where Z 5 is the quartic coupling parameter, and it is set to −2 along with Z 4 . In this scenario, t β and and m H 0 are taken as a free parameters and their ranges, as well as the other parameters, are given in Tab. III. • Low-m H mass scenario : As it is known, there are two CP-even states (h 0 /H 0 ) in 2HDM. In the previous scenario, it is assumed that the discovered scalar particle at the LHC in 2012 is the lighter CP-even Higgs (h 0 ) state. However, there is one other possibility that it could be the heavier CP-even Higgs (H 0 ) state. In that case, the coupling of the heavier CP-even Higgs to gauge bosons (c HV V ) will be scaled by a factor of c βα instead of s βα , and as in the previous case that also forces us to c βα ≈ ±1.
Due to the direct search limits, couplings of h 0 to vector bosons need to be suppressed heavily and even close to zero (s βα → 0). This scenario is analyzed in detail in Ref. [63] where the region 90 < m h < 120 GeV is not rejected by the LHC constraints (from h 0 → bb, τ τ ), which puts an upper limit on t β . According to the authors, it is also possible to set exact alignment (c βα = 1) with either Type-I or Type-II Yukawa couplings. Therefore, choosing a non-aligned value for c βα doesn't make much impact on the production of the cross section, but more on that is delivered in the results. In Tab. IV, the free parameters are depicted in the Hybrid base where Z 4 = Z 5 is taken so that in particular the T oblique parameter can not receive sizable contributions (Eq. (76) in Ref. [63]).  • Favored region in light of the recent experimental constraints: Taking into account all the recent updates particularly coming from the flavor physics, as the last scenario, we explored the region inspired by the results presented in Ref. [24]. Since the charged Higgs boson can contribute to flavor observables via the charged currents, flavor observables are quite significant. Motivated by the Ref. [24], the same masses for all the Higgs bosons are set m H 0 = m A 0 = m H ± , and that also satisfies the theoretical constraints. The s βα is set approaching the unity which guarantees that all the light Higgs self-couplings are close to the SM ones. In that case, the heavier CP-even Higgs boson (H 0 ) can not decay into W + W − and ZZ pairs. According to the authors, in Type-I, t β < 1 is strongly constrained by B(B 0 s → µ + µ − ) and the mass difference of the scalars. However, the mass of the scalars is not constrained on the large t β range compared to the Type-II. The cross section in Type-II, as in the Type-I, is calculated by setting the same masses for all the extra scalars. s βα is taken to unity, and due to the dominant constraints from B(B → τ ν) and B(B 0 q → µ + µ − ) high t β are excluded. The parameter region is defined in Tab. V, and more detailed discussion is given in all these constraints, the numerical analysis is carried out for unpolarized and polarized incoming beams.  , m H 0 = 425 GeV, and Type-II is assumed.

A. Non-Alignment Scenario
The cross section of the process (e + e − → H + H − Z 0 ) is computed for the parameters presented in the non-alignment scenario [63]. In Fig. 2, we plotted the cross section as a function of COM energy in 1 < √ s < 3 TeV range and for two distinct charged Higgs masses (m H ± ). We also explored the polarization configurations of the incoming beams.
As it is seen in Fig. 2, the cross section σ RL , where it represents the totally right-handed polarized e + -beam and totally left-handed polarized e − -beam, is always greater among other polarization cases. In Fig. 2 ( Fig. 2 (left). A large separation between m H ± and m A 0 is strongly constrained by the vacuum stability and perturbativity, for that reason, the calculation is carried in the limit of m A 0 ∼ m H ± in the non-alignment scenario. Therefore, any value of m A 0 is allowed by the EW precision tests.
In Fig. 2 (right), the production rate of H + H − Z 0 is plotted for various polarization configurations, and the Type-II is set. The unpolarized cross section is σ U U = 0.068 fb at √ s = 2.9 TeV. However, the polarized cross sections are σ pol-1 ≈ 0.137 fb and σ pol-2 ≈ 0.167 fb. Two-dimensional analysis of the cross section of the process is drawn in Fig. 3 as a function of the COM energy and the charged Higgs mass (m H ± ). In Fig. 3 (left), the nonalignment scenario with Type-I Yukawa couplings are set, on the (right) the same analysis with Type-II couplings are displayed. The charged Higgs mass is computed by Eq. 12, and we scanned in the allowed range 379 < m H ± < 690 GeV for this scenario. According to the couplings given in Tab. I, the m H ± dependence appears only in the CP-even to charged Higgs Higgses and fermions, these couplings do not affect the calculation even varying at this range. Moreover, even if we lose the exact alignment and set s βα = 0.9, the change in the cross section is around -2.5% at √ s = 1.5 TeV, and it falls at high COM energies. Finally, we get σ pol-1 ≈ 0.143 fb and σ pol-2 ≈ 0.175 fb for the polarized incoming beams. The Type-II Yukawa coupling structure is set in Fig. 6 (right), and due to the restrictions from the flavor observables charged Higgs mass is constrained from below (m H ± > 500 GeV).
The last result is given in Fig. 7 where the cross section at √ s = 3.0 TeV plotted as a function of t β and m H ± . On the left side, the Type-I Yukawa structure is set, and on the right, the type-II Yukawa structure is used. In Type-I, the cross section gets high for small charged Higgs masses, and it gets low as usual for high m H ± values. However, as it can be seen clearly in both of the Fig 7, t β does not affect the cross section. That means couplings given in Tab. I which have t β dependence indirectly through s β , c β and directly in ζ L do not influence the cross section in this special decoupling and exact alignment limit (s βα = 1).
Comparing the two plots at m H ± = 500 GeV shows that the cross section in Type-I is higher   for the low m H region, there will be in total 12 jets at the final state and also the decay products of the Z 0 boson (2 jets or 2 leptons). Apparently, this scenario produces many jets at the final state, and it is a challenge to reconstruct the W-boson, so the charged Higgses.

C. Decay channels of H ± in the favored region
In this scenario, the charged Higgs boson decays via two channels H + → tb and H + → τ ν τ for both Type-I and -II. In Fig. 11, it can be seen that the B(H + → tb) is the dominant decay channel for each Yukawa structure. Then, the same final state will be obtained with the non-alignment scenario Type-II. The subsequent decays of t → W b and W ± → qq(lν l ) forms the final state of the charged Higgs boson. Then, considering the Z → qq, the final state of the process is 4 jets + 4 b-tagged jets + 2 jets. Type-II Yukawa couplings are set.

D. Differential distributions
In this section, the differential cross sections are calculated as a function of the kinematical variables for each scenario, and comparison is performed. The computation is carried at √ s = 3 TeV, s βα or c βα is taken from the corresponding table given for the each scenario. Type-I (Type-II) in the favored region, respectively. The t β = 10 is set for all the scenarios.
In Fig. 12 (left), the differential cross section as a function of the rapidity of the Z 0 boson is presented. It can be seen that the Z 0 boson is produced more likely in the central region compared to the high rapidities for all the scenarios. If a rapidity cut of |y Z | ≤ 2 on the Z 0 boson is applied, then depending on the scenario ≈ 75 − 88% of the events could be captured. That shows it could be a useful kinematical cut to eliminate more events in the background. In Fig. 12 (center), distribution for the differential cross section as a function of the transverse momentum of the Z 0 boson (p Z T ) is plotted for each scenario. According to the p Z T distribution, the Z 0 boson is produced more likely with small transverse momentum. That will most certainly affect the kinematical properties of the decay products of the Z 0 boson. Applying a cut of p Z T < 400 GeV reaps more than ≈ 70% of the events. At last, the differential rate as a function of the rapidity difference of the charged Higgses is plotted in Fig. 12 (right), and it shows that the charged Higgses are produced very likely with a small rapidity difference, and ≈ 95% of the events fall in |∆y H − H + | ≤ 2. Moreover, the normalized rate is the same in all the scenarios.  (left): dσ/d|y Z | distribution, (center): dσ/dp Z T distribution, (right): dσ/d|∆y H − H + | distribution.

E. SM Backgrounds
The production of H + H − Z 0 in e + e − collider has a small cross section as it is presented in the previous section. Regarding the weakness of the signal, one may wonder if such a process could be extracted from the SM background. The decay chains of the charged Higgs in each scenario showed that the b-quark identification is vital in the reconstruction of the signal.
Besides, the number of jets and the b-tagged jets at the final state in real-world are not fixed due to the parton-branching of the quarks and the efficiency of the b-tagging as well as the mis-identification of them. By considering these points, there are many background channels which will shadow the process. Therefore, a large background is expected mainly from e + e − → qq, and where later quarks will hadronize to jets. Moreover, e + e − → qqZ, e + e − → tt, and e + e − → ttZ will contribute to the main background. On the other hand, the following processes could also contribute to the backround. These are e + e − → ttbbbb along with Z 0 boson or multiple numbers of light jets (u, d, c, or s type quark), e + e − → ttZH, e + e − → ttZZ, e + e − → ttHH, e + e − → ttbbZ, and e + e − → W + HW − HZ.
Hadronic decays of the W ± /Z 0 -boson along with the jets could mimic the final state of the process. It looks like the top quark is essential for reconstructing the charged Higgs as well as in the elimination of the various background channels. Besides, it is also possible to allow the leptonic decays of W ± /Z 0 bosons. However, since the branching ratios are small, the multiplication of B(W → lν l ) 2 · B(Z → ll) ≈ 4.4 · 10 −3 reduces the cross section by ∼ 1/225. Eventually, exploring the hadronic decays of W ± /Z 0 bosons gives more events in the detector. Therefore, jet finding algorithms will determine whether the process could be (P e + = +0.60). The option of upgrading the incoming electron and the positron beam to be polarized has the power to enhance the potential of the machine. If it is assumed that CLIC could produce a total of integrated luminosity of L ∼ 3 ab −1 at √ s = 3 TeV [30], then there could be more than ∼ 2.7 · 10 3 events assuming polarization.
The decay channels of the charged Higgs boson in each scenario are also investigated with the help of 2HDMC. The analysis shows that some channels come forward such as H + → tb, H + → W + h 0 , and H + → W + H 0 . The subsequent hadronic decays of the top quark, W ± /Z 0 bosons, and h 0 /H 0 bosons have higher branching ratio compared to the leptonic ones. In that case, the final state of the process contains 6 jets + 4 b-tagged jets in the non-alignment scenario with Type-I, and 6 jets + 8 b-tagged jets in Type-II. Moreover, the other two scenarios have the same final state though following different decays, and there are 6 jets + 4 b-tagged jets in the end. Unfortunately, since the number of jets in the final state is high, it will be hard to reconstruct the process, and that is the most significant disadvantage. High efficiency in b-tagging and reconstructing the W ± /Z 0 boson, then reconstructing the top quark is vital for the charged Higgs detection. Further, charged Higgs pair has higher production rate compared to the extra Z 0 boson in the final state.
Considering the background, a detailed Monte Carlo study is required to determine the significance and the acceptance of the signal in a detector. Some differential distributions for the charged Higgs and Z 0 boson are presented, but the best selection cuts with a higher elimination of the background signals require a full detector simulation and maybe a help of the artificial neural networks.
LHC experiment confirmed the existence of a neutral Higgs boson [38,[40][41][42], hereafter, the discovery of another scalar and even a charged one at the future colliders would be clear

VIII. ACKNOWLEDGEMENT
The numerical calculations reported in this paper were partially performed at TUBITAK ULAKBIM, High Performance and Grid Computing Center (TRUBA resources) and also at the computing resource of fencluster (Faculty of Science, Ege University). Ege University