Searching for Heavier Higgs Boson via Di-Higgs Production at LHC Run-2

The LHC discovery of a light Higgs particle $h^0$ (125GeV) opens up new prospect for searching heavier Higgs boson(s) at the LHC Run-2, which will unambiguously point to new physics beyond the standard model (SM). We study the detection of a heavier neutral Higgs boson $H^0$ via di-Higgs production channel at the LHC (14TeV), $H^0 \to h^0h^0 \to WW^*\gamma\gamma$. This directly probes the $Hhh$ cubic Higgs interaction, which exists in most extensions of the SM Higgs sector. For the decay products of final states $WW^*$, we include both pure leptonic mode $WW^* \to \ell\bar{\nu}\bar{\ell}\nu$ and semi-leptonic mode $WW^* \to q\bar{q}'\ell\nu$. We analyze signals and backgrounds by performing fast detector simulation for the full processes $pp \to H \to hh \to WW^*\gamma\gamma \to \ell\bar{\nu}\bar{\ell}\nu\gamma\gamma$ and $pp \to H \to hh \to WW^*\gamma\gamma \to \ell\nu q\bar{q}'\gamma\gamma$, over the mass range $M_H=250-600$GeV. For generic two-Higgs-doublet models (2HDM), we present the discovery reach of the heavier Higgs boson at the LHC Run-2, and compare it with the current Higgs global fit of the 2HDM parameter space.


Introduction
Most extensions of the standard model (SM) require an enlarged Higgs sector, containing more than one neutral Higgs states. After the LHC discovery of a light Higgs particle h 0 (125 GeV) [1][2], a pressing task of the on-going LHC Run-2 is to search for additional new Higgs boson(s), which can unambiguously point to new physics beyond the SM.
Such an enlarged Higgs sector [3] may contain additional Higgs doublet(s), or Higgs triplet(s), or Higgs singlet(s). For instance, the minimal supersymmetric SM (MSSM) [4] always requires two Higgs doublets and its nextto-minimal extension (NMSSM) [5] further adds a Higgs singlet. The minimal gauge extensions with extra SU (2) or U(1) gauge group [6] [7] will invoke an additional Higgs doublet or singlet. The minimal left-right symmetric models [8] include an extra product group SU(2) R ⊗ U(1) B−L , and thus requires a Higgs bidoublet plus two Higgs triplets. For the demonstration in our present LHC study, we will consider generic two-Higgs-doublet models (2HDM) [9] under the SM gauge group. To evade constraints of flavor changing neutral current (FCNC), it is common to impose a discrete Z 2 symmetry on the 2HDM. For different model settings of Higgs Yukawa interactions, the 2HDMs are conventionally classified into type-I, type-II, lepton-specific, neutrino-specific, and flipped 2HDMs [9]. The current study will focus on the conventional type-I and type-II 2HDMs.
For the heavier Higgs state H 0 with mass above twice of the light Higgs boson h 0 , M H > 2M h 250 GeV, the di-Higgs decay channel H → hh is opened and becomes significant, in addition to the other SM-like major decay modes H → WW, ZZ . Hence, the LHC can search for the di-Higgs production channel pp → H → hh [4,6,10,11]. ATLAS analyzed the decay channel hh → bbγγ at the LHC (8 TeV) run and found a 2.4σ excess at M(bbγγ) 300 GeV [12]. CMS performed similar searches for this channel and derived limits on the parameter space [13]. An analysis of this channel at 14 TeV runs with high luminosity 1000 fb −1 was done for 2HDM [14]. Another study considered the SM plus a heavy singlet scalar via H → hh → bbWW * → bb ν ν channel for 14 TeV runs with 3000fb −1 luminosity [15]. We note that it is possible to increase the sensitivity of H 0 searches by studying and combining more decay channels of the di-Higgs bosons.
In this work, we perform systematical study of H 0 production via a new decay channel of di-Higgs bosons, pp → H → hh → WW * γγ . For the final state weak bosons, we will analyze both pure leptonic mode WW * → ν¯ ν and semi-leptonic mode WW * → qq ν . Since a SM-like Higgs boson h 0 (125GeV) has decay branching fractions Br[h → bb, WW * , ZZ * ] (58%, 22.5%, 2.77%), we see that the di-Higgs decay mode hh → WW * γγ (with pure leptonic or semi-leptonic WW * decays) has the advantage of much cleaner backgrounds than hh → bbγγ , while Br[h → WW * ] is only smaller than Br[h → bb] by about a factor of 2.6 . Hence, we expect that the hh → WW * γγ mode should have comparable sensitivity to hh → bbγγ mode, and is more sensitive than hh → bbWW * mode.
This letter paper is organized as follows. In section 2, we present the production and decays of the heavier Higgs boson H 0 in 2HDM type-I and type-II. Then, in section 3, we systematically analyze the signals and backgrounds for the reaction pp → H → hh → WW * γγ , including both pure leptonic and semi-leptonic decay modes of the WW * final state. In section 4, we further analyze the LHC probe of the parameter space for 2HDM-I and 2HDM-II, and compare it with the current Higgs global fit. Finally, we conclude in section 5.

