Multi-photon production in the Type-I 2HDM

This paper presents a study of a possible contribution to a Higgs boson signal in the $hh\to \gamma\gamma\gamma\gamma$ channel due to $H\to hh$ decays, in the framework of the CP-conserving 2-Higgs Doublet Model Type-I (2HDM-I), where the heavier of the two CP-even Higgs bosons defined herein, $H$, is the SM-like Higgs state observed with a mass of 125 GeV at the Large Hadron Collider (LHC). We perform a broad scan of the 2HDM-I parameter space, in presence of both up-to-date theoretical and experimental constraints, in order to extract the interesting regions yielding such a signal. Then, after validating our numerical framework against public experimental analyses carried out at the LHC, we proceed to assess its scope in constraining and/or extracting the $gg\to H\to hh\to \gamma\gamma\gamma\gamma$ signal in presence of a sophisticated Monte Carlo (MC) simulation. We find that, over a substantial region of the 2HDM-I parameter space presently un-accessible, the LHC will be able to establish such a potential signature in the next 2--3 years.


I. INTRODUCTION
In carrying out the latter, we borrow from existing experimental results. The ATLAS collaboration carried out searches for new phenomena in events with at least three photons at a CM energy of 8 TeV and with an integrated luminosity of 20.3 fb −1 . From the non-observation of any excess, limits are set at 95% CL on the rate of the relevant signal events in terms of cross section multiplied by a suitable BR combination [15] σ BSM × β ≤ 10 −3 σ SM , where β = BR(H → AA) × BR(A → γγ) 2 , σ BSM is the Higgs production cross section in a possible BSM scenario and σ SM is the same, but for the SM Higgs. The above constraint sets an upper limit on β as provided the Higgs state H in the context of new physics phenomena is the SM-like Higgs boson of mass 125 GeV. In particular, we will validate a numerical toolbox that we have created to carry out a Monte Carlo (MC) analysis against results published therein for the case of H → AA decays and extrapolate them to the case of H → hh ones, which constitute the dominant four photon signal in our case. The plan of the paper is as follows. In Section 2, we introduce 2HDMs [3][4][5] in general and describe in particular our construct (the 2HDM-I [13]), including dwelling on the theoretical (see [16][17][18][19][20][21]) and experimental (see later on) constraints placed upon its parameter space. (Herein, we also comment on the fermiophobic limit of the 2HDM-I and its experimental status.) Section 3 is devoted to present our numerical results for the (inclusive) four photons cross section and to motivate the selection of our Benchmark Points (BPs). Then we move on to describe the numerical tools we have used and the MC analysis carried out, including illustrating our results for the exclusive cross section in Section 4. We then conclude in Section 5. Some technical details of our calculations are presented in Appendix A.

