Search for neutral Higgs bosons within type-I 2HDM at future linear colliders

In this study, we investigate observability of the neutral scalar (H) and pseudoscalar (A) Higgs bosons in the framework of the type-I 2HDM at SM-like scenario at a linear collider operating at s=\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}=$$\end{document} 500 and 1000 GeV. The signal process is e-e+→AH→ZHH→jjbb¯bb¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e^- e^+ \rightarrow A H \rightarrow ZHH\rightarrow jj b{\bar{b}}b{\bar{b}}$$\end{document} where jj is a di-jet resulting from the Z boson decay and bb¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b{\bar{b}}$$\end{document} is a b-jet pair from scalar Higgs boson decay. Several benchmark scenarios with different mass hypotheses are studied. The assumed signal process is mainly motivated by the large decay rate of A→ZH\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A\rightarrow ZH$$\end{document} (at high tanβ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tan \beta $$\end{document} and enough Higgs boson mass splitting) and H→bb¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H\rightarrow b{\bar{b}}$$\end{document} in the type-I. After event generation the detector response is simulated based on the SiD detector at the ILC. The ISR and beamstrahlung effects are also included for each center of mass energy assuming unpolarized beams. The simulated events are analyzed to obtain candidate mass distributions of the Higgs bosons. Among standard model backgrounds, the top quark pair production process is the main background but is under control. Results indicate that, in all considered scenarios, both Higgs bosons H and A are observable with signals exceeding 5σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$5\sigma $$\end{document} with possibility of mass measurement. Final results are shown as 5σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma $$\end{document} contours in (mH,mA\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_H,m_A$$\end{document}) or (mA,tanβ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m_A,\tan \beta $$\end{document}) space.


Introduction
The Standard Model (SM) of elementary particles has emerged as a significant achievement explaining plenty of phenomena. The existence of the Higgs boson, as one of the most striking predictions of the SM, attracted much attention even before its experimental verification [1,2]. The prediction of the Higgs boson was a direct consequence of the assumed scalar structure which was chosen to be the simplest possible one. To be specific, the SM assumes a single SU (2) doublet with four degrees of freedom as the Higgs a e-mail: majid.hashemi@cern.ch b e-mail: hosseinhaqiqat@gmail.com field leading to the prediction of a single Higgs boson [3][4][5][6][7][8]. However, there are some important motivations, namely the SM inability to explain the universe baryon asymmetry [9], supersymmetry [10], axion models [11], etc., which provide possibility of an extended scalar structure which leads to the prediction of existence of additional Higgs bosons. Assuming a scalar structure based on two SU (2) Higgs doublets, as one of the simplest possibilities, leads to the Two-Higgs-Doublet model (2HDM) [12][13][14][15][16][17][18][19] which provides interesting environment and phenomenological features. The two assumed Higgs doublets in the 2HDM carry eight degrees of freedom, three of which are absorbed by the three of the electroweak gauge bosons and the remaining five degrees of freedom finally lead to the prediction of five Higgs bosons. The lightest Higgs boson (h) is assumed to be the same as the discovered SM Higgs boson (the SM-like scenario) and the four additional Higgs bosons are thought of as yet undiscovered particles. These Higgs bosons include a neutral scalar (H ), a neutral pseudoscalar (A) and two charged (H ± ) Higgs bosons. The present paper concentrates on the H and A Higgs bosons and investigates observability of these particles in several benchmark points in the parameter space of the 2HDM at SM-like scenario.
The 2HDM, in its general formulation, predicts tree level flavor-changing neutral currents (FCNCs) which is not in agreement with experimental observations. However, imposing the discrete Z 2 symmetry, the tree level FCNCs are well avoided and four 2HDM types with different Higgsfermion coupling scenarios which naturally conserve flavor are obtained. Different types with different coupling scenarios provide different interesting characteristics [18]. The Type-X 2HDM was studied earlier to investigate the observability of the H Higgs boson with the help of large enhancements the Higgs-lepton coupling provides at high tan β values and led to promising results [20]. In this study, the Type-I 2HDM is chosen as the theoretical framework and the signal process through which the Higgs bosons observability is investigated is assumed to be e − e + → AH → Z H H → j jbbbb where j j is a di-jet resulting from the Z boson decay and bb is a b quark pair resulting from the scalar Higgs H decay. The pseudoscalar Higgs boson A undergoes the decay mode A → Z H which is dominant due to the SM-like assumption and also the non-zero mass splitting assumed between the A and H Higgs bosons. Scenarios with equal Higgs masses were also considered earlier leading to promising results [21]. The assumed signal process is also motivated by the large enhancement the scalar Higgs decay mode H → bb receives in low tan β regime. Such an enhancement is due to the H -fermion coupling in the Type-I which depends on cot β. Moreover, the A-Z -H vertex depends on sin(β − α) which is set to unity because of the SM-like assumption. Hence, the process e − e + → AH → Z H H is independent of tan β and thus, the scalar Higgs decay mode H → bb may receive large enhancements at low tan β values without any destructive effect on the Z H H production process.
Several benchmark points with different mass hypotheses are assumed and investigation of the Higgs bosons observability is done separately for each scenario. Masses of the Higgs bosons are chosen from intermediate and heavy mass regions and thus, the center-of-mass energy of 1000 GeV provides observation possibility in a relatively large portion of the space parameter. However, the center-of-mass energy of 500 GeV is also studied in this study since this energy is easily accessible to future linear colliders and it will be shown that a few scenarios can also be observed at this center-of-mass energy. The LHC can easily provide the center-of-mass energies considered here. However, since linear colliders provide a cleaner environment with less background processes and underlying events, a linear collider is assumed in this study.
Event generation is performed for each scenario independently and both beams are assumed to be unpolarised. The beamstrahlung effects [22] are taken into account and the detector response is simulated based on the SiD detector at the International Linear Collider (ILC) [23]. Simulated signal and relevant background events are analyzed to reconstruct Higgs bosons by first performing jet clustering and btagging and then applying proper selection cuts to enrich the signal events. Identifying proper bb and j jbb combinations and computing their invariant masses, candidate mass distributions of the Higgs bosons are obtained at the assumed integrated luminosities. The integrated luminosity is set to 500 f b −1 for all of the benchmark scenarios except for one scenario for which the integrated luminosity of 1000 f b −1 is assumed (for search for the A Higgs boson). It will be shown that, in all of the scenarios under consideration, both of the Higgs bosons H and A are observable with signals exceeding 5σ . Moreover, reconstructing masses of the Higgs bosons, it will be shown that masses of both of the Higgs bosons are measurable. In what follows, a brief introduction to the 2HDM and its different types is provided, and then analysis and results will be discussed.