2HDM Setup and Parameter Space
For the present phenomenological study, we consider the 2HDM [9] as the minimal extension of the SM Higgs sector. We set the Higgs potential to have CP conservation, and the two Higgs doublets H 1 and H 2 have hypercharge Y = + 1 2 , under the convention Q = I 3 + Y . It is desirable to assign a discrete Z 2 symmetry to the Higgs sector, under which the Higgs doublet H 1 ( H 2 ) is Z 2 even (odd). With these, the Higgs potential can be written as where the masses and couplings are real, and we have allowed a soft Z 2 breaking mass term of M 2 12 . The minimization of this Higgs potential gives the vacuum expectation values (VEVs), The two doublets jointly generate the electroweak symmetry breaking (EWSB) VEV v 246 GeV, via the relation v = (v 2 1 + v 2 2 ) 1/2 , where v 1 = v cos β and v 2 = v sin β . Thus, the parameter tan β is determined by the Higgs VEV ratio, tan β = v 2 /v 1 . The two Higgs doublets contain eight real components in total, Three imaginary components are absorbed by (W ± , Z 0 ) gauge bosons, while the remaining five components give rise to the two CP-even neutral states (h 0 1 , h 0 2 ), one CP-odd neutral state A 0 , and two charged states H ± . The mass eigenstates (h, H) of the neutral CP-even Higgs bosons are given by diagonalizing the mass terms in the Higgs potential (1). They are mixtures of the gauge eigenstates (h 1 , h 2 ) , where α is the mixing angle. Among the two neutral Higgs bosons, h is the SM-like Higgs boson with mass M h 125 GeV, as discovered at the LHC Run-1 [1] [2], and H is the heavier Higgs state. We will systematically study the LHC discovery potential of H state in the present work. The Higgs potential (1) contains 8 parameters in total, three masses and five couplings. Among these, we redefine 7 parameters as follows: the EWSB VEV v, the VEV ratio tanβ , the mixing angle α, and the mass-eigenvalues (M h , M H , M A , M H± ). We may choose the 8th parameter as the Higgs mass-parameter M 2 12 . Note that once we fix the mass spectrum of the 5 Higgs bosons as inputs, we are left with only 3 independent parameters (α , tanβ) and M 2 12 . The current LHC data favor the parameter space of the 2HDM around the alignment limit [9], under which cos(β − α) = 0 . This limit corresponds to the light Higgs state h to behave as the SM Higgs boson with mass 125 GeV. For practical analysis, we fix M h 125 GeV by the LHC data and vary the heavier mass M H within the range of 250 − 600 GeV. We consider the Higgs states A and H ± to be relatively heavy, within the mass-range M A , M H± = 0.3 − 2 TeV for simplicity. We will scan the parameter space and analyze the LHC production and decays of H in the next section. The heavier neutral Higgs boson H has gauge couplings with (W ± , Z 0 ) and Yukawa couplings with quarks and leptons, which depend on the VEV ratio tanβ and mixing angle α . The gauge couplings of H with V (= W, Z) differ from the SM Higgs coupling by a scaling factor cos(β−α) , The Yukawa interactions of H with fermions can be expressed as follows, where the dimensionless coefficient ξ f H differs between the Type-I and Type-II of 2HDM, as summarized in Table 1.
where in the second step we have used the relation M 2 A + λ 5 v 2 = 2M 2 12 /sin2β . In the SM, the cubic Higgs coupling G sm hhh = −3M 2 h /v . We define a coupling ratio, ζ = G Hhh /G sm hhh , which characterizes the relative strength of the Hhh coupling as compared to the h 3 Higgs coupling of the SM. Under alignment limit cos(β − α) → 0 , the trilinear scalar coupling (6) takes the asymptotical form, In Fig. 1  by the factor ζ 2 . In Fig. 1, the red points present the viable parameter space consistent with vacuum stability, unitarity and perturbativity bounds of the Higgs potential [9]. We also take into account the 3σ constraints from the current Higgs global fit (cf. Sec. 4). The electroweak precision data also constrain the parameter space of the 2HDM. It was found that in the 2HDM the charged Higgs mass satisfies, −600 GeV < M H± − M 3 < 100 GeV and M H± > 250GeV [16], where M 3 is the mass of the heaviest neutral scalar. In the case with exact Z 2 (M 12 = 0), the potential could be valid up to the scale ∼ 10TeV [17], while for the present case of a softly broken Z 2 , the bound is much more relaxed, and the theory can be valid up to the Planck scale. For the analysis of In the following analysis, we will consider the same range of the 2HDM parameter space unless specified otherwise.