II. THE 2HDM-I AND ITS FERMIOPHOBIC LIMIT
The 2HDM scalar potential The most general 2HDM scalar potential which is SU (2) L ⊗ U (1) Y invariant with a softly broken Z 2 symmetry can be written as where φ 1 and φ 2 have weak hypercharge Y = +1 while v 1 and v 2 are their respective Vacuum Expectation Values (VEVs). Through the minimisation conditions of the potential, m 2 11 and m 2 22 can be traded for v 1 and v 2 and the tree-level mass relations allow the quartic couplings λ 1−5 to be substituted by the four physical Higgs boson masses and the neutral sector mixing term sin(β − α), where β is defined through tan β = v 2 /v 1 and α is the mixing angle between the CP-even interaction states. Thus, in total, the Higgs sector of the 2HDM has 7 independent parameters, which include tan β, sin(β − α) (or α), m 2 12 and the four physical Higgs boson masses. As explained in the Introduction, the 2HDM possesses two alignment limits: one with h SM-like [7,8] and an other one with H SM-like [9,10]. In the present study, we are interested in the alignment limit where H is the SM-like Higgs boson discovered at CERN, which implies that cos(β − α) ≈ 1. Then, we take m h < m H /2 ≈ 62.5 GeV, so that the decay channel H → hh would always be open.
From the above scalar potential one can derive the following triple scalar couplings needed for our study: where g is the SU (2) gauge coupling constant. We have used the notation s x and c x as short-hand for sin(x) and cos(x), respectively. It is clear from the above couplings that Hhh is proportional to c β−α which is close to unity in our case and hence the BR(H → hh) would not be suppressed. Moreover, in the exact fermiophobic limit α ≈ ±π/2 becomes proportional to m 2 12 . The vertex HAA has two terms, one proportional to c β−α and the other proportional to s β+α which is close to c β in the fermiophobic limit α ≈ π/2. Finally, the coupling hH ± H ∓ can be large so as to contribute sizably to the h → γγ decay rate.
Fermiophobic limit of the 2HDM-I In general, in the 2HDM, both Higgs doublets can couple to quarks and leptons exactly as in the SM. However, in such case one has tree level FCNCs which would lead to large contribution to B-physics observables in conflict with data. In order to avoid this, the 2HDM needs to satisfy a discrete Z 2 symmetry [11,12] which guarantees the absence of this phenomenon. Several type of 2HDMs exist depending on the Z 2 charge assignment of the Higgs doublets [12]. In our study, we will focus on the 2HDM-I. In this model only the doublet φ 2 couples to all the fermions as in the SM while φ 1 does not couple to any of the fermions.
The Yukawa interactions in terms of the neutral and charged Higgs mass eigenstates in a general 2HDM can be written as: where V ud is the top-left entry of the Cabibbo-Kobayashi-Maskawa (CKM) matrix and P L and P R are the left-and right-handed projection operators, respectively. In the 2HDM-I, we have From the above Lagrangian (7) and (8), it is clear that for α ≈ ± π 2 the tree level coupling of the light CP-even Higgs h to quarks and leptons are very suppressed. Hence h is fermiophobic in this limit [22]. Note that, in the 2HDM-I, the CP-odd Higgs coupling to fermions is proportional to cot β and hence would not be fermiophobic for any choice of tan β. Since we are interested in the case where H is SM-like (cos(β − α) ≈ 1) and m h ≤ m H /2 ≈ 62.5 GeV such that the decay H → hh is open, the main decays of the lightest Higgs state h are into the tree level channels h → V * V * and h → Z * A when m A < m h , otherwise the one loop h → Z * γ and h → γγ ones dominate. The 1 → 4 decays h → V * V * (→ ff f f ) have two sources of suppression, the phase space one and the fact that hV V ∝ sin(β − α) ≈ 0, while the 1 → 3 decay h → Z * γ(→ ff γ) is both loop and phase space suppressed. Therefore, the decays h → γγ and h → Z * A (m A < m h ) are expected to compete with each other and dominate in the fermiophobic limit. In fact, it is well known than h → γγ is dominated by the W ± loops which interfere destructively with top and charged Higgs loops. In the limit where sin(β − α) → 0, the W ± loops vanish and only the top and charged Higgs ones contribute. When cos α vanishes, the h state, with mass ≤ 62 GeV, becomes fermiophobic and consequently the BR(h → γγ) can become 100% if h → Z * A is not open. In contrast, the coupling hZA is proportional to cos(β − α), which is close to unity in our scenario, therefore, when h → Z * A is open for m A < m h , it dominates over h → γγ.
Fermiophobic Higgs bosons have been searched for at LEP and Tevatron. The LEP collaborations used e + e − → Z * → Zh followed by the decay h → γγ and set a lower limit on m h of the order 100 GeV [23][24][25][26]. At the Tevatron, both Higgs-strahlung (pp → V h, V = W ± , Z) and vector boson fusion (qq → q q h) have been used to search for fermiophobic Higgs decays of the type h → γγ [27], with similar results to those obtained at LEP. Note that both LEP and Tevatron assumed a full SM coupling for hV V (V = Z, W ) which would not be the case for the CP-even Higgs h in the 2HDM-I where hV V ∝ sin(β − α) can be very suppressed, as explained. Therefore one could imagine a scenario with a very light h state (m h 60 GeV) which has escaped LEP and Tevatron limits due to suppression in the coupling hV V . In addition, the LEP, OPAL and DELPHI collaborations have searched for fermiophobic Higgs decays through e + e − → Z * → Ah with h → γγ and A decaying mainly into fermions and set a limit on [23,24]. Note that this limit will depend on the coupling ZhA ∝ cos(β − α) and hence becomes weaker for cos(β − α) 1. However, a very light h with m h ≤ 60 GeV is sill allowed if the CP-odd is rather heavy. We refer to [28] for more detail on these aspects. Finally, following phenomenological studies in [29], CDF at Tevatron also studied qq → H ± h, which can lead to 4-photon final states for H ± → W ± h and h → γγ [30], however, CDF limits are presented only for the exactly fermiophobic scenario and are not readily extendable to our more general setup.

