Search for a lighter Higgs boson in Two Higgs Doublet Models

We consider present constraints on Two Higgs Doublet Models, both from the LHC at Run 1 and from other sources in order to explore the possibility of constraining a neutral scalar or pseudo-scalar particle lighter than the 125 GeV Higgs boson. Such a lighter particle is not yet completely excluded by present data. We show with a simplified analysis that some new constraints could be obtained at the LHC if such a search is performed by the experimental collaborations, which we therefore encourage to continue carrying out light diphoton resonance searches at s=13\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s}=13 $$\end{document} TeV in the context of Two Higgs Doublet Models.


Introduction
After the discovery of a Higgs boson at the LHC in 2012 [1,2], many studies, both from the theoretical and experimental side, have considered extensions of the Standard Model (SM) with an enlarged scalar sector. Concerning this scalar sector of physics beyond the Standard Model (BSM), most studies have considered the possibility of new scalars heavier than the 125 GeV Higgs boson which was discovered at the LHC. It is however possible to have a spectrum in which lighter scalars are present together with an SM-like Higgs boson at 125 GeV. Among these possibilities there are detailed BSM models as well as effective descriptions including only the extended scalar sector. Two Higgs Doublet Models (2HDMs) constitute one of the simplest possibilities, where the SM Lagrangian is extended by the addition of a second scalar doublet. Previous phenomenological studies describing the possibility of lighter Higgs bosons include [3][4][5][6][7][8][9], while for a recent study in supersymmetry (which naturally includes two doublets) we refer the reader to [10]. At masses below 125 GeV, the main search channel at the LHC is the di-photon decay channel [11,12]. This paper is organised as follows: in section 2 we describe the theoretical set-up adopted in our analysis; section 3 is dedicated to the study of present constraints

JHEP12(2016)068
Type I Type II Flipped Lepton Specific (Type Y) (Type X) Up-type quark Table 1. The different possible couplings between the SM fermions and the two scalar doublets in 2HDMs.
coming from flavour physics, electroweak precision tests, theoretical bounds, direct LEP constraints on the scalar sector and LHC limits given by a 125 GeV Higgs boson; section 4 contains the cross-section and branch ratio calculations for a light scalar Higgs boson, a study of the parameter space of the different 2HDMs and a comparison with the CMS low mass di-photon analysis at 8 TeV [11]; section 5 is dedicated to the study of the case where the lighter resonance is pseudo-scalar; finally we present our conclusions in section 6.

Two Higgs Doublet Models
We here briefly describe the theoretical framework of 2HDMs, see [13] for a general discussion. The 2HDMs are a simple extension of the Standard Model including two complex SU(2) doublets, φ 1 and φ 2 . In order to avoid flavour-changing neutral currents, one can introduce a Z 2 symmetry so that all fermions of a given electric charge couple to at most one Higgs doublet. These couplings can occur in different ways; the convention usually adopted is given in table 1. The most generic 2HDMs potential constrained by the Z 2 symmetry can be written as: where all the parameters are real. The parameter m 2 12 is responsible for a soft breaking of the Z 2 symmetry. The two scalar doublets acquire vacuum expectation values (vevs): with v ≡ v 2 1 + v 2 2 . After symmetry breaking we are left with five physical scalars: two neutral CP-even states h and H, one neutral CP-odd state A and two charged ones H ± . In order to move from the potential of eq. (2.1) to mass-eigenstates, one needs to introduce two angles: β, defined as tan β = v 2 v 1 , which rotates the two doublets in a basis where only one of them acquires a vev, and α which mixes the CP-even scalar states to give mass-eigenstates. The parameters of the potential can thus be translated into an equivalent set of parameters in JHEP12(2016)068 the physical basis: where v is set to the electroweak scale, and one of the masses of the CP-even states should be equal to the measured Higgs boson mass. The masses of the two CP-even states are ordered with m h < m H , where we will call h the light Higgs boson and H the heavy Higgs boson of the model. The couplings between the neutral Higgs bosons and the fermions and gauge bosons are summarised in table 2. In the rest of the study we will use the input parameters of the physical basis: we fix v = 246 GeV and the heavy Higgs boson H of the model is identified with the Higgs boson discovered at LHC, m H = 125 GeV, while the remaining six parameters are left free.