Heavier Higgs Boson H 0 : Decays and Production
Let us consider the decay modes of the heavier neutral Higgs boson H 0 . It is straightforward to infer the tree-level decay width for M H > 2M h ,  . We see that for M H < 250 GeV, the dominant decay channels are H → ZZ, WW, and for 250 GeV < M H < 350 GeV, the major decay channels include H → ZZ, WW, hh since the H → hh channel opens up. For M H > 350GeV, the H → tt channel is further opened, and will become dominant in 2HDM-II when cos(β−α) takes values around the alignment limit as shown in Fig. 2(b). But this situation can change when cos(β−α) becomes larger and falls into the allowed region which separates from the alignment region (cf. Fig. 9 in Sec. 4).
From Eq. (5) and Tabel 1, we see that the Yukawa coupling of the heavier Higgs boson H with tt has a scale factor ξ t H = sinα/sinβ relative to the SM Higgs Yukawa coupling. The major LHC production channel is the gluon fusion process gg → H . Other production processes include the vector boson fusion pp → Hqq , the vector boson associated production pp → HV, and the top associated production gg → Htt . The gluon fusion production cross section of H can be obtained from the corresponding SM cross section with a rescaling by H → gg partial width, where we will include all NLO QCD corrections to the gluon fusion cross section as done in the SM case [19]. We note that for 2HDM-I, Table 1 shows that the H Yukawa couplings with top and bottom quarks have the same structure as in the SM, so the H production cross section σ[gg → H] differs from the SM by a simple rescaling factor (sinα/sinβ) 2 . For the 2HDM-II, we see that the H coupling to b quarks differs from that of t quarks by a factor of tan β/tanα , which can enhance the b-loop contribution to gg → H production for large tan β region. Hence, the general relation (9) should be used. The uncertainty of the gluon fusion cross section is about 10% over the mass-range M H = 250 − 600 GeV [19], which is roughly the total uncertainty of signal and background calculations.
For the inclusive H production, we include the gluon fusion gg → H, and b-related processes bb → H, gb (gb) → Hb (Hb), and gg (qq) → Hbb . The production cross sections for these b-related processes are derived by rescaling a factor of (ξ d H ) 2 from their corresponding SM productions with the same Higgs mass. So we have the total inclusive cross section of pp → HX for the 2HDM, We present the inclusive H production rate for 2HDM Type-I and Type-II in Fig. 3(a)-(b). Multiplying the production cross section with decay branching fraction Br(H → hh → WW * γγ), we compute the signal rate in the channel 6 pp → HX with H → hh → WW * γγ . We summarize our results in Fig. 4 for 2HDM-I and 2HDM-II, respectively. In Fig. 3 In plots (c)-(d), the red curve (gg → H contribution) fully overlaps the black curve (summed total contribution). This is because the gluon fusion channel gg → H dominates the inclusive production cross section for low tan β region of the 2HDM. In general, Table 1 shows that for 2HDM-I, the H Yukawa couplings ( ξ u H = ξ d H = sin α/sin β ) are rather insensitive to tan β . Hence, in 2HDM-I the gluon fusion actually dominates the H production over full range of tan β 1, and the contributions of b-related sub-channels are always negligible. For 2HDM-II, the (up-type) H Yukawa coupling ξ u H = sin α/sin β is the same as 2HDM-I, and the down-type Yukawa coupling ξ d H ∝ 1/cos β = tan β/sin β is enhanced by a tan β factor relative to ξ u H . We find that for small tan β 3 , the gluon fusion channel still dominates the inclusive H production in 2HDM-II, and its cross section is larger than other b-related channels by a factor of O(10−100) for M H 300 GeV. The analysis of 2HDM-II in Sec. 4 also concerns the small tan β region [cf. Fig. 9(b)(d)]. Hence, in the following Sec. 3-4, we will focus our analysis on the Higgs production from gluon fusion channel, pp(gg) → H → hh → WW * γγ .