III. NUMERICAL RESULTS
As said previously, we are interested in the 2HDM-I for which we perform a systematic numerical scan for its parameter space. We have fixed m H to 125 GeV and assumed that 2m h < m H such that the decay H → hh is open. The other 2HDM independent parameters are varied as indicated in Tab. 1. We use the 2HDMC (v1.7.0) [31] public program to calculate the 2HDM spectrum as well as various decay rates and BRs of Higgs particles. The 2HDMC program also allow us to check several theoretical constraints such as perturbative unitarity, boundedness from below of the scalar potential as well as EW Precision Observables (EWPOs) which are all turned on during the scan. In fact, it is well known that EWPOs constrain the splitting between Higgs masses. In our scenario, since we ask that m H = 125 GeV and assume 2m h < m H , if we want to keep the CP-odd also light, it turns out that the charged Higgs boson would be also rather light, m H ± ≤ 170 − 200 GeV [32], as it can be seen from Tab. 1. Moreover, the code is also linked to HiggsBounds [33] and HiggsSignals [34] that allow us to check various LEP, Tevatron and recent LHC searches.
Once in the 2HDM-I the decay channels H → hh and/or H → AA are open, the subsequent decays of h and/or A into fermions, photons or gluons will lead either to invisible H decays that can be constrained by present ATLAS and CMS data on the Higgs couplings. In our study, we will use the fact that the total BR of the SM-like Higgs boson into undetected BSM decay modes is constrained, as mentioned, by BR BSM < 0.34 [14] where BR BSM will designate in our case the sum of BR(H → hh) and BR(H → AA).
In what follow, we will show our numerical results via three different scans (see Tab. 1). These results mainly concern the BR BSM , BR(H → hh) and the ensuing total cross section for four photons final states which is given by Note that, in writing the above cross section, we have used the narrow width approximation for the SM-like Higgs state H which is justified since the total width of H is of the order of few MeV (see Tab. 3). Furthermore, we are interested into multi-photon signatures coming from H → hh → 4γ and not from H → AA → 4γ decay. In the former case, in fact, because h can become totally fermiophobic, its BR(h → γγ) can in turn become maximal when h → Z * A is not open. In the latter case, of the CP-odd Higgs state A, couplings to fermions are proportional to 1/ tan β, which thus does not vanish, therefore the one loop decay A → γγ will be suppressed compared to the tree level ones A → f f and A → Z * h. We first show our results for σ 4γ without imposing constraint from ATLAS searches in events with at least three photons in the final state [15]. The results of scan-1 are shown in Fig. 1. In this scan we only allow H → hh to be open and deviate from the exact fermiophobic limit by taking α = ±π/2 ∓ δ where δ ∈ [0, 0.05]. It is clear that for some values of δ ≈ 0, one can have an exact fermiophobic Higgs with maximal BR for h → γγ. In this scenario, The BR(H → hh) can reach 17% in some cases. Thus, the four photon cross sections can become of the order of few pb when BR(h → γγ) is close to maximal and BR(H → hh) large. Here, the maximum cross section is reached for The output of scan-2, which is for the exact fermiophobic limit, α = π/2, is shown in Fig. 2. Here, we illustrate σ 4γ as a function of sin(β − α) in the left panel with m h coded with different colours on the vertical axis. The BR(h → γγ) as a function of sin(β − α) is depicted on the right panel of Fig. 2 with the BR(H → hh) on the vertical axis. The maximal value reached by BR(H → hh) in this scenario is again around 17%. Note that, in this case of exact fermiophobic limit, only W ± and H ± loops contribute to the h → γγ decay. In fact, in most cases, W ± loop contributions to h → γγ dominate over the H ± ones except for small sin(β − α), where W ± and H ± terms could become comparable and interfere destructively. In such a case, it may be possible that BR(h → γγ) will be suppressed and BR(h → W * W * ) slightly enhanced. This could explain the drop of the BR(h → γγ) up to 5.5 × 10 −1 . Note that, for large sin(β − α) ≈ −0.14, the off-shell decay h → V * V * can reach 3.5 × 10 −1 , 1 × 10 −1 for V = W and V = Z, respectively. It is interesting to see that, in this scenario, σ 4γ can be larger than 1 pb for a light m h ∈ [10, 50] GeV with significant BR(h → γγ).
In scan-3, we allow sin(β − α) ∈ [−0.35, 0] and the CP-odd Higgs state to be as light as 10 GeV. In this case, we specifically take into account constraints from the LEP measurement of the Z width. The results are illustrated in Fig. 3, where we show both σ 4γ and BR(h → γγ) as a function of sin(β −α). Note that, for any choice of tan β, one can tune sin(β − α) such that α becomes ±π/2 and then h is fermiophobic (in which case the previous discussion of scan-2 would apply again). Away from this fermiophobic limit, BR(h → bb) becomes sizeable and suppresses BR(h → γγ).
After quantifying the maximal size of the four photons cross section in the previous plots, we proceed to apply ATLAS limits coming from searches in events with at least three photons in the final state [15]. The results are shown in Fig. 4. The solid black line is the expected upper limit at 95% CL from ATLAS with 8 TeV centre-of-mass energy and 20.3 fb −1 luminosity. The green and yellow bands correspond, respectively, to a ±1 σ and ±2 σ uncertainty from the resonance search assumption. As it can be seen from this plot, for m h in the [10,62] GeV range, the ATLAS upper limit on σ H × BR(H → hh) × BR 2 (h → γγ) is 1 × 10 −3 σ SM . We also illustrate on this figure our projection for 14 TeV using 300 fb −1 luminosity (based on a MC simulation that we will describe below). The dots represent our surviving points from scan-1, scan-2 and scan-3 after passing all theoretical and experimental constraints. Most of the points with significant four photons cross section and large BR(H → hh) and/or BR(h → γγ) shown in the previous plots turn out to be ruled out by the aforementioned ATLAS upper limit [15]. It is clear from Fig. 4 (top-left and -right plots) that scenarios from scan-1 and scan-2 would be completely ruled out (or, conversely, be discovered) by our projection for the 14 TeV LHC run with 300 fb −1 luminosity while scenarios from scan-3 (bottom plot) would survive undetected. The maximal four photons cross section we obtain is of the order of 37 fb. It is interesting to note that, for scan-1 and scan-3, the remaining points still enjoy sizable BR(H → hh) while for the exact fermiophobic limit of scan-2 one can see from Fig. 4 (top-right) that the BR(H → hh) is less than 5 × 10 −3 . This limit is much stronger than the one from invisible SM-like Higgs decays discussed previously. This can be seen in Fig. 5, where we illustrate the correlation between BR(H → γγ) and BR(H → hh) for the three scans. Herein, one can verify that, for scan-1 and scan-3, BR(H → γγ) and BR(H → hh) are anti-correlated.
Based on the results of these three scans we have selected a few Benchmark Points (BPs), which are given in Tab. 2. These BPs can be seen in Fig. 4 ), BR(h → Z * A), BR(A → γγ) and the four photons cross section σ 4γ in fb. In fact, for these BPs, we take into account all theoretical constraints as well as the LEP and LHC constraints implemented in the HiggsBounds code plus the limits from ATLAS on multi-photons final states [15], as explained in the introduction, see Eqs. (1) and (2). It is also interesting to see from Tab. 3 that the BR(A → γγ) is always suppressed and cannot be used to generate multi-photon finale states. Finally, it can also be seen from this table that, even for small BR(H → hh) ≈ 10 −3 but with maximal BR(h → γγ), one can still get a large σ 4γ ≈ 38 fb.
Before ending this section, we would like to comment on charged Higgs and CP-odd Higgs boson searches. As mentioned previously, the charged Higgs and CP-odd Higgs states are rather light in our scenarios. LHC limits on light charged Higgs states produced from top decay and decaying to H± → τ ν, cb in the 2HDM-I can be evaded by advocating the dominance of the H ± → W ± A or H ± → W ± h BRs (see [28] for more details). On the one hand, the LHC searched for a CP-odd Higgs state decaying via A → ZH [35][36][37] and A → τ + τ − . In our scenario, the BR(A → ZH) will suffer two suppressions: one coming from the coupling AZH, which is proportional to sin(β−α) ≈ 0, and the other one coming from the fact that A → Zh would dominate over A → ZH since h is lighter than 125 GeV and the coupling ZAh is proportional to cos(β − α) ≈ 1. On the other hand, ATLAS and CMS searches for a CP-odd Higgs state decaying to a pair of τ leptons [38,39], when applied to the 2HDM-I, only exclude small tan β ≤ 1.5 for m A ∈ [110, 350] GeV [40]. This can be understood easily from the fact that A couplings to a pair of fermions in the 2HDM-I are proportional to 1/ tan β, hence both the production gg → A and the decay A → τ + τ − are suppressed for large tan β values. Moreover, in our scenario, BR(A → τ + τ − ) would receive an other suppression from the opening of the A → Z * h channel. Note also that LEP limits on a light h and a light A are implemented in the HiggsBounds code through limits on the processes e + e − → Zh and e + e − → hA.

