LHC Search of New Higgs Boson via Resonant Di-Higgs Production with Decays into 4W

Searching for new Higgs particle beyond the observed light Higgs boson h(125GeV) will unambiguously point to new physics beyond the standard model. We study the resonant production of a CP-even heavy Higgs state $H^0$ in the di-Higgs channel via, $gg\to H^0\to h^0h^0\to WW^*WW^*$, at the LHC Run-2 and the high luminosity LHC (HL-LHC). We analyze two types of the $4W$ decay modes, one with the same-sign di-leptons ($4W\to\ell^\pm\nu\ell^\pm\nu 4q$) and the other with tri-leptons ($4W\to\ell^\pm\nu\ell^\mp\nu\ell^\pm\nu 2q$). We perform a full simulation for the signals and backgrounds, and estimate the discovery potential of the heavy Higgs state at the LHC Run-2 and the HL-LHC, in the context of generical two-Higgs-doublet models (2HDM). We determine the viable parameter space of the 2HDM as allowed by the theoretical constraints and the current experimental limits. We systematically analyze the allowed parameter space of the 2HDM which can be effectively probed by the heavy Higgs searches of the LHC, and further compare this with the viable parameter region under the current theoretical and experimental bounds.


Introduction
Since the LHC discovery of the light Higgs boson h 0 (125 GeV) in 2012 [1,2], both ATLAS and CMS collaborations have much improved the measurements on its mass and couplings which behave fairly standard-model-like [3][4][5]. But, so far its self-interactions have not yet been tested at the LHC. The cubic Higgs coupling of h 0 can be directly probed via the di-Higgs production at hadron colliders [6][7][8][9][10][11], though it would be much harder. At the LHC (14 TeV), the di-Higgs production cross section in the standard model (SM) is small. But, most extensions of the SM contain an enlarged Higgs sector and predict new cubic interaction H 0 h 0 h 0 between a heavier Higgs state H 0 and the light Higgs pair h 0 h 0 . For M H > 2M h , the di-Higgs cross section can be significantly enhanced via the resonant production pp → H 0 → h 0 h 0 , which simultaneously serves as an important discovery channel of new Higgs boson. Such an extended Higgs sector may include additional new singlets, doublets, or triplets under the SM gauge group SU(2) L ⊗ U(1) Y , or under an enlarged gauge group with extra SU (2). 1 Among these, the two-Higgs-doublet model (2HDM) [14] is a minimal extension by adding the second Higgs doublet to the SM. It is a necessary ingredient of the minimal supersymmetric SM (MSSM) [15] and its next-to-minimal extension (NMSSM) [16]. It is common to impose a discrete Z 2 symmetry on the 2HDM for preventing the tree-level flavor changing neutral currents. There are at least four kinds of model setup due to the different assignments of fermion Yukawa couplings with each Higgs JHEP06(2018)090 2 Heavy Higgs boson H 0 in 2HDM: decays and production In this section, we will first define the model setup of the 2HDM and its parameter space in section 2.1. Then, in section 2.2, we analyze the relevant theoretical constraints and the current direct/indirect experimental limits on the 2HDM parameter space. Finally, in section 2.3, we present the decays and production of the heavy Higgs boson H 0 at the LHC. With these, we set up three benchmarks for our LHC studies in the subsequent sections.

The model setup
The 2HDM [14] is a minimal extension of the SM Higgs sector. To avoid the tree-level flavor changing neutral currents, it is common to impose a discrete Z 2 symmetry on the Higgs sector, with the Higgs doublets H 1 and H 2 being Z 2 odd and even, respectively. Thus, the CP conserving Higgs potential under Z 2 can be written as, where all parameters are real [14,38,48,49]. The potential V respects the Z 2 symmetry except the mixing mass-term of M 2 12 which provides a soft breaking of Z 2 . This potential contains eight free parameters from the start, including three mass parameters (M 2 11 , M 2 22 , M 2 12 ) and five quartic couplings (λ 1 , λ 2 , λ 3 , λ 4 , λ 5 ). The vacuum is determined by the potential minimum, with the Higgs vacuum expectation values (VEVs), H j = (0, v j / √ 2). Both Higgs fields contribute to the electroweak symmetry breaking, with their VEVs obeying the condition v = (v 2 1 +v 2 2 ) 1/2 246 GeV. Defining v 1 = v cos β and v 2 = v sin β, we see that the VEV ratio is described by the parameter tanβ = v 2 /v 1 . The two Higgs doublets contain eight real components in total, The Higgs VEVs satisfy the extremal conditions, These are equivalent to the Higgs tadpole conditions h 0 1 = h 0 2 = 0, and determine the two VEVs (v 1 , v 2 ) or (v, tan β) as functions of the eight parameters in the Higgs potential (2.1).
With the three massless would-be Goldstone bosons (π ± j , π 0 j ) eaten by the weak gauge bosons (W ± , Z 0 ), the physical spectrum consists of five states: two CP-even neutral scalars (h 0 1 , h 0 2 ), one pseudoscalar A 0 , and a pair of charged scalars H ± . The CP-even sector -3 -