Higgs Signal and Background Simulations
In this section, we compute the Higgs signals and backgrounds at the LHC (14 TeV). We perform systematical simulations by using MadGraph5 package [21] for the process, pp(gg) → H → hh → WW * γγ, via gluon fusion channel. The parton-level Higgs production cross section σ(gg → H) is derived from Eq. (9), including NLO QCD corrections. We illustrate the signal Feynman diagram by the left plot of Fig. 5. For signal process, we generate the model file using FeynRules [22], containing Hhh vertex and the effective ggH vertex. We compute signal and background events using MadGraph5/MadEvent [21]. Then, we apply Pythia [23] to simulate hadronization of partons and adopt Delphes [24] to perform detector simulations.
For the final state WW decays, we will study both the pure leptonic mode WW → ν ν and the semi-leptonic mode WW → qq ν. The W decay branching fractions to eν and µν equal 10.8% and 10.6%, respectively, while Figure 5: LHC production process gg → WWγγ . The left diagram shows the signal production via gg → H → hh → WWγγ , and the right diagram illustrates an irreducible background process gg → qq νγγ . Table 2: Signal and background cross sections of pp → WW * γγ → ν νγγ and pp → WW * γγ → qq νγγ processes at the LHC (14 TeV) after each set of cuts. The signal significance(Z 0 ) is computed for the LHC (14 TeV) runs with 300 fb −1 integrated luminosity. We input the heavier Higgs mass M H = 300 GeV, and set the sample signal cross section as σ(pp → H → hh → WW * γγ) = 5 fb . From the 3rd to 5th columns, we show the signals and backgrounds after imposing each set of cuts. The "Selection + Basic Cuts" are choosing according to Eqs. (11)- (12). In the pure leptonic mode, we impose the Final Cuts M T ( νν), M( ), M T ( ννγγ), ∆φ( ), ∆R( ), and ∆R(γγ). In the semi-leptonic mode, we add the Final Cuts P T (γ) M T (qq ν), and ∆R(γγ). that of W → τν is about 11.3% [25]. The dijet branching ratio of W → qq equals 67.6% [25]. Thus, the inclusion of semi-leptonic mode will be beneficial. Since τ leptons can decay into e, µ , the detected final state e, µ will include those from the τ decays. For M h = 125 GeV, the branching fraction of h → γγ in the SM equals 2.3 × 10 −3 [18].
In the following, we will first present the analyses for M H = 300 GeV in Sec. 3.1-3.2, and then for heavier masses M H = 400, 600 GeV in Sec. 3.3.
3.1. Pure Leptonic Channel: hh → WW * γγ → ν νγγ For pure leptonic channel, we have hh → WW * γγ → ν νγγ . Although this channel has an event rate about two orders of magnitude lower than that of hh → bbγγ mode, it has much cleaner background as compared to bbγγ final state. After imposing simple cuts, we find that the backgrounds can be substantially reduced. We follow the ATLAS procedure for event selections. To discriminate the Higgs signal from backgrounds, we set up preliminary event selection by requiring two leptons (electron or muon) and at least two photons in the final state, In the first step of event analysis, we need to prevent the potential double-counting, i.e., the reconstructed objects are required to have a minimal spatial separation [26]. The two leading photons are always kept, but we impose the following criteria [26]: (i) electrons overlapping with one of those photons within a cone ∆R(e, γ) < 0.4 are rejected; (ii) jets within ∆R(jet, e) < 0.2 or ∆R(jet, γ) < 0.4 are rejected; (iii) muons within a cone of ∆R(µ, jet) < 0.4 or ∆R(µ, γ) < 0.4 are rejected. After this, we apply the basic cuts to take into account the detector conditions, which are imposed as follows, P T (γ), P T (q) > 25 GeV, P T ( ) > 15 GeV, |η(γ)|, |η(q)|, |η( )| < 2.5 .
Next, we turn to the background analysis for pure leptonic mode. Besides the ν νγγ and γγ backgrounds, there are additional reducible backgrounds from Higgs bremstrahlung, vector boson fusion, and tth production. The cross section of the former two processes are fairly small and thus negligible for the present study. The tth associate production, with tt → WWbb , can be important because the diphoton invariant-mass cut does not effectively discriminate the signal process. But, this background can be suppressed by imposing b-veto [27]. The production cross section for tth in the SM is σ(pp → tth) = 0.6113 pb [28]. The latest b-veto efficiency of ATLAS is, (b−veto) = 22% [29]. Thus, we estimate the cross section for this background process, where W → ν includes = e, µ, τ. We see that imposing the b-veto has largely suppressed the tth background. We note that the tth background is much smaller than the ν νγγ background before kinematic cuts, while after all the kinematic cuts it could be non-negligible. So we will include both for the present background analysis.  Another potential background may arise from the Higgs pair production pp → hh in the SM [30,31,32]. Our signal process pp → H produces on-shell Higgs boson H with decays H → hh , which has much larger rate as well as rather different kinematics from the non-resonant di-Higgs production in the SM. (Since our signal has on-shell H production, we find that its interference with the SM-type non-resonant hh production is negligible after kinematical cuts.) For instance, we can further suppress this SM di-Higgs contribution by imposing a cut on the transverse mass of di-Higgs bosons.
We also consider a reducible background from the Zh associate production. The SM cross section of this process pp → Zh at the LHC is σ(pp → Zh) = 0.761 pb [18]. Hence, this background gives σ(pp → Zh → γγ) = 0.175 fb before any cuts. Because the Zh background must have the invariant mass M( ) of final state di-leptons peaked at M Z 91.2 GeV, we can efficiently kill this background by applying a narrow cut on M( ) , which has little effect on the signal rate. In the present analysis, we choose, M( ) ∈ (M Z − 5Γ Z , M Z + 5Γ Z ), where Γ Z 2.5 GeV is the total width of Z boson. Other reducible backgrounds come from the fake events in which quark and/or gluon are misidentified as photons. These backgrounds include ν νqγ, ν νgγ, ν νqq, ν νqg, and ν νgg . For our analysis, we adopt the fake rates used by ATLAS detector [33], With such small fake rates, we find that these reducible backgrounds are negligible. In summary, with the above considerations of the SM backgrounds, we will compute the irreducible backgrounds with final state ν νγγ , and the reducible backgrounds including the γγ final state, the tth associate production, the Zh associate production, and the SM di-Higgs production.
In Fig. 6, we present the distributions of relevant kinematical variables for the pure leptonic channel, including both signals and backgrounds. In this figure, we show the signal distributions at the LHC (14TeV) with 300 fb −1 integrated luminosity for M H = (300, 400, 600) GeV by (red, green, blue) curves as well as the backgrounds (black curves). Here we have input the sample cross section σ(pp → H → hh → WW * γγ) = (5, 3, 1) fb for M H = (300, 400, 600) GeV, respectively. In the following, we will analyze how to effectively suppress the SM backgrounds by implementing proper kinematical cuts.
From Fig. 6(a)-(b), we first impose kinematical cuts on the diphotons invariant-mass M γγ and the missing energy E / T of final state neutrinos, 120 GeV < M γγ < 130 GeV, The missing energy cut can also sufficiently remove the γγ background. Then, inspecting Fig. 6(c)-(f), we apply the kinematical cuts on the azimuthal angle ∆φ and opening angle ∆R for the final state di-leptons and di-photons, respectively, Here, from the distributions of Fig. 6(c), we find that the ∆φ(γγ) cut is not effective for Higgs mass M H = 300 GeV. So we do not implement this cut. For the transverse mass cut [25], we consider the transverse mass M T for the νν system with two leptons and missing energy, which should be no larger than the Higgs mass M h 125 GeV. All the final state leptons/neutrinos are nearly massless, so the transverse energy of each final state equals its transverse momentum E T,i | P T,i | , (i = 1, 2, 3), where i = 1, 2 denote two leptons 1,2 and i = 3 denotes the system of two neutrinos. Thus, we have With this and inspecting Fig. 6(g), we implement the transverse mass cut, From Fig. 6(h), we will further impose the transverse mass cut for the full final state ννγγ , 60GeV < M T ( ννγγ) < 320 GeV .
The kinematical cuts for the cases of M H = 400 GeV and 600 GeV will be discussed in Sec. 3.3.
We summarize the results in Table 2 for both signal and backgrounds. For demonstration, we first input the heavier Higgs mass M H = 300 GeV, and set the sample signal cross section σ(pp → H → hh → WW * γγ) = 5 fb for the LHC (14 TeV). In Table 2, we also show the significance of signal over backgrounds after each set of kinematical cuts at the LHC Run-2 with 300 fb −1 integrated luminosity. When the event number is small, we can use the median significance(Z 0 ) (instead of S / √ B ), as defined in following [34], As shown in Table 2, after applying all kinematical cuts, we estimate the signal significance(Z 0 ) = 5.15 .