IV. SIGNAL AND BACKGROUND
As previously discussed, in Figs. (4-5, we have taken into account the constraints from the ATLAS collaboration reported in [1] from 8 TeV data. However, in order to project the sensitivity of the future LHC run at √ s = 14 TeV, we have to rescale these results. To determine the 'boost factors', for both signal and background processes, needed to achieve this, we resort to the MC tools. Specifically, we generate parton level events of both signal and background processes by using MadGraph 5 [41] and then pass them to PYTHIA 6 [42] to simulate showering, hadonisation and decays. We finally use PGS [43] to perform the fast detector simulations.
In order to pick out the relevant events, for our Run 1 analysis, we adopt the same selection cuts of the ATLAS collaboration given in [1], which read as follows.
• The two leading photons should have a P t (γ) > 22 GeV and the third one should have a P t (γ) > 17 GeV.  • The photons should be resolved in the range |η| < 2.37 and do not fall in the endcap region 1.37 < |η| < 1.52.
• The cone separation parameter ∆R(γγ) between a pair of photons should be larger than 0.4.
One interesting observation is that the kinematics of photons from the processes gg → H → hh → 4γ and that of gg → H → AA → 4γ are similar when m h = m A , which could be attributed to the fact that, although h and A have different parity, the differential cross sections of these two processes are both proportional to (k 1 · k 2 ) 2 (k 3 · k 4 ) 2 (plus permutations) after the sum over photon polarisations (the k i 's, for i = 1, 2, 3, 4, are the photon four-momenta). We provide an Appendix to demonstrate the details. In Fig. 6, we expose the similarity between the two processes by showing some kinematic spectra of gg → H → hh → 4γ and gg → H → AA → 4γ. In particular, we present the m 3γ spectrum (the invariant mass of the three leading P t -ordered photons) as well as the m 23 spectrum (the invariant mass of the 2nd and 3rd P t -ordered photons). Obviously, these spectra show no significant difference for these two processes gg → H → hh → 4γ and gg → H → AA → 4γ, except fluctuations from numerical simulation. Therefore, the experimental methods and results of multi-photon data from gg → H → AA → 4γ can also be applied to gg → H → hh → 4γ.
In order to establish LHC sensitivity to our signal process, we determine the scaling factors for both signal and backgrounds, necessary to map our own MC simulations onto the real data results of ATLAS. In doing so, we use the Leading Order (LO) cross section to determine such factors. Therefore the latter should encode the K-factors due to higher order corrections, the difference between real detector and the fast detector simulations, the mistaging rate of a jet as a photon and an electron as a photon as well, etc. In fact, when the rejection rate of a jet as a photon is considered [44], the fake rate could be around 10 −3 , so for the process γγ + j the scaling factor is dominantly determined by the fake rate while for the process γ + jj we expect the fake rate to be around 10 −6 , the scaling factors demonstrating then that significant background contribution is indeed due to fake rates.
The scaling factors for each process are listed in Tab. V, as mentioned, being all determined from the aforementioned ATLAS results at 8 TeV. We also compare the experimental line-shapes with those from our MC events, which are plotted in Fig. 7. Although the spectra from MC are slightly harder and noticeable differences appear in the bins with m 3γ < 50 GeV and m 2γ < 50 GeV, the total number of predicted events is close to the experimental ones.
By assuming the same scaling factors, we examine the boost factor in the LHC sensitivity for an increased collision energy of √ s = 14 TeV. The ensuing cross sections of the signal and background processes are given in Tab. VI. Since the signal production process gg → H has a larger boost factor when the collision energy increases from 8 TeV to 14 TeV due to the larger enhancement of gluon flux, as compared to the more varied background composition, it is natural to expect a better sensitivity for the future LHC runs, as readily seen in the table. For example, when the integrated luminosity of the LHC is assumed to be 300/fb, the boost factor in cross section (which is defined in the caption) is found to be 32.2 for thew signal and 25.7 for the background. This effect reflects then in the projected sensitivities shown in Fig. 4 for the LHC with √ s = 14 TeV and 300/fb of luminosity (blue lines).

CONCLUSIONS
In this paper, we built upon previous results of some of ours, which had extracted a region of parameter space of the 2HDM-I where very light h and A states, down to 15-20 GeV or so, can be found, when the H one is assumed to be the SM-like one discovered at the LHC in 2012. This spectrum is well compatible with all standard theoretical constraints (unitarity, vacuum stability, etc.) and all available experimental data (including flavour as well as Higgs data) and thus offers the possibility of testing Higgs cascade decays of the type H → hh and H → AA, compatibly with the total H width extracted by global fits to the 125 GeV Higgs data. Amongst the possible decays of the h and A states we concentrated here upon those yielding di-photons, the overall signature then being a 4γ one, primarily induced by gg → H creation. We do so as the 2HDM-I can develop, over the aforementioned region of parameter space, a (nearly) fermiophobic limit, so that h and A decays into fermions (chiefly, bb and τ + τ − ) are negligible. In fact, the availability of an ATLAS analysis performed on Run 1 samples of the LHC looking for these specific multi-photon signals allowed us, on the one hand, to validate our MC tools against the full detector environment of a multi-purpose  LHC experiment and, on the other hand, to project our finding into the future by extrapolating our results to a collider energy of 14 TeV and luminosity of 300/fb. This exercise revealed that the portion of 2HDM-I parameter space where the above phenomenology is realised, while being just below the current LHC sensitivity, is readily accessible at future stages of the LHC. To confirm or disprove its existence is of paramount importance as this would almost univocally point to a specific realisation of a generic 2HDM construct as such light and fermiophobic h and A states cannot be realised in alternative formulations of it.  In this appendix, we demonstrate in details that the process gg → H → hh → 4γ and the process gg → H → AA → 4γ have the same differential cross section, except an overall factors.
For the CP-even Higgs case, the matrix element of the process gg → H → hh → 4γ can be put as iM = iC(k 1 · k 2 η µν − k µ 2 k ν 1 ) * µ (k 1 ) * ν (k 2 )(k 3 · k 4 η ρσ − k ρ 4 k σ 3 ) × * ρ (k 3 ) * σ (k 4 )δ ab (p 1 ) · (p 2 ) , where p 1 and p 2 is the momentum of the initial gluons, k 1 − k 4 are momentum of 4 photons in the final state. Note that we have group the effective couplings of each vertices and the propagators together in C and just show their Lorentz structure and color indices. The squared amplitude with the average of the degree of freedom of initial states By using the polarisation sum formula * µ ν = −η µν , and the fact δ ab δ ab = 8, we arrive at Due to the on-shell conditions for the photons of the final state, we have k 2 1 = k 2 2 = k 2 3 = k 2 4 = 0, then we obtain the following squared matrix element For the CP-odd Higgs case, the matrix element of the process gg → H → AA → 4γ can be put as where D includes all effective couplings and the propagators. The squared matrix element with a sum over the