Two-Higgs-doublet model
Employing two SU (2) Higgs doublets, the potential where 1 and 2 are SU (2) Higgs doublets, is assumed as the Higgs potential in a general 2HDM. Using one additional Higgs doublet in this model leads to the prediction of five Higgs bosons. The lightest Higgs boson (h) is assumed to be the same as the observed Standard Model Higgs boson and the four additional Higgs bosons are thought of as yet undiscovered Higgs bosons. Additional Higgs bosons include a neutral scalar (H ), a neutral pseudoscalar (A) and two charged (H ± ) Higgs bosons. To respect experimental observations, the discrete Z 2 symmetry is imposed to avoid tree level flavor-changing neutral currents [14][15][16]. Consequently, it is implied that the parameters λ 6 , λ 7 and m 2 12 must be zero. However, letting m 2 12 be non-zero, Z 2 symmetry is softly broken in this model. Assigning a value to tan β which is defined as the ratio of the vacuum expectation values of the two Higgs doublets, the parameters m 2 11 and m 2 22 are obtained by the minimization conditions for a minimum of the vacuum. Setting λ 6 and λ 7 to zero to respect the discrete Z 2 symmetry and working in the "physical basis", tan β, mixing angle α, m 2 12 and physical masses of the Higgs bosons must be determined to specify the model completely [12]. As a result of the imposed Z 2 symmetry, Higgs-fermion coupling scenarios are limited. Table 1 provides Higgs coupling to up-type and down-type quarks and leptons in the allowed types which naturally conserve flavor. The types "X" and "Y" are also called "lepton-specific" and "flipped" respectively. Assuming sin(β − α) = 1 [12], we choose the model to be SM-like in order for the h Higgs boson to be thought of as the observed SM Higgs boson. Consequently, h coupling to fermions reduce to corresponding couplings in the SM Yukawa Lagrangian. Consequently, the neutral Higgs part of the Yukawa Lagrangian takes the form [12,24] where ρ X factors are provided in Table 2.
As seen in Table 2, coupling factors of different types differ dramatically leading to considerable differences in phenomenological features of different types [18]. According to Table 2, factors of the type-I acquire the same values (cot β) resulting in an interesting environment in both low and high tan β regions. Searching for the Higgs bosons in the type-I in this study is highly motivated by the large enhancement the H → bb channel receives at low tan β values.