Semi-leptonic Channel
The analysis of semi-leptonic channel WW * → qq ν is similar to that of the pure leptonic mode WW * → ν ν . But, there are nontrivial differences. One thing is that for each decay we need to specify which decay mode is from onshell W (qq or ν), since these two situations have different distributions. To illustrate this, we present the distribution of M qq in Fig. 7(a), where the green (blue) curve depicts the final state qq from on-shell (off-shell) W decays, and the red curve represents the actual distribution of M qq from WW * → qq ν . Fig. 7(a) shows that the M qq distribution from on-shell W decays (green curve) has event rate peaked around M qq = 70 − 80 GeV , while the M qq distribution from off-shell W decays (blue curve) is rather flat.
Our first step here is also to remove the pileup events, similar to Sec. 3.1. Then, we select the final states by imposing the preliminary cuts n j 2 , n γ 2 , n = 1 .
For jets we choose the leading and subleading pair, while for photons we choose the diphoton pair whose M γγ is closet to M h = 125 GeV. Then, we choose the basic cuts to be the same as in Eq. (12).
Next, we turn to the background analysis. The most important background for this channel comes from the SM irreducible background, pp → qq νγγ , whose cross section is about σ[qq νγγ] 31.6 fb. Another significant reducible background is the SM process pp → νγγ , which has a cross section σ[ νγγ] 143 fb . But this will be mainly rejected by the jet-selections in Eq. (21). For the tth background, we find that under b-veto its cross section is 0.0148 fb, as shown in Table 2. Single top associated Higgs production gives another background, σ[pp → th(th) + X] = 79.4 fb [35], where X represents single-jet or dijets in our simulation. We find that under b-veto this cross section of pp → th(th) + X → b νγγ + X reduces to about 0.013 fb . We also include the non-resonant di-Higgs production in the SM, which has much smaller event rate and rather different kinematics. Other potential SM backgrounds may include the reducible backgrounds such as qq νgg with gg misidentified as γγ . This is actually negligible due to the tiny g → γ misidentification rate shown in Eq. (14).
For the kinematic cuts, we choose the M γγ cut as in Eq. (15). The invariant-mass M qq should match the W mass. We depict the M qq distribution in Fig. 7. Plot-(a) depicts the decay mode with on-shell (off-shell) decays W (W * ) → qq by green (blue) curve, for M H = 300 GeV. The realistic decays of WW * → qq ν correspond to the red curve. In plot-(b), we present the M qq distribution for full signals of WW * → qq ν by (red, green, blue) curves for M H = (300, 400, 600) GeV. The black solid curve in each plot gives the full backgrounds. From Fig. 7, we choose the M qq cut, M qq < 250 GeV. (22) (GeV) Signal (300GeV)  Signal Signal (300GeV)  We present the distributions for other kinematical observables in Fig. 8, where we have input the sample cross section σ(pp → H → hh → WW * γγ) = (5, 3, 1) fb for M H = (300, 400, 600) GeV. From Fig. 8(a)-(b), we impose cuts on the diphoton invariant-mass M γγ and the missing energy E / T of final state neutrinos, 120 GeV < M γγ < 130 GeV, 10 GeV < E / T < 80 GeV.
We require E / T > 10 GeV to suppress certain reducible backgrounds, as also adopted in the ATLAS analysis. For instance, consider the background qqγγ+ j with j mistagged as a lepton, where j denotes a gluon or quark jet. Since Table 3: Signal and background cross sections of both pp → WW * γγ → ν νγγ and pp → WW * γγ → qq νγγ processes at the LHC (14 TeV) after each set of cuts. The signal significance(Z 0 ) is computed for the LHC (14 TeV) runs with 300 fb −1 integrated luminosity. We input the heavier Higgs mass M H = 400 GeV, and set the sample signal cross section σ(pp → H → hh → WW * γγ) = 3 fb . From the 3rd to 5th columns, we present the signals and backgrounds after imposing each set of cuts. In the pure leptonic mode, we impose the Final Cuts M T ( νν), M( ), M T ( ννγγ), ∆φ( ), ∆R( ), ∆φ(γγ), and ∆R(γγ). In the semi-leptonic mode, we add the Final Cuts P T (γ), M T (qq ν), ∆φ(γγ), and ∆R(γγ). it contains no neutrino in the final state, we can eliminate it by imposing the missing energy E / T cut. This is more like a basic cut. For the transverse momentum distribution of the leading photon shown in Fig. 8(d), we set the following cut, 60 GeV < P T (γ) < 150 GeV.
Then, we inspect the transverse mass distribution of qq ν final state, which arises from the decay products of h → WW * → qq ν . From Fig. 8(c), we impose the following cut, With Fig. 8(e)-(f), we have also examined possible cuts on ∆φ(γγ) and ∆R(γγ) distributions. We further impose, We summarize our results in Table 2. Here we present the signal and background cross sections after each set of cuts. We take an integrated luminosity of 300 fb −1 for the LHC (14 TeV), and derive the corresponding signal significance(Z 0 ). Under all cuts, we estimate the final significance of the signal detection to be 7.47 in the semileptonic channel qq νγγ , as shown in Table 2. For the heavier Higgs boson with mass M H = 400 GeV, from the distributions in Fig. 6, we choose the following kinematical cuts for the pure leptonic channel, 120 GeV < M γγ < 130 GeV, ∆φ(γγ) < 2.5 , ∆R(γγ) < 2.5 ,