Bounds on 2HDMs
As we briefly discussed in the previous section, one of the simplest modifications of the SM consists in incorporating two scalar doublets, imposing custodial symmetry in order to allow satisfying the electroweak precision tests.  Table 3. Experimental values of the oblique parameters with 1σ uncertainty and correlations between them [16].
reported results on the search for a light resonance in di-photon final states [11], giving the observed upper limit at 95% confidence level (C.L.) on the cross-section times branching ratio as a function of the mass of a light Higgs boson between 80 GeV and 110 GeV. In the following we list the different constraints we use to impose bounds on the model. We split them in three classes: indirect constraints, LEP constraints and LHC constraints.

Indirect constraints
The indirect constraints we apply on the 2HDMs parameter space include limits on the oblique parameters S, T and U [14] due to electroweak precision tests, flavour constraints and theoretical requirements due to ensure stability of the potential, unitarity and perturbativity.
The oblique parameters are computed in the model via the program 2HDMC [15] and compared to the experimental limits [16] at 2σ (see table 3 for a recap of the updated experimental values with 1σ uncertainties and the correlations between them).
The stability of the potential is needed in order to allow symmetry breaking with a stable vacuum, thus the potential of the theory needs to be bounded from below. This condition requires [13]: In addition we require to have tree-level perturbative unitarity for the scattering of Higgs bosons and the longitudinal parts of electroweak gauge bosons [17]. In order to trust perturbative calculations, we add a condition on the quartic Higgs bosons couplings The three conditions detailed above are also computed via the 2HDMC program.
Once the previous requirements are satisfied, the available parameter space is tested against flavour bounds. We look at the branching ratios BR(B → X s γ) and BR(B s → µ + µ − ), which obtain new contributions from the charged Higgs bosons and the neutral ones respectively and at the isospin asymmetry ∆ 0 (B → K * γ) and the ∆M d frequency oscillation which are sensitive to the presence of charged Higgs bosons. The value of each process is computed in the 2HDMs via the program SuperIso [18,19] and then compared to the experimental limits at 2σ. In order to take into account the theoretical uncertainties in the 2HDMs calculation, which are not evaluated in SuperIso, we add to the experimental  [20] 0.543 ± 0.091 ps −1 [25] 0.091 ps −1 Table 4. Values of the experimental and theoretical flavour constraints.
1σ uncertainty σ Exp of each process the 1σ theoretical uncertainty σ Th of this same process computed in the SM given by the most recent theoretical calculations. The combined error σ comb can then be obtained via: A summary of the results we use is available in table 4.

Direct LEP constraints
The HiggsBounds program [26][27][28][29] is a tool able to test a model against experimental data coming from LEP, Tevatron and the LHC. The program can be interfaced with 2HDMC which will give appropriate inputs to HiggsBounds. In our analysis we use HiggsBounds version 4.2.1 with the LEP experiment constraints only, in order to impose LHC constraints separately. 2HDMC gives HiggsBounds a parton-level input for the three scalar Higgs bosons and the two charged ones. The exclusion test at 2σ is then performed on the five physical scalars of the theory. HiggsBounds returns a binary result indicating if the specific model point has been excluded at 95% C.L. or not.