Signal process
The process chain e − e + → AH → Z H H → j jbbbb is assumed as the signal process in this study and the Type-I 2HDM is assumed as the theoretical framework to benefit from possible enhancements in low tan β regime. j j is a pair of jets resulting from the Z boson hadronic decay Z → j j. After the pseudoscalar Higgs boson undergoes the decay channel A → Z H, both of the scalar CP-even Higgs bosons experience the decay mode H → bb which receives a large enhancement due to the cot β factor in the Higgs-fermion coupling factors as shown in Table 2 and thus, is dominant in low tan β regime as long as the scalar Higgs mass m H is below the threshold of the on-shell top quark pair production. The initial e − e + collision is assumed to occur at the centerof-mass energies of 500 and 1000 GeV at a linear collider.

Benchmark points
Several benchmark points with different mass hypotheses are assumed in the parameter space of the 2HDM as shown in Table 3. Benchmark points corresponding to the two assumed center-of-mass energies are simulated and analyzed separately. The scalar Higgs boson mass m H is assumed to vary in range 150-250 GeV and the mass splitting between the H and A Higgs bosons is assumed to range from 50 to 100 GeV. tan β is set to 10 in all scenarios resulting in a considerable enhancement in the assumed scalar Higgs boson decay channel.

Theoretical constraints
To ensure that the considered scenarios are consistent with theoretical constraints, potential stability [25], perturbativity and unitarity [26][27][28][29] of each scenario is checked with the use of 2HDMC 1.7.0 [30,31] and the allowed range for m 2 12 parameter is provided in Table 3.

Experimental constraints
The assumed masses of the Higgs bosons are consistent with results of 86 analyses as checked by HiggsBounds 4.3.1 [32] and HiggsSignals 1.3.0 [33]. The experimental constraint [34,35], based on the measurement performed at LEP [36], limits the deviation of the parameter ρ = m 2 W (m Z cos θ W ) −2 from its SM value. To respect this constraint, the Higgs bosons A and H ± are assumed to have the same masses in all of the benchmark points. This is because of the demonstration provided in [37,38] that concludes the deviation of the ρ parameter from its SM value is negligible if any of the conditions is met. As reported in [39], the ATLAS collaboration has excluded mass regions m A = 310−410, 335−400, 350-400 GeV for m H = 150, 200, 250 GeV respectively at tan β = 10 in the type-I. Both CMS and ATLAS also put the limit m A > 350 on the pseudoscalar Higgs mass for tan β < 5 in the type-I [40,41]. Another study by ATLAS collaboration excludes the mass range m H = 170 − 360 GeV for tan β < 1.5 in this type [42]. It is obvious from Table 3 that the assumed scenarios are consistent with the current constraints resulting from direct LHC observations. In Fig. 1 the selected benchmark points are shown in the (m H ,m A ) plane together with the current LHC exclusion contour at tan β = 10 [39]. What is important here is that the selected points are chosen in a way to cover regions where the mass difference between the scalar and pseudo-scalar Higgs bosons is below the Z boson  The selected benchmark points and the current LHC excluded area [39] mass. Starting from any of the selected BPs, BR(A → Z H) increases by moving upward to higher values of m A and will be close to unity for m A − m H ≥ m Z . Therefore the only reason for the signal significance to decrease is the phase space saturation when heavy Higgs bosons are considered.

Flavour physics constraints
In the context of the type-II, flavor physics data [43,44] puts the constraints m H ± > 570 GeV (tan β > 2) and m H ± > 700 GeV (tan β < 2) on the charged Higgs mass. However, as indicated in [44], no condition limits the Type-I for tan β > 2. Moreover, it is shown by a review of LHC, LEP and Tevatron results that there is no exclusion around sin(β − α) = 1 for m H/A/H ± = 500 GeV in the Type-I 2HDM [45].