Analyses of Heavier Higgs Boson with 400 GeV and 600 GeV Masses
Comparing with the previous case of M H = 300 GeV, we find that the distributions ∆φ( ), ∆R( ), ∆φ(γγ), and ∆R(γγ) damp faster in the larger ∆φ and ∆R regions, as shown in Fig. 6. This is because the di-Higgs bosons Table 4: Signal and background cross sections of both pp → WW * γγ → ν νγγ and pp → WW * γγ → qq νγγ processes at the LHC (14 TeV) after each set of cuts. The signal significance(Z 0 ) is computed for the LHC (14 TeV) with an integrated luminosity of 3 ab −1 . We input the heavier Higgs mass M H = 600 GeV, and set the sample signal cross section σ(pp → H → hh → WW * γγ) = 1 fb. From the 3rd to 5th columns, we present the signals and backgrounds after imposing each set of cuts. In the pure leptonic mode, we impose the Final Cuts M T ( νν), M( ), M T ( ννγγ), ∆φ( ), ∆R( ), ∆φ(γγ), and ∆R(γγ). In the semi-leptonic mode, we add the Final Cuts P T (γ), M T (qq ν), ∆φ(γγ), and ∆R(γγ).  Table 3, where we set a sample signal cross section σ(pp → H → hh → WW * γγ) = 3 fb . In this case, we derive a signal significance(Z 0 ) = 4.05 after all the kinematical cuts. We also note from Fig. 4 We summarize cut efficiency of the final state qq νγγ for M H = 400 GeV in Table 3. We derive a significance Z 0 = 6.22 after all kinematical cuts.
Next, for the heavier Higgs H with mass M H = 600 GeV , the distributions of pure leptonic mode are shown in Fig. 6. From these, we set up the following kinematical cuts, 120 GeV < M γγ < 130 GeV, E / T > 25 GeV, The cut efficiency for M H = 600 GeV is summarized in Table 4. With these, we summarize the cut efficiency of qq νγγ final state for M H = 600 GeV in Table 4. Since the typical production cross section with M H = 600 GeV becomes significantly smaller over the parameter space, we take a sample input σ(pp → HX) × Br(H → hh → WW * γγ) = 1 fb , and consider an integrated luminosity of 3 ab −1 at the LHC (14 TeV). Hence, from Table 4, we can estimate the significance Z 0 = 7.76 and Z 0 = 9.29 for channels WW * γγ → ν νγγ and WW * γγ → qq νγγ , respectively. Besides, from Fig. 4(a)-(b) we see that for M H = 600 GeV, the cross section σ(pp → HX)×Br(H → hh → WW * γγ) in 2HDM-I can reach up to 3 fb, while this cross section in 2HDM-II is below about 0.2 fb. Thus, the significance for probing the 2HDM-II with M H = 600 GeV will be rescaled accordingly. In the following Sec. 4, we will give a general analysis of the significance by scanning the parameter space of 2HDM-I and 2HDM-II without assuming a sample cross section.