LHC Higgs boson constraints
What we call "LHC constraints" are restrictions coming from experimental results on the 125 GeV Higgs boson, i.e. the 2HDM heavy Higgs boson H, in our case. To implement such limits, we use the exclusion contours in the plane of the signal strength for each individual production mode µ VBF/VH vs. µ ggh/tth given by the combined ATLAS and CMS experiments at Run 1 [30]. Assuming a Gaussian profile for the likelihood L at 68% C.L., each exclusion contour for a specific decay channel Y obeys the following equation: where µ ggH/ttH,Y and µ V BF/V H,Y are the data best fit values and a Y , b Y and c Y are the parameters of the ellipse. These five parameters fully describe the ellipse. We fit the ellipses for each decay channel Y = {W W, ZZ, γγ, τ τ, bb} and hence obtained the parametrisation for each of them (see [31] for more details). We then compute the χ 2 Y value in the 2HDM JHEP12(2016)068 using equation (3.3). For this, we assume the following relations: table 2).
Combining the log-likelihood ratios, we obtain: with p j the set of free parameters on which the function depends and p j their value minimising the χ 2 function. According to Wilks's theorem, the ∆χ 2 function follows a χ 2 distribution with a number of degrees of freedom equal to the number of free parameters. In our case, we have six degrees of freedom (κ 2 H→bb , as BR 2HDM H→ZZ is linked to BR 2HDM H→W W ). This choice in the free parameters implies that we assume there is no correlation between the kappas and the branching ratios, which is correct as long as the deviation of the branching ratio is not too large with respect to the Standard Model values. A point in the 2HDM parameter space passing the LHC constraints, therefore, has a ∆χ 2 value lower than 12.85, which is the value at 95% C.L. for a 6 degrees-of-freedom χ 2 distribution.

Search for a lighter scalar Higgs boson in the 2HDMs
A light resonance decaying into two photons is being searched for by CMS [11] in the range of mass between 80 and 110 GeV. In this section we will explore the possibility that the signal may be given by the light scalar state in the 2HDMs. To compare with the experimental sensitivity (in particular at 8 TeV), we need to compute the expected production cross-sections in the different production modes and branching ratios into the observed final states. In the following subsections we illustrate the procedure we followed and the results used in the present work to obtain restrictions on the parameter space of the various 2HDMs. In 4.1 we discuss the calculation method used for cross-sections times branching ratios. Then we apply, in section 4.2, the present bounds coming from the three sets of constraints defined in the previous section in order to define the available parameter space. We finally test the sensitivity of the CMS low mass di-photon analysis at the LHC Run 1, in section 4.3, in the available parameter space for the four types of 2HDMs. To do so, we rely on a scan on the six free parameters in the physical basis.