JHEP06(2018)090
involves a generic mass-mixing between (h 0 1 , h 0 2 ), and the mass eigenstates (h 0 , H 0 ) are given by the orthogonal rotation with mixing angle α, Given the current LHC data, it is most natural to identify lighter Higgs state h 0 as the observed Higgs boson of mass 125 GeV [1,2]. The heavier state H 0 is a brandnew Higgs boson beyond the SM, and can have sizable di-Higgs decays H → hh for the mass-range M H > 2M h . The current LHC measurements show that the Higgs boson h 0 (125 GeV) behaves rather SM-like. The favored 2HDM parameter space is then pushed to the region around alignment limit cos(β −α) ∼ 0. As mentioned above, the Higgs potential (2.1) contains eight parameters in total, including three mass parameters and five dimensionless self-couplings. We can reexpress the eight parameters in terms of four Higgs masses (M h , M H , M H ± , M A ), the combined VEV v = v 2 1 +v 2 2 , the VEV ratio tan β = v 2 /v 1 , the mixing angle α, and the mixing mass parameter M 2 12 . Imposing the experimental inputs v 246 GeV and M h 125 GeV, we note that the Higgs sector is described by six parameters in total: the VEV ratio tan β, the mixing angle α, heavy Higgs masses (M H , M A , M H ± ), and the mass-mixing parameter M 12 . Thus, we can express the five dimensionless Higgs couplings λ j as follows, which are consistent with [48]. Then, the masses M 2 11 and M 2 22 can be solved from eqs. (2.3) and (2.5), so they are not independent parameters.
For the later numerical analyses in this section and in section 4, we will consider the 2HDM parameter space in the following ranges, , and all the mass parameters are in the unit of GeV. Here we choose the value of M ± to be consistent with the bound of the weak radiative B-meson decays [50]. It was found that the B-decay constraint is quite weak for 2HDM-I, since it only requires M H ± > 450 GeV for tan β 1 and quickly drops below the LEP -4 -
The 2HDM type-I and type-II are defined according to their different assignments for the Yukawa sector under Z 2 symmetry. In the 2HDM type-I, all the SM fermions are defined as Z 2 even, thus only the Higgs doublet H 2 joins Yukawa interactions and generates all the fermion masses. For the 2HDM type-II, all the right-handed down-type fermions are assigned as Z 2 odd, while all other fermions are Z 2 even. Thus, the 2HDM-II has the Higgs doublets H 2 and H 1 couple to the up-type and down-type fermions, respectively. Under the Z 2 assignments, the H 0 Yukawa couplings for 2HDM-I and 2HDM-II can be expressed in the form,