Probing 2HDM Parameter Space at the LHC
In this section, we study the probe of 2HDM parameter space by using the LHC Run-2 detection of the heavier Higgs state H 0 via pp(gg) → H → hh → WW * γγ (Sec. 3), as well as the current global fit for the lighter Higgs boson h 0 (125GeV) at the LHC Run-1. For the present analysis, we will convert the collider sensitivity (Sec. 3) into the constraints on the parameter space of 2HDM-I and 2HDM-II. As we showed in Fig. 3(c)-(d) and explained in the last paragraph of Sec. 2, the inclusive Higgs production cross section σ(pp → HX) is always dominated by the gluon fusion channel gg → H in the small tan β region, while other b-related channels are negligible. (For 2HDM-I, this feature actually holds for full range of tan β 1. ) Hence, for the present analysis, we will use Higgs production via gluon fusion pp(gg) → H → hh → WW * γγ (Sec. 3) to probe the 2HDM parameter space.
We combine the significance(Z 0 ) from both pure leptonic channel WW * γγ → ν¯ νγγ and semi-leptonic channel WW * γγ → qq νγγ at the LHC Run-2 with 300 fb −1 integrated luminosity. For this analysis, the relevant massparameters of the 2HDM are (M H , M A , M 12 ) . For demonstration, we will take the sample inputs, M H = 300, 400GeV and (M A , M 2 12 ) = (500GeV, −(180GeV) 2 ) . With these, we have two remaining parameters in the 2HDM: the mixing angle α and the VEV ratio tan β . In Fig. 9, we impose projected sensitivity of the LHC Run-2 by requiring significance(Z 0 ) > 5 . From this, we derive the red contours in the parameter space of cos(β − α) − tan β plane, for 2HDM-I [plots (a) and (c)] and for 2HDM-II [plots (b) and (d)]. The plots (a)-(b) correspond to M H = 300 GeV and plots (c)-(d) correspond to M H = 400 GeV . This means that the LHC Run-2 with an integrated luminosity L = 300 fb −1 can probe the red contour regions in each plot of Fig. 9 with a significance(Z 0 ) > 5 . It gives a discovery of the heavier Higgs boson H (with 300 GeV or 400 GeV mass) in the red regions of the 2HDM parameter space.
In Fig. 9, we further present the global fit for the lighter Higgs h (125GeV) by using existing ATLAS and CMS Run-1 data, where the 2 σ and 3 σ contours of the allowed parameter space are shown by the green and yellow shaded regions, respectively. As we checked, our LHC global fit of the 2HDM is consistent with those in the literature [36]. From this fit, we see that the parameter space favored by the current global fit is around the alignment limit of 2HDM with | cos(β−α)| 0.55 for 2HDM-I and | cos(β−α)| 0.15 for 2HDM-II. But, 2HDM-II still has an extra relatively narrow parameter region starting from tan β 2 . Fig. 9(a) has input M H = 300 GeV for 2HDM-I. In this plot, the Z 0 > 5 region overlaps a large portion of the parameter space favored by the current LHC global fit. But, in Fig. 9(b) for 2HDM-II, the situation is different because the overlap becomes smaller in the region cos(β − α) < 0 , and gets enlarged for cos(β − α) > 0 . For the case of M H = 400 GeV in Fig. 9(c), the probed parameter space of 2HDM-I has sizable reduction, especially for the region of cos(β−α) −0.05 , in comparison with Fig. 9(a) of M H = 300 GeV . This is because the signal rate decreases as H becomes heavier [cf. Fig. 4(a)]. On the other hand, for 2HDM-II, Fig. 9(d) shows that the Z 0 > 5 contours significantly shrink for M H = 400 GeV. This is because the signal rate for 2HDM-II drops more rapidly as Higgs mass rises to M H = 400 GeV in the small tan β region [cf. Fig. 4(b)]. In this case, we see that the LHC Run-2 with 300 fb −1 integrated luminosity has rather weak sensitivities to the parameter space (shown by red contours), and the red contours no longer overlap with the favored region by the current LHC global fit (yellow and green contours). We further analyze the probe from the upcoming High Luminosity LHC (HL-LHC) with 3 ab −1 integrated luminosity. We find that the HL-LHC can significantly extend the discovery reach of the parameter space of 2HDM-II, as demonstrated by the pink contour regions (Z 0 > 5) of Fig. 9(d).