Constraints from precision measurements of the couplings
The SM Higgs boson couplings with fermions and gauge bosons have been measured within reasonable uncertainties at LHC. However, future experiments such as HL-LHC [46][47][48] and ILC [49,50] are expected to achieve more precision on those observables. In models with extended Higgs sectors, the couplings of SM Higgs with other particles may in general deviate from the corresponding SM predictions due to loop corrections which involve extra Higgs bosons.
Although the current results are consistent with SM predictions within the uncertainties, a deviation may be found at future colliders where higher precision measurements are expected. Therefore calculating radiative corrections to the couplings of the SM Higgs boson with SM particles can fingerprint extended Higgs sectors [51,52]. Such calculations can then be used to set limits on the masses of the extra Higgs bosons of the new model [53]. In order to achieve this goal, coupling constants have been calculated at loop level including radiative corrections. Assuming a 2HDM type-II as the theoretical framework and choosing different scenarios for the coupling constant uncertainties, upper limits have been obtained on the masses of the extra Higgs bosons as a function of tan β [53].
Results of the analysis reported in [53] can be translated to a 2HDM type-I, however, since in type-I, the Higgs couplings with fermions are proportional to cot β, no serious constraint is expected at high tan β values which are assumed in this work.

Decay channels
The assumed decay channel A → Z H in the signal process receives an enhancement due to the assumption sin(β −α) = 1 which is needed for the SM-like scenario, since the A-Z -H vertex depends on sin(β − α) in the Type-I 2HDM. Moreover, the assumed non-zero mass splitting between A and H Higgs bosons in the range 50-100 GeV facilitates the possibilities for this decay channel. On average, we obtain BR(A → Z H) 0.72 for the assumed benchmark points of Table 3 using 2HDMC 1.7.0. The A-Z -H vertex appears twice in the signal process. First in the production process e + e − → Z * → H A and then in the assumed decay mode for the pseudoscalar Higgs boson. Therefore, due to the assumption sin(β−α) = 1, the signal process is independent of tan β as long as the scalar Higgs decay mode H → bb is not considered. This fact provides an opportunity for the signal to benefit from the enhancement received by the decay channel H → bb at low tan β values, since lowering tan β value doesn't result in any decrease in the cross section of the process e + e − → Z * → AH → Z H H. Setting the value of 10 for tan β for all of the scenarios, on average, we obtain BR(H → bb) 0.61 for the assumed benchmark points using 2HDMC 1.7.0.
The b quark pairs resulting from the H decays annihilate into hadronic jets and are identified by first performing jet reconstruction and then performing a proper b-tagging algorithm. The identified b-jets are then used to reconstruct the H Higgs bosons. Although choosing the hadronic decay channel Z → j j for the Z boson may give rise to more errors in the final results due to the uncertainties arising from jet reconstruction and b-tagging algorithms, the hadronic decay channel is chosen since the large branching ratio of this channel (BR Z → j j 0.69) may fully compensate for the potential arising errors. Reconstructing the Z boson by the identified jets, the A Higgs boson is also reconstructed with the help of the reconstructed H Higgs boson.
Background processes relevant to the considered signal process include W ± pair production, Z /γ production, Z pair production and top quark pair production. Cross sections of the signal and relevant background processes are computed at the center-of-mass energies of 500 and 1000 GeV by CompHEP 4.5.2 [54] and are provided in Tables 3  and 4. According to Table 3, observing scenarios with heavier Higgs masses must be more challenging since the Higgs masses and cross section behave oppositely to each other.