JHEP06(2018)090
violate the perturbative unitarity bounds for large tanβ. We will discuss this in more detail later when we show the results of the parameter scan in figure 2. The existing electroweak precision data will constrain the one-loop contributions induced by the Higgs-gauge couplings via oblique corrections [66]. Around the alignment limit, we can expand the 2HDM contributions to the oblique parameters [67][68][69] as follows, where F (x, y) = 1 2 (x+y)− xy x−y ln x y , and the function B 22 is given by where ∆(x, y) = 2(x + y) − (x − y) 2 − 1. The leading order contributions to the oblique corrections (2.9) only involve the masses of new Higgs bosons (H 0 , A 0 , H ± ). This is because the couplings of trilinear vertices involving two new Higgs bosons (W ± -H ∓ -H 0 , W ± -H ∓ -A 0 , Z 0 -H 0 -A 0 , and Z 0 -H + -H − ) either contain the factor sin(β − α) or have no (α, β)-dependence, while the cubic vertices with only one new Higgs boson (W ± -H ∓ -h 0 , Z 0 -A 0 -h 0 , and H 0 -V -V ) are suppressed by cos(β−α). Besides, the other cubic vertices h 0 -V -V (V = W ± , Z 0 ) have couplings proportional to sin(β −α), and lead to the suppression factor sin 2 (β−α) − 1 = − cos 2 (β−α) after subtracting the corresponding SM contributions.
In figure 1(a) we present the 2HDM predictions of (S, T ) by scanning the 2HDM parameter space, where the (red, green, blue) dotted regions correspond to H 0 mass M H = (300, 400, 500) GeV, respectively. As a comparison, we also show the 95% C.L. contour from the precision constraints [70] in the same plot (with U = 0). In figure 1(b), we present the U parameter prediction of the 2HDM over the mass-range M H = (0.3 − 0.8)TeV. We find that in the 2HDM, the oblique contribution to T is much larger than S, and the S and U parameters are fairly small in the relevant parameter region, namely, |T | |S| ∼ |U | < 2 × 10 −2 . Hence, figure 1 shows that the nontrivial constraint mainly comes from the T parameter. Since the current electroweak precision data constrain the T parameter to be quite small, especially for small |S| 2 × 10 −2 as restricted by the ellipse contour  in figure 1(a), this requires the mass of H ± to be fairly degenerate with that of H 0 or A 0 . (For the numerical analyses in figure 1, we have used the exact one-loop formulas for the oblique parameters [69]. We also compared our results with ref. [69] and the 2HDMC code [71] for consistency checks.) Next, we further derive the existing constraints on the 2HDM parameter space by making a global fit for the LHC Run-1 and Run-2 measurements on the light Higgs boson h (125 GeV). Since this Higgs state h (125 GeV) is fairly SM-like, the new physics corrections to h couplings are tightly constrained and other new Higgs states need to be significantly heavier. Hence, we may regard the h 0 global fit bounds as indirect constraints on the new Higgs states, in contrast to the constraints from their direct searches.
In our analysis, we perform the Higgs fit by minimizing χ 2 with the inputs of signal strengths µ i from the LHC experimental fits (including their errors and correlations). For the Run-1 Higgs data of the LHC (7+8 TeV), we have used the combined analysis of ATLAS and CMS [5]. For the current Higgs data from the Run-2 of the LHC ( We make use of these Run-1 and Run-2 Higgs data from ATLAS and CMS collaborations, and perform a global fit to derive the current LHC constraints on the 2HDM parameter space. We present the 2σ contours (red curves) and 3σ contours (blue curves) in the cos(β − α) − tan β plane for 2HDM-I in figure 2(a) and for 2HDM-II in figure 2(b). We further incorporate the theoretical requirements from the Higgs stability and perturbative unitarity, and present the combined constraints on the allowed parameter regions marked by red dots (in box shape) at the 2σ level and by blue dots (in circle shape) at the 3σ level.  In figure 2, we are making a uniform parameter scan according to eq.(2.6), which shows that the likelihood for the red and blue dots becomes smaller when tanβ increases. 2 We see that without special fine-tuning of the parameter space, the current fit favors the relatively low tanβ region, i.e., tanβ 9 for 2HDM-I and tanβ 7 for 2HDM-II. This is largely due to the unitarity bounds (2.8a)-(2.8b) in which the Higgs coupling λ 1 is enhanced by 1/cos 2 β ∼ tan 2 β [cf. eq.(2.5a)]. The other Higgs couplings (2.5c)-(2.5e) are also enhanced by the factor sin2α/sin2β ∝ tanβ and 2M 2 12 /sin2β ∼ M 2 12 tanβ for large tanβ. If we select a special parameter space obeying the condition we find that eq.(2.5a) reduces to λ 1 = (sin 2 α/cos 2 β)(M 2 h /v 2 ), and the enhancement factor 1/cos 2 β ∼ tan 2 β can be removed by taking the alignment limit cos(β − α) → 0. The alignment limit gives α = β − π 2 , under the convenition [84], β ∈ [0, π 2 ] and β − α ∈ [0, π]. In this limit, the above condition (2.11) becomes M 2 12 = 1 2 sin 2β M 2 H , and the Higgs couplings (2.5a)-(2.5e) reduces to (2.12d)

JHEP06(2018)090
We see that the resultant Higgs couplings above do not increase with tanβ and could more easily satisfy the unitarity bounds (2.8) only if the squared masses of the heavy Higgs bosons are all nearly degenerate M 2 For illustration, we make a separate scan on the special parameter region under the condition (2.11) [in addition to the default condition (2.6)] and for tanβ > 5. This is represented by the green dots in figure 2 and indeed reaches larger tanβ regions for both 2HDM-I and 2HDM-II, as expected. But we have to keep in mind that the condition (2.11) only represents a very small region of the generic parameter space defined in eq.(2.6), so it has low likelihood and is hard to be reached by the conventional uniform parameter scan. For consistency check, we have also made comparison with refs. [38,49] for the global Higgs fit. Figure 2 shows that the current LHC global fit of h (125 GeV) shifts the viable parameter space somewhat towards cos(β − α) < 0 region for 2HDM-I, while it pushes the allowed parameter range significantly to cos(β − α) > 0 side for 2HDM-II. The reasons are the following. We note that the current LHC data give tight constraints on the signal strengths µ(ggF+tth) and µ (VBF+V h). From the combined Run-1 data [5], we see that h 0 → W W * channel favors µ (VBF+V h) > 1 and h 0 → bb channel favors µ (VBF+V h) < 1, which can be both explained by a reduced hbb coupling (relative to its SM value), since the reduced hbb coupling can decrease Br(h 0 → bb) and enhance branching ratios of all other final states. Also, h 0 → ZZ channel prefers µ(ggF+tth) > 1 and thus an enhanced htt coupling. For the LHC Run-2, we note that ATLAS data [72] significantly favor µ (VBF)> 1 in both h 0 → γγ, ZZ * channels, while their combination has little effect on µ(ggF). The V h production at ATLAS Run-2 also prefers µ(V h) < 1 via h 0 → bb channel [76]. Thus, these features can be explained by a reduced hbb coupling. Besides, the CMS Run-2 data [77] mildly favor µ (VBF+V h) < 1 via h 0 → ZZ * and thus an enhanced hbb coupling, while h 0 → γγ channel at the CMS Run-2 slightly prefers µ(ggF+tth) > 1 [78] and so an enhanced htt coupling. Finally, inspecting the hbb and htt couplings in table 1 and expanding them around the alignment limit (with β − α ≡ π 2 − δ and δ as a small deviation), we find that up to the first order of δ, the hbb coupling equals 1+δcotβ in 2HDM-I and 1− δ tanβ in 2HDM-II, while the htt coupling equals 1 + δcotβ in both 2HDM-I,II. Hence, a reduced hbb coupling requires δ < 0 and thus cos(β − α) < 0 in 2HDM-I, pushing the viable parameter space towards the left-hand-side in figure 2(a); while for 2HDM-II, this requires δ > 0 and thus cos(β−α) > 0, shifting the parameter space towards the right-hand-side in figure 2(b). Also, a mildly enhanced htt coupling requires δ > 0 for both 2HDM-I,II, so its effect will partially cancel that of the htt coupling for 2HDM-I, and add together with the htt effect for 2HDM-II. This explains that the 2HDM-II has larger asymmetry over Finally, for this study, we will consider the upper bounds from the existing LHC Run-1 and Run-2 searches on a heavier neutral Higgs state H 0 with decays in various channels. These will put additional constraints on the 2HDM parameter space through various H 0 couplings. The LHC Run-1 searches include H 0 → hh → bbγγ from CMS [19], the combined searches of H 0 → hh → bbγγ, bbbb, bbτ τ , γγW W * from ATLAS [20], H 0 → ZZ from ATLAS [85], and the combined searches of H 0 → W + W − /ZZ from CMS [86]. The LHC Run-2 searches include H 0 → W W from ATLAS [87], H 0 → ZZ from ATLAS [  H 0 → τ τ from ATLAS [90] and from CMS [91], and H 0 → hh → bbτ τ from CMS [26]. With these, we summarize in table 2 the current upper limits (95% C.L.) from the LHC Run-1 and Run-2 direct searches on the production cross sections of the heavier Higgs boson H 0 in various decay channels, where the numbers do not contain the decay branching fractions of the final states V V , hh, and τ τ .
For the heavy Higgs state H 0 , we analyze its decay branching fractions in figure 3. We present the branching fractions over the mass range M H = (300 − 800) GeV, for tanβ = 2  Then, in figure 4, we present the current experimental constraints on the 2HDM parameter space in the plane of M H − cos(β−α). The parameter region with blue dots (circle shape) satisfy the theoretical conditions, the electroweak precision limits (2σ), and the LHC bounds (2σ) from the Higgs global fit of h (125 GeV) data. The red dots (square shape) present the parameter region obeying the existing LHC direct search limits (2σ) on the heavier Higgs boson H 0 in combination with the theoretical constraints. The electroweak precision tests mainly bound the oblique parameter T as shown in figure 1, and prefer the masses M H ± and M A to be fairly degenerate for M H 500 GeV. The present LHC global Higgs fit prefers h (125 GeV) to be quite SM-like, and favors the 2HDM parameter space around the alignment limit (cf. figure 2). Figure 4 shows that the allowed region of cos(β−α) (with blue dots) in 2HDM-I is more shifted to cos(β−α) < 0 as in plot-(a), while the region with blue dots in 2HDM-II is largely excluded on the cos(β −α) < 0 side, as in plot-(b). These features are consistent with figure 2. The current LHC direct search limits on the heavier Higgs state H 0 are reflected by the red dotted regions in figure 4. They are comparable to the bounds imposed by the LHC h (125 GeV) global fit (combined with the electroweak precision limits) for 2HDM-I, but they are significantly weaker for the case of 2HDM-II.

H 0 production in di-Higgs channel and benchmarks
For the present study, we focus on the productions and decays of the heavy Higgs state H 0 . We have summarized the H 0 Yukawa couplings for 2HDM-I and 2HDM-II in table 1. The gauge couplings of H 0 take the form , which differs from the SM Higgs-gauge coupling by a factor cos(β − α). With these, we determine the ggH vertex by rescaling SM contributions inside the loop accordingly. The ratios of -11 -JHEP06(2018)090 the decay widths with respect to the SM results are given as follows, Around the alignment limit the decay width Γ(H → V V ) is suppressed by cos 2 (β−α), while for down-type fermions the partial width Γ(H → f f ) is enhanced by a factor tan 2 β.
The main production mechanism of the neutral Higgs boson H 0 at the LHC is the gluon fusion production gg → H 0 . In the 2HDM-I, all Yukawa couplings rescale by a common factor (sin α/ sin β) with respect to the SM values. So the gluon fusion production is still dominated by the top-loop, and the cross section is rescaled by a factor (sin α/sin β) 2 . In the 2HDM-II, the up-type and down-type Yukawa couplings have different rescalings from their SM values, where the rescaling of down-type Yukawa couplings is proportional to 1/ cos β ∝ tan β . For tan β 1, the bottom-loop in the gluon fusion process may have visible contribution, while it becomes negligible for tan β = O(1). Figure 2(b) shows that the current constraints strongly favor tan β = O(1) for the 2HDM-II, and thus the bottomloop contribution is negligible [14]. We take the four-flavor scheme in the present study, and the relevant b-related production process is the bottom-pair associated production gg → Hbb [15,92,93], which is also negligible for the 2HDM-I and for the 2HDM-II [with small tan β = O(1)].
We can deduce the coupling of the cubic scalar vertex Hhh from the Higgs potential, where we expand the formula around the alignment limit in the second line. For the mass-range M H > 2M h , the tree-level decay width of H→ hh is The di-Higgs decay width is also suppressed by cos 2 (β −α) in the alignment limit, but it may receive enhancement from other mass-parameters M 2 12 and M 2 H in the Higgs potential.
-12 - To set viable benchmarks for searching the new Higgs boson H 0 via resonant di-Higgs production in hh → W W * W W * channel, we will implement the theory constraints and the current experimental bounds on the 2HDM parameter space, as we have analyzed in section 2.2. The requirements of vacuum stability and perturbative unitarity strongly favor the small tanβ region [except the alignment limit cos(β − α) ∼ 0] [31,38,49]. For small tanβ, the bottom-loop contribution in the gluon fusion production can be safely ignored [14]. Besides, the bottom-pair associated production is subdominant, especially for the experimental searches aiming at gluon fusion production without making extra btagging. According to the analysis in section 3, we generate events for pp → H 0 bb in the four-flavor scheme. With the selection cuts aiming at gluon fusion production, in particular the b-veto which helps to suppress the top-related backgrounds, we find H 0 bb contribution unimportant for the current study. Hence, in the following we will focus on the gluon fusion production, . We note that the existing constraints in figure 4 will set upper bounds on the resonant H 0 production cross sections at the on-going LHC Run-2 and the HL-LHC. In figure 5, we plot the red dots (square shape) to show the viable parameter region allowed by the LHC limits (2σ) of the existing H 0 direct searches combined with the theoretical requirements. For comparison, the blue dots (circle shape) represent the parameter space obeying the theoretical requirements, the indirect electroweak precision limits (2σ), and the -13 -

JHEP06(2018)090
LHC bounds (2σ) from the global fit of h (125 GeV). Inspecting figure 5, we see that the allowed H 0 production cross sections in the 2HDM-I and 2HDM-II are not much different over the wide mass range of M H = (300 − 1000) GeV. But, the distribution of the blue dots for 2HDM-II [plot-(b)] is relatively sparser than that for 2HDM-I [plot-(a)], due to the stronger constraints on the 2HDM-II by the current LHC global fit of h (125 GeV) (cf. figure 4). As a side remark, for the neutral pseudo-scalar Higgs boson A 0 , it has no A-h-h vertex and A-V -V vertex (V = W, Z) at tree level. Hence its dominant decay channel is [15,94], where the final state fully differs from that of our diHiggs channel H → hh → 4W and thus does not affect our current LHC analysis.
Based upon our analyses of the existing indirect and direct experimental bounds on the 2HDM (combined with theoretical constraints), we will systematically study the direct probe of the heavy Higgs boson via gg → H → hh → W W * W W * channel at the LHC (14 TeV) in the following section 3. For this, we set up three benchmark scenarios for the mass M H and the cross section σ(gg → H)×Br(H→ hh → W W * W W * ) as follows, (M H , σ×Br) = (300GeV, 60fb), (400GeV, 40fb), (500GeV, 12fb), (2.18) which will be denoted by (H300, H400, H500) for short.
For an illustration, we further present three explicit parameter samples to realize the benchmarks (2.18). It is consistent with all the theoretical and experimental constraints discussed above. We show this sample in table 3

Analyzing signals and backgrounds at the LHC
In this section, we perform systematical Monte Carlo analysis for the resonant neutral Higgs H 0 signal gg → H 0 → h 0 h 0 → W W * W W * and its main backgrounds at the LHC (14 TeV). We set up the signal process model by FeynRules [95] with ggH 0 and H 0 h 0 h 0 vertices. We generate the events by MadGraph5 package [96] at parton level, and then process them by Pythia [97] for hadronization and parton shower. Finally, we use Delphes 3 [98] for detector simulations. Here, for the resonant diHiggs production, we include the new vertices  precise form factor from the top-loop, which only depends on the masses M H and M t . We include the same K-factor for the NLO QCD corrections as in the SM-type Higgs production gg → H 0 [99]. As consistency checks, we have used the SusHI package [100] to recompute the resonant production cross section σ(gg → H 0 ) and get full agreement. We also note that the cross section of the resonant on-shell H 0 production (followed by the cascade decay H 0 → h 0 h 0 ) overwhelms the nonresonant SM contribution [101]. We study two major decay channels of the final state W W * W W * : (i). ± ν ± ν 4q with same-sign di-leptons (SS2L); (ii). ± ν ∓ ν ± ν 2q with tri-leptons (3L). For W boson, its branching fractions of leptonic decays W → eν, µν, τ ν are 10.8%, 10.6%, and 11.3%, respectively, while its hadronic decay W→ qq has branching fraction 67.6% [102]. For this analysis, we include the detected e and µ from τ decays as well. These two decay channels of W W * W W * have branching fractions about 9.6% and 9.2%, respectively. Although they are quite small (less than 10%), we note that requiring the detection of the same-sign dileptons or the tri-leptons in the final state can significantly reduce the QCD backgrounds and enhance the signal sensitivity.

Final state identification
To analyze the signal sensitivity, we apply the ATLAS procedure to identify the final states in both SS2L and 3L decay channels. Jets, leptons, and transverse missing energy are selected by the following cuts, Electrons with 1.37 < |η( )| < 1.52 are rejected in order to remove the transition region of electromagnetic calorimeter of ATLAS. After the trigger, we require the leading lepton passing the trigger requirement p T ( ) > 25 GeV. The b-taging algorithm based on p T of jet is implemented in Delphes [103]. The reconstructed objects in the final state have to be well separated spatially to prevent the potential double-counting. We implement the following criteria [104]: (i). Any -15 -

Analysis of same-sign di-lepton decay channel
With the identification of the final state particles as in section 3.1, our analysis of the samesign di-leptons (SS2L) channel further requires the sub-leading lepton obeying p T ( ) > 25 GeV to reduce the fake backgrounds (as will be described in the following) and n = 2 (same-sign), n j 3 . (3. 2) The jets arising from the off-shell W boson decays could be soft and the requirement of n j 3 provides an optimal significance. The above defines the basic event selection for the SS2L channel.
The main prompt backgrounds that contribute to the same-sign di-leptons include , ZZ (with Z → ) and W Z (with W → ν and Z → ). The background processes with a pair of top quarks (namely, ttW , ttZ, and tth) can be efficiently rejected by b-veto. With the basic event selection and b-veto, the ttZ background becomes negligible. The diboson backgrounds, W Z, ZZ, and Zh, can be suppressed by requiring exactly two same-sign leptons as in eq. (3.2). Hence, we can safely ignore ZZ and Zh backgrounds given their small cross sections, and only include W Z channel for the background estimate.
For the SS2L channel, backgrounds with fake leptons from jets or charge misidentifications (QmisID) can also be significant. Jet faking leptons mainly come from W +jets final state and semi-leptonic mode of tt final state (which both have large cross sections). For the samples of fake electrons, we assign a weight to each event in the following way. (i). We generate W +jets background (including two or more jets) and the tt background (in the semi-leptonic decay mode). (ii). For each selected event of lepton+jets, we loop over all possible jets that could fake an electron with a probability P = 0.0048×exp[−0.035×p T (GeV)] as a function of jet p T [105]. (iii). We sum over all fake rates and divide it by two to account for the same sign fakes with the selected leptons. (iv). We randomly choose one jet to be the fake electron according to the fraction of fake rates, and rescale the jet's energy to its 40% as that of the fake electron [105]. Since fake electrons are usually soft, we find that the selection cut p T (e) > 25 GeV helps to significantly suppress these backgrounds and makes them comparable to other prompt backgrounds. With the upgrade Level-0,1 Muon Trigger for ATLAS detector [106], the contribution of fake muons with p T (µ) > 25 GeV is small and can also be safely ignored. The QmisID mainly comes from the Z+jets and the pure leptonic mode of tt, with one charge misidentified lepton. Pseudo-events are generated in a way similar to that for the fake electrons, with a weight 0.0026 assigned to each event [106]. These backgrounds can be significantly suppressed by Z-veto, i.e., |M ( ) − M (Z)| > 10 GeV. With these, we find that the contribution due to QmisID is negligible.   Table 4. σ × Br for the signal process (with three benchmarks) and for the major backgrounds in the SS2L decay channel before and after pre-selections (PreS), where the last two categories denote backgrounds with fake leptons and charge-misidentifications. For each background (with a cited reference), the cross section includes the QCD corrections (K-factor), while the rest of cross sections are computed by MadGraph5 at the leading order.
After all pre-selection cuts, including the basic event selection, b-veto, and Z-veto, we obtain σ ×Br as shown in the last column of table 4. We will make further optimization to increase the significance. In figure 6, we present the distributions of kinematic variables after all pre-selections, each of which has its own advantage to help the discrimination of signals from backgrounds. To maximize the background rejection, we train a boosted decision tree (BDT) discriminant, implemented in a TMVA package, using all the input variables [113] for each benchmark. By scanning the median significance Z [114] (cf. eq. (3.5)) as a function of the BDT response cuts, we select the BDT-optimized working point before the statistical fluctuations becoming important. As an illustration, figure 7(a) shows the distribution of the output BDT response with an integrated luminosity L = 1000 fb −1 for the case of M H = 400 GeV. In figure 7(b), we describe the median significance Z as a function of We input an integrated luminosity of (300, 1000, 3000)fb −1 for M H = (300, 400, 500) GeV.  Table 5. BDT optimization in the same-sign di-lepton (SS2L) channel for the three benchmarks with different integrated luminosity L (in the unit of fb −1 ). We present, after the optimal BDT cuts, the selected event numbers for the signals, backgrounds, and median significance Z.  Table 6. σ × Br for the signal process (with three benchmarks) and the major backgrounds in the 3L decay channel before and after the pre-selections (PreS), where the last category denotes backgrounds with fake leptons.

Analysis of tri-lepton decay channel
The analysis of the tri-lepton (3L) channel is similar to that of the SS2L channel. With the final state identification in section 3.1, we require the sub-leading electron to satisfy p T (e) >20 GeV, and the sub-leading muon to obey p T (µ) > 25 GeV. These will help to significantly suppress the fake electrons and fake muons. Furthermore, we implement a stronger cut / E T > 20 GeV, since the third neutrino in the 3L channel can contribute to the reconstructed transverse missing energy with less back-to-back cancellation than the case of the SS2L channel. The events need to further satisfy, n = 3 (total charge = ±1), n j 2 . (3. 3) The requirement on the total charge can enhance the Z-veto efficiency. Among the three leptons, 0 denotes the lepton with the opposite charge to the others, 1 for the one closest to 0 , and 2 for the remaining one. Thus, we define the Z-veto by |M ( 0 1 )−M Z | > 20 GeV and |M ( 0 2 )−M Z | > 10 GeV. The background estimation for the 3L channel also has several differences. Besides those discussed in the analysis of SS2L channel, the prompt backgrounds include ZZ final state (which might not be negligible here) and tri-boson processes W W W and W W Z. For the fake backgrounds, only jet-faked leptons from the final state Z+jets and the final state tt in pure leptonic mode need to be considered due to their large cross sections. We generate samples with fake electrons in the same way as described for the analysis of SS2L channel. The fake backgrounds from Z+jets can be largely suppressed by Z-veto. We show σ×Br after the pre-selections in table 6.
-20 -  Table 7. BDT optimization in the tri-lepton (3L) channel for the three benchmarks with different integrated luminosity L (in the unit of fb −1 ). We present, after the optimal BDT cuts, the selected event numbers for the signals and backgrounds (in which the Z+jets background is negligible and not listed here), as well as the median significance Z.

JHEP06(2018)090
In figure 8, we present the distributions of useful kinematical variables after the preselection. Using these we further implement the BDT optimizations in table 7. By definition, 0 and 1 come from the decays of one Higgs boson h, while 2 and the two closest jets arise from the decays of the other Higgs boson h. The invariant masses, M ( 0 1 ) and M ( 2 jj), are then controlled by the light Higgs mass M h . The ∆R distances, ∆R( 0 , 1 ) and ∆R( 2 , jj), tend to be smaller for increasing M H and are useful for the optimization. In summary, ∆R( 0 , 1 ) and ∆R( 2 , jj) have better powers of background suppression for the case of larger M H , while the rest are important for the case of lower M H .
Using the same procedure, we can implement the BDT optimization in the analysis of 3L channel for the three benchmarks. Our results are summarized in table 7. The dominant prompt background W Z gets efficiently suppressed by the valuable variables constructed for 0 and 1 . Impressively, the fake leptons backgrounds become almost negligible after the optimal BDT cut. With these, we have achieved higher sensitivity for the 3L channel since it has less backgrounds than the SS2L channel.
In passing, we also compare the multi-variable-analysis via the BDT approach with the traditional cut-based method, and show how the BDT approach can improve the sensitivities. For this comparison, we will use the same set of optimization variables in the two analyses. For the decay channels 2LSS and 3L, we list the set of variables as follows,  As an illustration, we optimize the significance with the two most sensitive variables, p sum T and ∆R( 1 , j) for the 2LSS channel, and p sum T and ∆R( 0 , 1 ) for the 3L channel. We summarize the results in table 8 for the benchmarks (H300, H400, H500), where we apply the optimal cuts p sum T < (300, 500, 600) GeV and ∆R( 1 , j) < (2.5, 1.5, 1.0) to the 2LSS channel, and p sum T < (300, 500, 600) GeV and ∆R( 0 , 1 ) < (2.0, 1.5, 1.2) to the 3L channel. We see that the BDT optimization gains about (16 − 64)% increase of significance over the cut-based method.   Table 8. Comparison of the significance S/ √ S +B for the cut-based analysis and the BDT approach. In the 4th and 8th columns, the quantity R (BDT/Cut) equals the ratio of the significance of the BDT approach over that of the cut-based analysis.

Combination of the SS2L and 3L decay channels
According to the above systematical analyses of the signals and backgrounds in the SS2L and 3L decay channels for the three Higgs benchmarks, we will study the combined sensitivity of both the SS2L and 3L channels for detecting the new Higgs boson H 0 in this subsection. To handle the relatively small event number, we will use the median significance Z [114]. For estimating the sensitivity of future experiments, we use the following median significance for the discovery reach under the background-only hypothesis [114], while for the exclusion, we use the formula under the background with signal hypothesis [114], For the case of B S, we can expand the formulas (3.5) and (3.6) in terms of S/B or S/(S +B), and find that they both reduce to the form of Z = S/ √ B or Z = S/ √ S +B , as expected.
Using our results from the SS2L and 3L decay channels as presented in table 5 and  table 7, we combine the significance Z from table 5 and table 7, which is summarized in table 9 for the three Higgs benchmarks in eq.(2.18). In the fourth row of table 9, we derive the combined significance Z comb of the benchmarks (H300, H400, H500) under sample inputs of the corresponding integrated luminosity L = (300, 1000, 3000) fb −1 , respectively. In the second row of this table, we also present the required minimal integrated luminosity L 5σ min to reach the significance Z comb = 5 for each benchmark, by combining the LHC searches in both SS2L and 3L channels. We see that given an integrated luminosity of 446 fb −1 , the LHC (14 TeV) can reach a discovery of the heavier Higgs state H 0 with mass up to 400 GeV via the di-Higgs channel H 0 → h 0 h 0 → 4W .
In passing, we have also compared the significance of our diHiggs decay channel hh → 4W with another channel hh → bbW W * in [30,40] which studied the SM extension with a real singlet scalar S. After proper rescaling, we find that for M H < 400 GeV, the 4W channel has better sensitivity than the bbW W channel [30,40] for detecting the heavy Higgs H 0 , while for M H 400 GeV, the bbW W channel study in [40] has higher sensitivity.

Probing the parameter space of 2HDM
In this section, we analyze the probe of the 2HDM parameter space by searching the heavy Higgs resonant production in the di-Higgs channel gg → H 0 → hh → W W * W W * at the LHC (14 TeV). For this, we shall combine the analyses of both SS2L and 3L channels in section 3. We further take into account both the theoretical constraints and the current experimental bounds as studied in section 2.
In figure 9, , respectively. The blue dots (square shape) present the allowed parameter region satisfying the theoretical constraints and the indirect experimental bounds (including the electroweak precision limits and the LHC global fit of the SM-like Higgs boson h (125 GeV), as we discussed in section 2). The red dots (circule shape) represent the parameter space which can be probed by the direct heavy Higgs searches at the LHC with a significance Z 2, including the existing heavy Higgs search bounds and our study of gg → H → hh → 4W searches (combined with the theoretical constraints). For the heavy Higgs searches via 4W channel, we derive the expected sensitivity from the LHC Run-2 with an integrated luminosity L = 300 fb −1 .
From figures 9(a) and 9(b), we see that for both 2HDM-I and 2HDM-II with M H = 300 GeV, the direct searches of H 0 at the LHC Run-1 and Run-2 can substantially probe the parameter space towards the alignment limit (represented by the regions with red dots), and largely cover the viable parameter space allowed by the current indirect searches (shown by the regions with blue dots). We note that the red dots can cover sizable regions around the alignment limit cos(β −α) = 0. This is because the existing direct searches of H 0 at the LHC Run-2 via H → τ τ channel [90,91]  We also see that the direct heavy Higgs searches can probe the viable parameter region up to tan β 5 for 2HDM-I and tan β 3 for 2HDM-II. Given the lower capability of the LHC Run-2 for detecting H 0 in the mass range M H 400 GeV, we expect that more sensitive direct probes should be achieved at the HL-LHC.
Next, we further extend our above analysis to the HL-LHC with an integrated luminosity L = 3000 fb −1 . We present this in figure 10, including all three sample inputs of the JHEP06(2018)090 Here the blue dots present the viable parameter space allowed by the current indirect constraints, which are the same as in figure 9. For the bounds from LHC direct searches of the heavy Higgs boson H 0 , we identify the parameter space with red dots in the same way as in figure 9, except that we include the HL-LHC probe by the resonant di-Higgs production in the 4W channel with L = 3000 fb −1 . In the relatively low mass range such as M H = 300 GeV, we see that the improvements are  . The parameter region with red dots (circle shape) can be probed by the LHC direct searches of the heavier Higgs boson H, including the existing heavy Higgs search bounds and our study of gg → H→ hh → 4W searches at the HL-LHC (14 TeV) with L = 3000 fb −1 (combined with the theoretical constraints) for all cases. As in figure 9, the blue dots (square shape) satisfy the theoretical constraints, the electroweak precision limits and the Higgs global fit of h(125GeV). All these bounds are shown for significance Z 2.
-26 -JHEP06(2018)090 been significantly expanded. For the case of M H = 500 GeV, figures 10(e)-(f) demonstrate that the projected sensitivities (regions with red dots) via the direct heavy Higgs searches, in comparison with the current indirect bounds (regions with blue dots). Since the current direct heavy Higgs search limits are rather weak for M H = 500 GeV, in figures 10(e)-(f) the regions with red dots are mainly probed by the direct search process gg → H → hh → 4W at the HL-LHC. We see that the projected sensitivities in the plots (e)-(f) are somewhat weaker than the case of M H = 400 GeV in the plots (c)-(d), but they are still comparable. This shows the power of the HL-LHC runs for probing the 2HDM parameter space with higher masses of the heavy Higgs boson H 0 . It is clear that the HL-LHC (L = 3000 fb −1 ) has significantly increased the sensitivity to probing the 2HDM parameter space via the H → hh → 4W channel.
We may compare the current analysis of gg → H → hh → W W * W W * with our previous study for gg → H → hh → γγW W * (including both the pure leptonic and semi-leptonic decays of the W W * final state) [31]. We find that our combined sensitivity in the W W * W W * channel is comparable to that of the γγW W * channel. Our present BDT optimization analysis in section 3 does help to improve the sensitivity as shown in table 8. Even though the W W * W W * channel is not the most sensitive search channel for detecting of the heavier Higgs boson H 0 via resonant dihiggs production, it is important to study all possible di-Higgs production channels which will allow a combined analysis of the H 0 discovery. In passing, we note that for other extended Higgs sectors, there could be a new singlet scalar S 0 beyond the SM Higgs doublet [40,116] or beyond the two Higgs doublets [47]. The scalar particle S 0 can mix with the Higgs states (h 0 , H 0 ), so it will couple to quarks tt and bb, and to gauge bosons W W and ZZ. Thus, the present analysis can be also generalized to the resonant channels gg → S → hh → 4W and gg → H → SS, Sh → 4W .

Conclusions
It is a pressing task for the LHC Run-2 and the upcoming HL-LHC to search for new Higgs state(s) beyond the light Higgs boson h (125 GeV), which generally exists in all extended Higgs sectors and would point to new physics beyond the standard model (SM) without ambiguity. The resonant di-Higgs production is an important channel to search for the heavy neutral Higgs state H 0 , which also directly probes the cubic Higgs interaction Hhh.
In section 2, we analyzed the viable parameter space for the 2HDM Type-I and Type-II, which satisfies the theoretical constraints and the current experimental limits. These are described in figures 1-4. Then, we presented in figure 5 the resonant di-Higgs production with σ(gg → H)×Br(H→ hh) as a function of the Higgs mass M H at the LHC (14 TeV), under the theoretical constraints and the current indirect and direct experimental limits. With these, we further set up three Higgs benchmark scenarios as in eq.(2.18) for the subsequent LHC analyses.
In section 3, we performed systematical Monte Carlo analysis for the resonant di-Higgs production via gg → H 0 → hh → W W * W W * channel at the LHC (14 TeV) by using Delphes 3 fast detector simulations. We studied the decay channels with the same-sign di-leptons (SS2L) final state in section 3.2 and with the three leptons (3L) final state in -27 -JHEP06(2018)090 section 3.3, where the QCD backgrounds can be efficiently suppressed. The top quark and Z boson induced backgrounds can be suppressed by b-veto and Z-veto, respectively. We analyzed the event distributions of different kinematical variables in figure 6 for the SS2L channel and in figure 8 for the 3L channel. We further optimized the signal significance by using the BDT method for each Higgs benchmark as presented in figure 7 and tables 4, 6. We found that the 3L channel has a higher significance than the SS2L channel. In section 3.4, we derived the combined significance for both SS2L and 3L channels, as summarized in table 9 for the three Higgs benchmarks. This table shows that to discover the new Higgs state H 0 at the LHC (14 TeV) with a 5σ significance, the required minimal integrated luminosity is L 5σ min = (222, 446, 2954) fb −1 for M H = (300, 400, 500) GeV, respectively. In section 4, we systematically analyzed the LHC probe of the 2HDM parameter space by using the combined searches of both SS2L and 3L decay channels of the 4W final state. In figure 9, we presented the parameter space (red dots) which can be probed by the LHC direct searches of the heavy Higgs boson H 0 , including the existing H 0 search bounds and our study of gg → H 0 → hh → 4W searches at the LHC Run-2 with L = 300 fb −1 (combined with the theoretical constraints). In figure 10, we further demonstrated that with an integrated luminosity of 3000 fb −1 at the HL-LHC, the probed parameter space of the 2HDM-I and 2HDM-II can be significantly expanded towards the alignment region, especially for the higher Higgs mass range M H 400 GeV. We find that considerable parts of the parameter space are within the reach of the HL-LHC for the heavy Higgs boson mass up to M H = 500 GeV. Finally, we expect that extending our current study to the future high energy circular colliders pp(50−100TeV) [117,118] will further enhance the discovery reach of the new heavy Higgs boson H 0 with mass well above 1 TeV through this resonant di-Higgs production channel.  -33 -