Conclusion
After the LHC discovery of a light Higgs boson h 0 (125GeV) at Run-1, searching for new heavier Higgs state(s) has become a pressing task of the LHC Run-2. Such heavier Higgs state(s) exists in all extended Higgs sectors and can unambiguously point to new physics beyond the standard model (SM).
In this work, we systematically studied the heavier Higgs boson H 0 production with the new decay channel, pp → H → hh → WW * γγ , at the LHC Run-2. In section 2, we first analyzed the parameter space of the 2HDM type-I and type-II, including the Hhh cubic Higgs coupling (Fig. 1). We computed the decay branching fractions and production cross section of the heavier Higgs boson H at the LHC Run-2 over mass range M H = 250 − 600 GeV, as shown in Fig. 2-4. Then, in section 3, we analyzed both pure leptonic mode WW * → ν¯ ν and semi-leptonic mode WW * → qq ν . This channel has much cleaner backgrounds than the other process pp → H → hh → bbγγ . We computed signal and background events using MadGraph5(MadEvent). We applied Pythia to simulate hadronization of partons and adopted Delphes for detector simulations. We followed the ATLAS procedure for event selections and built kinematical cuts to efficiently suppress the SM backgrounds. We analyzed various kinematical distributions for pure leptonic and semi-leptonic decay channels in Fig. 6 and Figs. 7-8 for three sample inputs of Higgs mass M H = (300, 400, 600) GeV, respectively. In Table 2-4, we presented the signal and background rates of both channels under the kinematical cuts. In section 4, we combined the significance of pure leptonic and semi-leptonic channels, and analyzed the LHC Run-2 discovery reach of H as a probe of the parameter space in 2HDM-I and 2HDM-II (Fig. 9). For comparison, we further presented the current Higgs global fit of the LHC Run-1 data in the same plots. Finally, we note that it is hard to detect H with mass above 600 GeV at the LHC (14TeV) runs via di-Higgs production channel. We find it valuable to extend our present LHC study to the future high energy circular colliders pp(50−100TeV) [37], which are expected to further probe the heavier Higgs boson H with mass up to O(1−5) TeV range via pp → H → hh production channel.