Event generation
In order to take the beamstrahlung effects into account and simulate the detector response, event generation is performed in several steps for each benchmark scenario. Hard scat- tering part of the signal and background events (partonlevel part which doesn't include parton showers, hadronization, etc.) are generated using CompHEP 4.5.2. Simulation of the beamstrahlung effects is also performed by CompHEP with the use of the beam parameters provided in Table 5 with the assumption that both beams are unpolarised. To generate the remaining part of the events including multi-particle interactions, parton showers and hadronization, the SLHA (SUSY Les Houches Accord) files generated by 2HDMC 1.7.0 as well as the parton-level events generated by CompHEP are passed to PYTHIA 8.2.15 [56]. The SLHA files contain basic parameters of the Type-I 2HDM including coupling constants, branching ratios, etc.
Using the input files, PYTHIA performs further processing to complete the events. The generated events are then internally used by DELPHES 3.4 [57] to simulate detector response with the use of the DSiD detector card which is based on the full simulation performance of the SiD detector at the ILC. The anti-k t algorithm [58] from FASTJET 3.1.0 package [59,60] is used to perform jet reconstruction with the cone size R = ( η) 2 + ( φ) 2 = 0.4, where η = −ln tan(θ/2) and φ (θ ) is the azimuthal (polar) angle with respect to the beam axis. DELPHES output data is stored as ROOT files [61] and contains reconstructed jets and also b-tagging flags by which b-jets are identified.

Event selection and analysis
Analysis begins by counting the number of reconstructed jets satisfying the kinematic thresholds where p T is the transverse momentum. Based on the jet multiplicity distributions of Fig. 2 which are obtained for the two assumed center-of-mass energies, the selection cut where N jet is the number of jets, is applied. Using b-tagging flags, b-jets are identified and the b-jet multiplicity distributions of Fig. 3 are obtained for the b-jets satisfying the lower threshold p T b-jet ≥ 20 GeV . The condition where N b-jet is the number of b-jets, is then applied to rule out events with less than three b-jets. b-jets present in each events are analyzed to find the true b-jet pairs which originate from the H Higgs bosons decays. In events with three b-jets, the value of R bb , where R follows the definition b-jet multiplicity is satisfied. The b-jets are expected to make a more collinear pair at 1 TeV due to the higher energy of decaying particle (the Higgs boson). Therefore R bb tends to smaller values at 1 TeV. In events with at least four b-jets, all possible combinations of two b-jet pairs are considered and the difference between the invariant masses of the two b-jet pairs is computed for each combination, and the b-jet pairs of the combination with minimum invariant mass difference are identified as true pairs if the condition 8 is satisfied. Identifying true pairs, the selection cut where R follows the definition of Eq. (7), is satisfied. Again the condition is tighter at 1 TeV due to the more collinearity of the Z and H bosons at 1 TeV. In events with two reconstructed H Higgs bosons, R Z H is computed for each one of the two possible Z H combinations and the combination with smaller R Z H value is identified as the true combination. Having true combinations identified, the selection cut where N Z H is the number of identified true Z H combinations, is applied and the A Higgs boson is reconstructed using the identified Z H combination. Applying the selection cuts, event selection efficiencies are obtained for different signal and background processes at √ s = 500 and 1000 GeV and the results are provided in Table 6. According to the analysis, the Higgs bosons H and A are reconstructed after three and four selection cuts respectively. Thus, total efficiencies corresponding to the first three selection cuts and all of the four selection cuts are also provided.
Computing the invariant mass of the identified bb and j jbb combinations for events surviving the selection cuts, Higgs candidate mass distributions of Figs. 4 and 5 are obtained for √ s = 500 and 1000 GeV respectively. Distributions corresponding to different scenarios are shown separately for each one of the Higgs bosons. Signal and different background contributions are also shown separately. As seen, the tt process has the most contribution among the background processes, and is, however, well under control. In all of the distributions, signal contribution is seen as a significant excess of data on top of the total background. The distributions are normalized based on L × σ × , where L is the integrated luminosity, σ is the cross section and is the selection efficiency. The integrated luminosity is set to 500 f b −1 for all of the distributions except for the distribution of Fig. 5k which uses the integrated luminosity of 1000 f b −1 . Signal cross sections are obtained by multiplying the total cross sections provided in Table 3 by corresponding branching ratios of the decay modes A → Z H, Z → j j and H → bb. Background cross sections are taken from Table 4. Selection efficiencies for A mass distributions are taken from Table 6 and efficiencies corresponding to H mass distributions are obtained by computing the total number of signal reconstructed H Higgs bosons divided by twice the number of simulated signal events. Proper fit functions are fitted to the total background (B) and signal+total background (S+B) distributions of Figs. 4 and 5 and results are shown with associated error bars. Fitted curves show significant peaks near the generated Higgs masses. Fitting is performed by ROOT 5.34 [62]. The fit function used for B distributions is a polynomial function and the fit function used for S+B distributions is the combination of a polynomial and a Gaussian function. The Gaussian function covers the signal and the polynomial covers the total background. First, the polynomial function is fitted to the total background, and then the parameters of the fitted polynomial are used as input for S+B fit. The "Mean" parameter is one of the Gaussian fit function parameters and provides the center of the signal peak. The values obtained for the "Mean" parameter are considered as the reconstructed masses (m Rec. ) of the Higgs bosons and are provided in Table 7. The generated masses (m Gen. ) are also shown for comparison. According to Table 7, a difference is seen between the generated and reconstructed masses which can be explained by the uncertainties arising from jet clustering algorithm and jet mis-identification, jet mis-tag rate, fitting method and choice of the fit function, errors in energy and momentum of the particles, etc. Optimization of the jet clustering algorithm, b-tagging algorithm and fitting method may reduce the errors. However, since such optimizations lie beyond the scope of this paper, a simple off-set correction is applied to reduce the errors in this study. The off-set correction is applied as follows. According to Table 7a, on average, the reconstructed masses of the H and A Higgs bosons are 22.85 and 29.55 GeV smaller than the corresponding generated masses respectively. Calculating the same values for Table 7b, we obtain 17.85 and 36.17 GeV for the Higgs bosons H and A respectively. To reduce the errors, reconstructed H and A Higgs masses are increased by the same values and the results are provided in Table 7 as corrected

Signal significance
To assess the observability of the Higgs bosons in the considered scenarios, signal significance is computed for each of the candidate mass distributions of Figs. 4 and 5 by counting the number of signal and background candidate masses in the whole mass range. Computation of the signal significance is based on the integrated luminosity of 500 f b −1 for all of the distributions except for the distribution of Fig. 5k which uses the integrated luminosity of 1000 f b −1 . Computation results, namely total signal selection efficiency, number of signal (S) and background (B) candidate masses, signal to background ratio and signal significance at the center-ofmass energies of 500 and 1000 GeV are provided in Table 8.

Discovery contours
The results described in the previous section are limited to the selected benchmark points. In this section an extrapolation is made to a wider region in the parameter space. In order to do so, starting from each value of m H (in smaller increments than the chosen benchmark points with mass increments of 10 GeV), the m A is scanned from m H + 50 GeV up to a point where the signal significance decreases to 5σ as the minimum. This approach is first based on tan β = 10 but can be done for other values. Figure 6a shows the color plot for tan β = 10 keeping information of the significance at each point. The border of the blue region is obviously the limit of the 5σ . In Fig. 6b, the contour in m A vs tan β plane has been obtained with m H being 50 GeV lighter than m A at every point. As seen, higher tan β values result in no sizable change in the signal significances and the signal is only sensitive to low tan β values. Therefore the plot of Fig. 6a can be viewed as a contour for tan β > 10.

The signal sensitivity at CLIC 1.5 and 3 TeV
Since the signal cross section and the signal to background ratio are high enough to provide signal exceeding 5σ at √ s = 0.5 and 1 TeV, it is interesting to investigate the signal sensitivity of a higher energy collider like CLIC at 1.5 or 3 TeV to the signals studied in this analysis or even heavier Higgs bosons which may be out of the reach of ILC. To this aim the total cross section of the signal (each BP separately) and background processes were extrapolated to the CLIC region at 3 TeV. As Fig. 7a, b show, the signal and background cross sections decrease but the signal to background ratio remains almost the same as in 1 TeV (∼ 10 −3 ). Assuming the signal significance measure as S/ √ B, the new signal significance of each selected BP at 3 TeV would be S/B times √ B and since S/B shows no sizable change, the signal significance decreases like √ B. Therefore a higher energy collider may not improve these results. Other parameter space points involving heavier Higgs bosons will also have smaller cross sections than the "light" points selected here. Therefore their situation is not better than the selected BPs in the current analysis either. As a conclusion, it is expected that a 1 TeV ILC works better for the light and moderate Higgs masses unless the m H + m A reaches the kinematic threshold of the collider (i.e., √ s) in which case using a higher energy collider like CLIC is inevitable.

Conclusions
The signal process e − e + → AH → Z H H → j jbbbb based on two Higgs doublet model type I was staudied Table 8 Total signal selection efficiency ( T otal ), number of signal (S) and background (B) candidate masses, signal to background ratio, signal significance and the assumed integrated luminosity in different scenarios at the center-of-mass energies of (a) 500 and (b) 1000 GeV √ s = 500 GeV √ s = 1000 GeV    at a future linear collider like ILC at 0.5 and 1 TeV. The integrated luminosity was set to 500 f b −1 for all selected benchmark points. Detector simulation was performed by DELPHES using SiD detector card. Results indicate that all selected points are in the reach of ILC with the possibility of mass reconstruction for both scalar and pseudo-scalar Higgs bosons at the same time with a single analysis. The coverage in tan β direction reaches tan β = 60 with almost the same signal sensitivity for tan β higher than 10. The signal significance in m H , m A plane was also obtained by extending the analysis to higher values of m A . A comparison with the current LHC results shows that a very larger area can be covered by the current analysis with LHC contour embedded in a part of the whole 5 σ contour. Extrapolation to the CLIC center of mass energy of 1 and 3 TeV shows no gain in the signal significance with increasing the center of mass energy to higher values than 1 TeV unless forced by the kinematic requirements. Therefore the analysis is more relevant to ILC at 1 TeV.