Cross-sections and branching ratios
We use the program 2HDMC [15] version 1.7.0 to compute the branching ratios of the different Higgs bosons of the theory. As input, the program requires a numerical value for each of the JHEP12(2016)068 seven parameters of the physical basis and provides, as output, the total width, branching ratios and couplings at next to leading order (NLO) for each Higgs boson.
The cross-sections can also be computed via programs like SusHi [32]: however, the output is restricted to the gluon fusion and bb production modes while SusHi does not provide vector boson fusion production (VBF) nor associated production with gauge bosons (VH). In order to overcome this restriction, and to quicken the calculation, we compute the cross-sections using an approximation that we have briefly introduced in section 3.3 and that we denote in the following as the "kappa trick". Defining the generic parameter for a specific decay channel Y , we approximate the cross-sections as: The second equation has such a simple form because, as the couplings of the light scalar Higgs boson to the W and Z bosons are rescaled in the same way compared to the SM couplings (cf table 2 The SM cross-section is taken from the LHC Higgs Cross-Section Working Group [33]. The kappas are computed thanks to the output given by 2HDMC. Hence we are able to compute the cross-section times branching ratio of the two neutral scalar Higgs bosons using only the 2HMDC program, via equation (4.1).
It is pertinent at this point to comment on the level of validity of this approximation. The cross-section production in VBF and VH mode should not cause any problem as the leading effect arises at tree level, however for the gluon fusion mode a loop induced coupling is present and thus it is important to check the validity of the "kappa trick" for this production mode. Note indeed that for loop induced vertices the use of an effective kappa factor is not always appropriate and more general parameterisations exist (see for example [34,35]). In order to explore this issue and establish if this simple approximation could be used, we performe a comparison between the cross-sections in gluon fusion obtained via the program SusHi and the ones obtained with the "kappa trick". As 2HDMC only considers NLO corrections, we also ran SusHi at NLO. The SM inputs required by the two programs are set to the recommended values given by the Particle Data Group [24] summarised in table 5. The 2HDM inputs used are given in table 6: we chose to fix all the parameters except the mass of the light neutral scalars, whose cross-section we want to test. We use SusHi version 1.6.0 together with LHAPDF 6.1.6 [36]. The parton distribution functions used in the program are MMHT201468cl for LO and PDF4LHC15 mc for NLO and NNLO [37]. The renormalization and factorization scales µ R and µ F for the gluon fusion process are set to µ R = µ F = m φ /2 with φ = {h, H, A} [38]. The bb production mode proposed by SusHi is turned off.
Fixing six of the seven free parameters, we allow only m h to vary between 80 to 110 GeV with a step of 1 GeV between each point (see table 6). The results are plotted in figure 1. The dashed blue line corresponds to the cross-section computed with the "kappa trick", the dotted red line to the computation with SusHi and the solid green line to the deviation between the two, computed as:   Table 5. SM input parameters [24].
[80;110], step 1. 125 550 600 5 -0.2 30 Table 6. Input parameters in the 2HDM Type I for the comparison between SusHi and "kappa trick" cross-sections. The plot shows a deviation of less than 3% for the whole mass range and this deviation is stable upon modification of the values of the input parameters (see figure 16 in the appendix). As it stays within the range allowed by the uncertainties (theoretical, PDF and α s ) calculated by the LHC Higgs Cross-Section Working Group [33], we consider this test as having validated our method for the light Higgs boson. For completeness, we make a similar analysis for the heavy Higgs boson at 125 GeV (see figure 17 in the appendix) finding deviations less than 1% at m H = 125 GeV. In the rest of the study, therefore, we will use the "kappa trick" approximation to compute the cross-section of the light and heavy scalar Higgs bosons in the 2HDMs.

Constraining the 2HDMs parameter space
In this section we study the influence of the three sets of constraints defined in section 3 (indirect, LEP and LHC constraints) on the free parameters. For this purpose we generate a set of one million points for each of the four different types of model defined in table 1 with random values for each of the free parameters. The available ranges we use in the simulation are given in table 7. The range of variation for m h corresponds to the mass range available in the CMS di-photon analysis. The lower bound of 80 GeV for m H ± comes from the bound obtained at the LEP experiment [39]. The ranges for m A and m 2 12 , although not totally general, are the result of previous quick scans that we will not show in this paper and which eliminate areas with a very low density of points passing the three sets of constraints (indirect, LEP and LHC constraints).
Once the points are generated, we impose the three kinds of constraints detailed above: indirect ones, direct LEP and LHC Higgs boson ones.
In figure 2, all the generated points are plotted in the plane m A vs. m H ± . The upper left panel corresponds to Type I, the upper right to Type II, the lower left to the Flipped JHEP12(2016)068 model and the lower right to Lepton Specific model. The points passing only the indirect constraints are plotted in green, those passing indirect and LEP constraints are in blue and those passing indirect, LEP and LHC constraints are in red. We will use these same conventions in the rest of this section.
Firstly we can see that the m A and m H ± masses are very correlated: when m A and m H ± grow, the indirect constraints force them to be near the black line corresponding to m A = m H ± . This is due to the T parameter which is very sensitive to these two masses and enforces them to be close to each other. Looking only at the red points, those which pass the three sets of constraints we defined previously, we can see that the two masses are bounded. In Type I, we find that most of the red points lie in the ranges m A ∈ [60 GeV; 650 GeV] and m H ± ∈ [80 GeV; 630 GeV]. In Type II and Flipped, the two masses are much more constrained m A ∈ [400 GeV; 650 GeV] and m H ± ∈ [430 GeV; 630 GeV]: this is due to the fact that the down-type quarks couple now to the φ 1 doublet instead of the φ 2 doublet as in Type I, thus the BR(B → X s γ) flavour limit imposes a very strong constraint on the mass of the charged Higgs bosons (see figure 18 in the appendix). Associated with the T parameter constraint, it imposes also the bounds on the pseudo-scalar mass. The Lepton Specific case is very similar to Type I as the couplings of the down-type quark are the same. We find m A ∈ [80 GeV; 630 GeV] and m H ± ∈ [80 GeV; 630 GeV] to be the preferred regions. We should remark that these bounds are not absolute and that there may be red points exceeding these bounds. However, our simulation shows that the bulk of the allowed points are inside the ranges, so that we decided to use them in order to increase the statistics of our scan.
Looking now at the plane tan β vs. sin(β − α) (shown in figure 3) we can constrain in the same way the tan β parameter. If it is difficult to impose an upper limit in all types as we lack statistics for high values of tan β and we see a few red points up to the upper value, nevertheless we can impose a lower bound of tan β > 1.2 for the four different types.
The bounds on sin(β − α) can be more easily seen in the plane m h vs. sin(β − α) (shown in figure 4): we see that m h is not constrained as red points span the whole range of masses. For sin(β − α), the allowed range is close to zero, which is consistent with our choice of m H = 125 GeV: as sin(β − α) 0, we have cos(β − α) 1 which means that the couplings of the heavy Higgs boson H to the gauge bosons are close to the SM ones. We are therefore close to the alignment limit [8]. We find that the preferred ranges are sin Finally, looking at the plane m 12 vs. sin(β −α) we can constrain the last free parameter (see figure 5). We cannot put any lower bound on m 2 12 but we find m 2 12 < (100 GeV) 2 in the four different types.
The previous results show that the range of the free parameters can be further limited in order to increase the statistics of the allowed points. In addition to this, as we are interested in checking the sensitivity to a lighter Higgs boson at LHC Run 1 in the diphoton decay channel, we can further restrict the areas of interest to where the red points correspond to relatively high values of cross-section times branching ratio to two photons. The minimum value of the CMS observed upper limit [11] is 0.032 pb in the gluon fusion JHEP12(2016)068 channel, obtained for m h = 103 GeV and 0.019 pb in the VBF/VH channel, obtained for m h = 100.5 GeV. Keeping these values in mind, we can look at the predicted 2HDM crosssection times branching ratio values as a function of sin(β − α). We plot the results for the gluon fusion production mode in figure 6 and for VBF/VH production mode in figure 7. The red dotted line corresponds to the minimum value of the CMS observed upper limit for each of the production modes. If all the red points are below this line, it means that CMS was not sensitive to a lighter Higgs boson in this particular channel at LHC Run 1.
The first important result we can extract from these figures is that in the Type II, Flipped and Lepton Specific models, CMS had no sensitivity to a lighter Higgs boson at LHC Run 1 in the h → γγ decay channel, neither in the gluon fusion nor in the VBF/VH production mode. Therefore we will not carry on with these types any further. Looking at the results for Type I, we can see that there is no sensitivity in the gluon fusion channel. However, in the VBF/VH channel, we find red points above the dashed line. As the value of the CMS observed upper limit depends on the mass of the light Higgs boson considered, the dashed line represented on the plots is not an absolute bound. Some of the red points above it can be de facto below the CMS observed limit, but it is a good indication of the potential capability of the channel for some exclusion. We can therefore expect to have some sensitivity in the VBF/VH channel.
We can exploit        are close to the CMS analysis limit sensitivity. We choose a lower bound at 0.01 pb to select the points, which corresponds to sin We can similarly work with tighter ranges for the parameters tan β and m 2 12 (see figure 8). We choose tan β ∈ [2; 12] and m 2 12 ∈ [−(100 GeV) 2 ; +(100 GeV) 2 ]. After having defined the allowed parameter region, and the more promising region with respect to the di-photon search, we are now ready to perform a second "focused" simulation and make a detailed comparison with the sensitivity of the CMS search at 8 TeV.

Comparison with the CMS low mass di-photon analysis
We thus perform a new scan with one million points, this time for Type I only, using the restricted parameter ranges we found in the previous section (see table 8). We remind the reader that for m A and m H ± the new range results only from the three sets of constraints (the indirect, LEP and LHC constraints) we imposed. For the parameters sin(β − α), tan β and m 2 12 it results from our choice to restrict the scan to areas with large value of σ VBF/VH × BR h→γγ (above 0.01 pb), as explained in section 4.2.
The resulting points of this second scan are plotted in figure 9 in the plane σ ×BR h→γγ in the gluon fusion production mode (left panel) and the VBF/VH production mode (right panel) vs. m h , superimposed on the public exclusion limits of CMS collaboration. For convenience only the red points, i.e. the points passing all of the indirect, LEP and LHC constraints, are plotted here. The results confirm our expectation from figures 6 and 7 that there is no sensitivity in the gluon fusion production mode but many points are above the CMS observed limit in the VBF/VH production mode for a light Higgs boson with mass below 105 GeV.
As the points above the observed CMS upper limit are excluded at 95% C.L., we can expect to exclude some new region in the parameter space thanks to this analysis. To illustrate this point, in figure 10 we plot the points resulting from the previous scan (see JHEP12(2016)068 Figure 9. Points generated in the 2HDM Type I passing indirect, LEP and LHC constraints, superimposed on the results of the CMS 8 TeV low-mass di-photon analysis [11] in the gluon fusion production mode (left panel) and the combined VBF and VH production mode (right panel). The dashed line corresponds to the expected upper limit on σ × BR h→γγ at 95% C.L., with 1 and 2 sigma errors in green and yellow respectively. The solid line is the observed upper limit at 95% C.L. table 8) and passing the three sets of constraints in the plane tan β vs. sin(β − α) (left panel) and in the plane tan β vs. m h (right panel). The violet points have a value of σ VBF/VH × BR h→γγ below the CMS observed upper limit for the corresponding mass; the orange points have a value of σ VBF/VH × BR h→γγ above the CMS observed upper limit and are consequently excluded by the experiment.
The left panel shows that most of the orange points cluster in an exclusion band in the region tan β ∈ [3; 6], sin(β − α) ∈ [−0.27; −0.14]. However, we cannot conclude that the whole orange band is excluded as we have many free parameters: the plot shows in fact a projection of a five-dimensional space on a plane. Therefore, we can have multiple points with a same value of tan β and sin(β − α) but with different values for the other free parameters, producing violet and orange points at the same position in this specific plane. Hence the orange band in the left plot of figure 10 cannot be taken as an absolute exclusion area.
In order to illustrate this point, we produce two additional plots, shown in figure 11, in the plane tan β versus sin(β − α) with all the other free parameters fixed. We choose m h = 87 GeV, m H = 125 GeV, m 12 = 30 GeV and perform this scan for two different values of the mass of the pseudo-scalar and charged Higgs bosons: m A = m H ± = 80 GeV (left panel) and m A = m H ± = 500 GeV (right panel). As before, we only consider points passing the indirect, LEP and LHC constraints. The color code is the same as in figure 10. The exclusion zone does not have the same shape in the two different scans and we can see that the violet points in the left panel are in orange in the panel on the right. It means that we are able to exclude some region in the plane tan β vs. sin(β − α) but the shape and extent of the exclusion zone depends on the value of the other free parameters.

Search for a light pseudo-scalar Higgs boson in the 2HDMs
In the previous section we have seen that values of the pseudo-scalar A masses below 110 GeV are allowed in Type I and Lepton Specific models. It is thus natural to ask if the di-photon resonant signal may be due to the decays of the pseudo-scalar instead of the light scalar h. In this section we will pursue this possibility, limiting ourselves to the same configuration studied above, i.e. fixing the mass of the heavy Higgs boson H to m H = 125 GeV. The constraints on the free parameters of the model coming from JHEP12(2016)068 indirect, LEP and LHC constraints obtained in section 4.2 are also valid in the case of a pseudo-scalar. We can therefore focus on the predicted cross-sections for the pseudo-scalar.
As the kinematic behaviour of the two photons coming from the decay of a pseudoscalar particle is very similar to the the one coming from a scalar particle [40], we can directly apply the CMS study as for the scalar case to constrain a possible light pseudoscalar. The pseudo-scalar A does not couple at tree level to the W and Z bosons, therefore we will only focus on the gluon fusion production mode. Note also that the mass of the other light scalar h is left free, and in principle it can also contribute to the signal at the same time as the pseudo-scalar. To simplify the analysis, however, we will not consider the possible bounds coming from h in this case (as the available parameter space we discuss in the following gives very small cross-section times branching for the pseudo-scalar A which can not be probed at present).

Computation of the cross-section value
The production cross-section of the pseudo-scalar is different from the one for the scalar case. It is clear that for example the effective vertex with the gluons will be different due to different couplings and to the absence of couplings with the gauge bosons. However the "kappa trick" technique used for a scalar can be used here too: where the label SM indicates that the couplings of the pseudo-scalar are set to be equal to the SM couplings of the Higgs boson (except for the different CP properties). However the values of σ SM ggA are not available from the LHC Higgs Cross-Section Working Group and we cannot assume that they are the same as those of the cross-section of the SM scalar Higgs boson σ SM ggh . Furthermore the program 2HDMC does not supply the value of Γ SM A→gg .
JHEP12(2016)068 Figure 13. Production cross-section in gluon fusion mode computed at NNLO by SusHi for an SM scalar particle (in violet) and for an SM pseudo-scalar particle (in green).
We resolve the first issue by obtaining the values of the production cross-section in the gluon fusion mode for a pseudo-scalar with SM-like couplings from SusHi for a discrete set of values and then interpolating between the obtained values to obtain a smooth function. Figure 13 shows the significant difference between the cross-section obtained from SusHi in the gluon fusion production mode at NNLO for an SM scalar particle (in violet) and for an SM-like pseudo-scalar particle (in green) plotted as a function of the mass of the spin-0 particle.
The second issue can be overcome by using an analytical computation. The pseudoscalar A couples to the quarks as g Aqq = g A q ×i mq v ×iγ 5 with g A q = 1 in the SM-like case and g A q = tan β or cot β in the 2HDM case (see table 2). The decay width of a pseudo-scalar A into two gluons can therefore be written at LO as [41]: with τ q ≡ m 2 A /4m 2 q and A A f the fermionic amplitude defined as: The NLO corrections in the heavy top limit can be written as an additional factor to the LO width [41]. Considering only the top and the bottom quarks in the loop, we can then compute the parameter κ 2 g at NLO in Type I: (5.4) Using the values of σ SM ggA from SusHi and the analytic value of κ 2 g given above we are now able to compute the value of σ 2HDM ggA for any possible value of the free parameters. In order to check the validity of the method, we compare the cross-section values obtained with the "kappa trick" method with the ones given by SusHi. We give the results for m h = 87 GeV, m H = 125 GeV, m H ± = 500 GeV, tan β = 8, sin(β − α) = −0.2 and m 12 = 30 GeV in figure 14. In the left panel the mass m A ranges from 60 GeV to 1000 GeV, while the right panel is a zoom in the mass region of interest for the study of a light pseudoscalar. The dashed blue line corresponds to the cross-section computed with the "kappa trick", the dotted red line to the one computed with SusHi and the green solid line to the deviation between the two methods. At low mass the deviation is below 10%, which is low enough with respect to the current uncertainties to be used in an analysis. Above m A = 120 GeV the deviation grows significantly and is about 24% at m A = 1000 GeV. This is due to the NLO corrections in the "kappa trick" which only consider corrections in the infinite top mass approximation. As m A grows, this approximation becomes invalid and the cross-section value diverges from SusHi's results.

Comparison with the CMS low mass di-photon analysis.
The constraints on the free parameters coming from indirect, LEP and LHC constraints obtained in section 4.2 remain valid for the study of a light pseudo-scalar in the scenario where the heavy scalar is identified with the SM-like one at 125 GeV. We can therefore perform a new scan using these bounds, with the additional constraint that the pseudoscalar must have a mass between 80 GeV and 110 GeV in order to fit with the available range of the CMS analysis. The range of variation for the free parameters are given in table 9. We restrict ourselves to Type I only in the gluon fusion production mode.
As for the scalar study we apply the indirect, LEP and LHC constraints. The resulting points are plotted in red in the plane σ gg→A × BR A→γγ vs. m A and superimposed on the CMS results in figure 15.  Table 9. Range of variation for the free parameters used in the study of the pseudo-scalar A. Figure 15. Points generated in the 2HDM Type I passing indirect, LEP and LHC constraints, superimposed on the CMS 8 TeV low-mass di-photon analysis [11] in the gluon fusion production mode. The dashed line corresponds to the expected upper limit at 95% C.L.. The solid line is the observed upper limit at 95% C.L..

JHEP12(2016)068
We can see that the points are well below the CMS observed upper limit on production cross-section times branching ratio at 95% C.L.. We therefore conclude that CMS had no sensitivity to a light pseudo-scalar during the LHC Run 1 in the di-photon final state.

Conclusions
The search for an extended Higgs sector is ongoing at the LHC and represents one of the most important avenues for probing the possible structure of physics beyond the Standard Model. In the simplified setting of Two Higgs Doublet Models, we have explored current constraints from flavour, precision electroweak tests and direct collider searches. We have tested the possible reach of the CMS experiment at the LHC Run 1 for a second Higgs particle lighter than the 125 GeV Higgs boson already discovered. We have explored in detail the different production modes (gluon fusion, vector boson fusion, associated production with a gauge boson) and the subsequent decay to two photons for the light boson. We have found that some sensitivity in these last two production modes is expected even simply recasting an existing Run 1 CMS analysis. A lighter (than the 125 GeV Higgs boson) neutral scalar or pseudo-scalar particle is not completely excluded by present bounds and searches. Out of the four types of 2HDMs, in the low-mass region for a neutral scalar, only Type I has in its parameter space points with large enough cross-section times branching JHEP12(2016)068 ratio to allow detection or exclusion in the gamma gamma decay channel by this analysis. We have applied this analysis also to the case of a light neutral pseudo-scalar, for which however cross-section times branching ratio in the γγ channel is below reach at present. It is however interesting to perform such a low mass analysis (even possibly for lower masses than those considered at Run 1) at 13 TeV for the LHC in Run 2 as the increased sensitivity to lower cross-section values will allow to further explore and constrain or possibly discover new scalar or pseudo-scalar neutral particles and in any case allow a better understanding of an extended Higgs sector.

A Extra numerical results
We list in this appendix some extra numerical results in the form of plots used for the validation of the analysis.  gg→H for the heavy Higgs boson as a function of m H computed with the "kappa trick" (dashed blue line) and with SusHi (dotted red line) and the deviation between the two (solid green line). As the approximation of the infinite mass for the top quark becomes false, the results with "kappa trick" move away from SusHi's results. Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.