Showcasing HH production: Benchmarks for the (HL-)LHC

TL;DR

This paper develops benchmark scenarios for BSM di-Higgs production at the LHC, emphasizing non-SM signatures and multi-fermion channels to guide future experimental searches for new physics.

Contribution

It introduces signature-driven benchmark points for resonant and non-resonant BSM di-Higgs production, focusing on multi-Higgs final states and their unique kinematic features.

Findings

Proposes benchmark points for BSM di-Higgs signatures

Highlights the importance of multi-Higgs final states in BSM searches

Provides guidance for experimental analyses targeting new physics

Abstract

Current projections suggest that the LHC will have only limited sensitivity to di-Higgs production in the Standard Model (SM), possibly even after the completion of its high luminosity phase. Multi-Higgs final states play a fundamental role in many extensions of the SM as they are intrinsically sensitive to modifications of the Higgs sector. Therefore, any new observation in multi-Higgs final states could be linked to a range of beyond the SM (BSM) phenomena that are not sufficiently addressed by the SM. Extensions of the Higgs sector typically lead to new phenomenological signatures in multi-Higgs final states that are vastly different from the SM expectation. In this work, we provide a range of signature-driven benchmark points for resonant and non-resonant BSM di-Higgs production that motivate non-SM kinematic correlations and multi-fermion discovery channels. Relying on theoretically well-motivated assumptions, special attention is devoted to the particular case where the presence of new physics will dominantly manifest itself in multi-Higgs final states.

Tables

Table 1. Table 1: The four Yukawa types of the ℤ 2 subscript ℤ 2 \mathbb{Z}_{2} -symmetric 2HDM defined by the Higgs doublet that couples to each kind of fermions.

	$u$-type	$d$-type	leptons
type I (T1)	${\Phi }_2$	${\Phi }_2$	${\Phi }_2$
type II (T2)	${\Phi }_2$	${\Phi }_1$	${\Phi }_1$
lepton-specific	${\Phi }_2$	${\Phi }_2$	${\Phi }_1$
flipped	${\Phi }_2$	${\Phi }_1$	${\Phi }_2$

Table 2. Table 2: Input parameters for the C2HDM scan, all parameters varied independently between the given minimum and maximum values. The two minimum values of the charged Higgs mass range refer to the scan in the C2HDM T1 and T2, respectively. For more details, see text.

	$t_{\beta }$	${\alpha }_{1,2,3}$	$\text{Re}\left(m_{12}^2\right)$ [TeV2]	$m_{H^{\pm }}$ [TeV]	$m_{H_i\ne h}$ [TeV]
min	0.8	$-\frac{\pi }{2}$	0	0.15/0.59	0.01
max	20	$\frac{\pi }{2}$	0.5	1.5	1.5

Table 3. Table 3: Input parameters for the NMSSM scan, all parameters varied independently between the given minimum and maximum values.

	$t_{\beta }$	$\lambda$	$\kappa$	$M_1$	$M_2$	$M_3$	$A_t$	$A_b$	$A_{\tau }$	$m_{{\tilde{Q}}_3}$	$m_{{\tilde{L}}_3}$	$A_{\lambda }$	$A_{\kappa }$	${\mu }_{\text{eff}}$
				in TeV
min	1	0	-0.7	0.1	0.2	1.3	-6	-6	-3	0.6	0.6	-2	-2	-5
max	50	0.7	0.7	1	2	7	6	6	3	4	4	2	2	5

Table 4. Table 4: Maximum cross section values in fb for LO gluon fusion into Higgs pairs, σ ​ ( g ​ g → H i ​ H j ) 𝜎 → 𝑔 𝑔 subscript 𝐻 𝑖 subscript 𝐻 𝑗 \sigma(gg\to H_{i}H_{j}) , in the C2HDM T1 and T2, with an exclusion luminosity ≥ 64 absent 64 \geq 64  fb -1 that satisfy all theoretical and experimental constraints described above.

$H_i{H}_j$/model	T1	T2
$hh$	794	34.2
$h{H}_{\downarrow }$	49.17	11.38
$h{H}_{\uparrow }$	17.65	10.84
$H_{\downarrow }{H}_{\downarrow }$	3196	0.18
$H_{\downarrow }{H}_{\uparrow }$	12.58	0.11
$H_{\uparrow }{H}_{\uparrow }$	7.10	0.18

Table 5. Table 5: NMSSM: Maximum cross section values in fb for LO gluon fusion into Higgs pairs, σ ​ ( g ​ g → H i ​ H j ) 𝜎 → 𝑔 𝑔 subscript 𝐻 𝑖 subscript 𝐻 𝑗 \sigma(gg\to H_{i}H_{j}) with an exclusion luminosity ≥ 64 absent 64 \geq 64  fb -1 that satisfy all theoretical and experimental constraints described above.

$H_i{H}_j$	NMSSM
$hh$	67
$h{H}_{\downarrow }$	26
$h{A}_{\downarrow }$	493
$h{H}_{\uparrow }$	25
$H_{\downarrow }{H}_{\downarrow }$	4114
$H_{\downarrow }{H}_{\uparrow }$	1.20
$H_{\uparrow }{H}_{\uparrow }$	0.09
$A_{\downarrow }{A}_{\downarrow }$	15894

Table 6. Table 6: T1: Lines 1-8: The input parameters of the benchmark points T1BP1 - 6 . Lines 9 to 13: The derived 3rd neutral Higgs boson mass, the CP-odd admixtures R i ​ 3 2 superscript subscript 𝑅 𝑖 3 2 R_{i3}^{2} and the exclusion luminosity ℒ excl subscript ℒ excl {\cal L}_{\text{excl}} . Lines 14 and 15: The NLO QCD gluon fusion h ​ h ℎ ℎ hh production cross section at s = 14 𝑠 14 \sqrt{s}=14  TeV and the corresponding K 𝐾 K -factor.

	T1BP1	T1BP2	T1BP3	T1BP4	T1BP5	T1BP6
$m_{H_1}$ [GeV]	125.09	125.09	125.09	119.73	125.09	62.67
$m_{H_2}$ [GeV]	130.24	131.52	233.86	125.09	265.60	125.09
$m_{H^{\pm }}$ [GeV]	169.99	282.75	164.87	185.41	307.47	164.35
$\text{Re}\left(m_{12}^2\right)$ [GeV2]	679	12376	11473	7522	11435	130
${\alpha }_1$	1.300	1.249	1.268	1.276	1.246	-0.145
${\alpha }_2$	-0.075	-0.032	0.00262	1.494	$7.125\cdot {10}^{-3}$	-0.0536
${\alpha }_3$	1.306	1.570	-0.809	-1.460	-1.478	-0.0650
$\tan \beta$	4.05	3.23	3.32	5.30	5.54	8.26
$m_{H_3}$ [GeV]	132.95	290.17	234.51	211.43	279.70	138.01
$R_{13}^2$	$5.558\cdot {10}^{-3}$	$1.027\cdot {10}^{-3}$	$6.863\cdot {10}^{-6}$	0.994	$5.077\cdot {10}^{-5}$	$2.870\cdot {10}^{-3}$
$R_{23}^2$	$0.926$	$0.999$	0.523	$5.819\cdot {10}^{-3}$	0.991	$4.212\cdot {10}^{-3}$
$R_{33}^2$	0.068	1.217$\cdot {10}^{-6}$	0.477	$7.140\cdot {10}^{-5}$	$8.611\cdot {10}^{-3}$	0.993
${ℒ}_{\text{excl}}$ [fb-1]	11500	1641	1365	2628	1082	2579
${\sigma }_{hh}^{\text{NLO}}$ [fb]	36.52	36.59	37.88	38.21	897.74	37.26
$K$-factor	1.95	1.95	1.95	1.95	2.06	1.95

Table 7. Table 7: C2HDM T1: T1BP1 - 6 : The rates ( x ​ x ) ​ ( y ​ y ) H i ​ H j 𝑥 𝑥 subscript 𝑦 𝑦 subscript 𝐻 𝑖 subscript 𝐻 𝑗 (xx)(yy)_{H_{i}H_{j}} into the di-Higgs states [ H i ​ H j ] delimited-[] subscript 𝐻 𝑖 subscript 𝐻 𝑗 [H_{i}H_{j}] normalized to the SM from NLO QCD di-Higgs production through gluon fusion at s = 14 𝑠 14 \sqrt{s}=14  TeV in the final state ( x ​ x ) ​ ( y ​ y ) 𝑥 𝑥 𝑦 𝑦 (xx)(yy) : [ p ​ p → H i ​ H j → ( x ​ x ) ​ ( y ​ y ) ] / [ p ​ p → H SM ​ H SM → ( x ​ x ) ​ ( y ​ y ) ] delimited-[] → 𝑝 𝑝 subscript 𝐻 𝑖 subscript 𝐻 𝑗 → 𝑥 𝑥 𝑦 𝑦 delimited-[] → 𝑝 𝑝 superscript 𝐻 SM superscript 𝐻 SM → 𝑥 𝑥 𝑦 𝑦 [pp\to H_{i}H_{j}\to(xx)(yy)]/[pp\to H^{\text{SM}}H^{\text{SM}}\to(xx)(yy)] .

	T1BP1	T1BP2	T1BP3	T1BP4	T1BP5	T1BP6
$\left(b\overline{b}\right){\left(b\overline{b}\right)}_{H_i{H}_j}$ [fb]	$\left[hh\right]$: 1.05	$\left[h{H}_{\downarrow }\right]$: 1.69	$\left[h{H}_{\downarrow }\right]$: 0.026	$\left[h{H}_{\downarrow }\right]$: 0.66	$\left[hh\right]$: 23.80	$\left[H_{\downarrow }{H}_{\downarrow }\right]$: 145
$\left(b\overline{b}\right){\left(\tau \overline{\tau }\right)}_{H_i{H}_j}$ [fb]	$\left[hh\right]$: 1.04	$\left[h{H}_{\downarrow }\right]$: 1.70	$\left[h{H}_{\downarrow }\right]$: 0.027	$\left[h{H}_{\downarrow }\right]$: 0.66	$\left[hh\right]$: 23.51	$\left[H_{\downarrow }{H}_{\downarrow }\right]$: 124
$\left(b\overline{b}\right){\left(\gamma \gamma \right)}_{H_i{H}_j}$ [fb]	$\left[hh\right]$: 0.91	$\left[h{H}_{\downarrow }\right]$: 0.97	$\left[h{H}_{\downarrow }\right]$: 0.019	$\left[h{H}_{\downarrow }\right]$: 0.41	$\left[hh\right]$: 24.32	$\left[H_{\downarrow }{H}_{\downarrow }\right]$: 0.29

Table 8. Table 8: T1 heavy spectrum: Lines 1-8: The input parameters of the benchmark points T1BP1_H , T1BP2_H and T1BP5_H . Lines 9 to 13: The derived 3rd neutral Higgs boson mass, the CP-odd admixtures R i ​ 3 2 superscript subscript 𝑅 𝑖 3 2 R_{i3}^{2} and the exclusion luminosity ℒ excl subscript ℒ excl {\cal L}_{\text{excl}} . Lines 14 and 15: The NLO QCD gluon fusion h ​ h ℎ ℎ hh production cross section at s = 14 𝑠 14 \sqrt{s}=14  TeV and the corresponding K 𝐾 K -factor.

	T1BP1_H	T1BP2_H	T1BP5_H
$m_{H_1}$ [GeV]	125.09	125.09	125.09
$m_{H_2}$ [GeV]	407.30	364.98	397.67
$m_{H^{\pm }}$ [GeV]	410.24	436.72	428.96
$\text{Re}\left(m_{12}^2\right)$ [GeV2]	17152	39827	17992
${\alpha }_1$	1.406	1.291	1.379
${\alpha }_2$	-5.946$\cdot {10}^{-3}$	3.132$\cdot {10}^{-3}$	$5.435\cdot {10}^{-3}$
${\alpha }_3$	0.244	-1.550	-0.685
$\tan \beta$	9.20	3.40	8.28
$m_{H_3}$ [GeV]	425.13	401.33	403.92
$R_{13}^2$	$3.535\cdot {10}^{-5}$	$9.809\cdot {10}^{-6}$	$2.950\cdot {10}^{-5}$
$R_{23}^2$	$0.059$	0.999	0.400
$R_{33}^2$	0.941	4.330$\cdot {10}^{-4}$	0.600
${ℒ}_{\text{excl}}$ [fb-1]	2461	1792	1590
${\sigma }_{hh}^{\text{NLO}}$ [fb]	206.53	43.81	400.60
$K$-factor	1.98	1.98	1.99

Table 9. Table 9: C2HDM T2: Lines 1-8: The input parameters of the benchmark points T2BP1 and T2BP2 . Lines 9 to 13: The derived 3rd neutral Higgs boson mass, the CP-odd admixtures R i ​ 3 2 superscript subscript 𝑅 𝑖 3 2 R_{i3}^{2} and the exclusion luminosity ℒ excl subscript ℒ excl {\cal L}_{\text{excl}} . Lines 14 and 15: The NLO QCD gluon fusion h ​ h ℎ ℎ hh production cross section at s = 14 𝑠 14 \sqrt{s}=14  TeV and the corresponding K 𝐾 K -factor.

	T2BP1	T2BP2
$m_{H_1}$ [GeV]	125.09	125.09
$m_{H_2}$ [GeV]	858.09	814.56
$m_{H^{\pm }}$ [GeV]	835.85	894.84
$\text{Re}\left(m_{12}^2\right)$ [GeV2]	252703	227697
${\alpha }_1$	1.141	1.042
${\alpha }_2$	$-5.268\cdot {10}^{-4}$	$6.184\cdot {10}^{-4}$
${\alpha }_3$	1.198	-1.157
$\tan \beta$	2.16	1.71
$m_{H_3}$ [GeV]	858.65	814.94
$R_{13}^2$	$2.775\cdot {10}^{-7}$	$3.824\cdot {10}^{-7}$
$R_{23}^2$	$0.867$	0.832
$R_{33}^2$	0.133	0.162
${ℒ}_{\text{excl}}$	2664	2016
${\sigma }_{hh}^{\text{NLO}}$ [fb]	37.82	38.02
$K$-factor	1.95	1.95

Table 10. Table 10: NMSSM: Lines 1-14: The input parameters of the benchmark points NMBP1 - 4 . Lines 15 to 20: The derived Higgs boson masses. Line 21: The exclusion luminosity. Lines 22-23: The NLO QCD gluon fusion h ​ h ℎ ℎ hh production cross section at s = 14 𝑠 14 \sqrt{s}=14  TeV and its corresponding K 𝐾 K -factor.

	NMBP1	NMBP2	NMBP3	NMBP4
$M_1$ [GeV]	638	457	608	313
$M_2$ [GeV]	1254	386	546	569
$M_3$ [GeV]	4169	6345	6778	3485
$A_t$ [GeV]	2456	5134	1092	532
$A_b$ [GeV]	-2213	-2908	-4015	2009
$A_{\tau }$ [GeV]	1443	-667	2370	354
$M_{{\tilde{Q}}_3}$ [GeV]	1293	3175	2574	3581
$M_{{\tilde{L}}_3}$ [GeV]	1147	1276	790	1188
$\tan \beta$	1.96	1.87	1.68	1.49
$\lambda$	0.55	0.50	0.60	0.54
$\kappa$	0.43	0.47	0.33	0.27
$A_{\lambda }$ [GeV]	-55	33	425	-416
$A_{\kappa }$ [GeV]	373	358	-672	667
${\mu }_{\text{eff}}$ [GeV]	-293	-299	321	-327
$m_{H_1}$ [GeV]	124.34	124.11	101.13	39.52
$m_{H_2}$ [GeV]	335.27	409.70	125.88	125.64
$m_{H_3}$ [GeV]	530.39	465.57	627.95	634.32
$m_{A_1}$ [GeV]	487.34	406.66	608.57	580.42
$m_{A_2}$ [GeV]	540.58	553.09	624.77	631.97
$m_{H^{\pm }}$ [GeV]	520.47	426.32	621.81	628.97
${ℒ}_{\text{excl}}$ [fb-1]	1370	118	1192	101
${\sigma }_{hh}^{\text{NLO}}$ [fb]	69.29	131.83	43.62	42.31
$K$-factor	1.97	1.97	1.96	1.96

Table 11. Table 11: NMSSM NMBP1 - 4 : The rates ( x ​ x ) ​ ( y ​ y ) H i ​ H j 𝑥 𝑥 subscript 𝑦 𝑦 subscript 𝐻 𝑖 subscript 𝐻 𝑗 (xx)(yy)_{H_{i}H_{j}} into the di-Higgs states [ H i ​ H j ] delimited-[] subscript 𝐻 𝑖 subscript 𝐻 𝑗 [H_{i}H_{j}] normalized to the SM from NLO QCD di-Higgs production through gluon fusion at s = 14 𝑠 14 \sqrt{s}=14  TeV in the final state ( x ​ x ) ​ ( y ​ y ) 𝑥 𝑥 𝑦 𝑦 (xx)(yy) : [ p ​ p → H i ​ H j → ( x ​ x ) ​ ( y ​ y ) ] / [ p ​ p → H SM ​ H SM → ( x ​ x ) ​ ( y ​ y ) ] delimited-[] → 𝑝 𝑝 subscript 𝐻 𝑖 subscript 𝐻 𝑗 → 𝑥 𝑥 𝑦 𝑦 delimited-[] → 𝑝 𝑝 superscript 𝐻 SM superscript 𝐻 SM → 𝑥 𝑥 𝑦 𝑦 [pp\to H_{i}H_{j}\to(xx)(yy)]/[pp\to H^{\text{SM}}H^{\text{SM}}\to(xx)(yy)] .

	NMBP1	NMBP2	NMBP3	NMBP4
$\left(b\overline{b}\right){\left(b\overline{b}\right)}_{H_i{H}_j}$ [fb]	$\left[hh\right]$: 2.35	$\left[hh\right]$: 4.77	$\left[h{H}_{\downarrow }\right]$: 2.33	$\left[H_{\downarrow }{H}_{\downarrow }\right]$: 541.60
$\left(b\overline{b}\right){\left(\tau \overline{\tau }\right)}_{H_i{H}_j}$ [fb]	$\left[hh\right]$: 2.31	$\left[hh\right]$: 4.70	$\left[h{H}_{\downarrow }\right]$: 2.27	$\left[H_{\downarrow }{H}_{\downarrow }\right]$: 432.18
$\left(b\overline{b}\right){\left(\gamma \gamma \right)}_{H_i{H}_j}$ [fb]	$\left[hh\right]$: 2.10	$\left[hh\right]$: 3.79	$\left[h{H}_{\downarrow }\right]$: 1.16	$\left[H_{\downarrow }{H}_{\downarrow }\right]$: 7.11

Table 12. Table 12: NMSSM: Lines 1-14: The input parameters of the benchmark points NMBP5 , 6 , 8 ( NMBP7 = NMBP5 ). Lines 15-20: The derived Higgs boson masses. Line 21: The exclusion luminosity. Lines 22-23: The NLO QCD gluon fusion h ​ h ℎ ℎ hh production cross section at s = 14 𝑠 14 \sqrt{s}=14  TeV and its corresponding K 𝐾 K -factor.

	NMBP5	NMBP6	NMBP8
$M_1$ [GeV]	455	842	321
$M_2$ [GeV]	1741	1510	749
$M_3$ [GeV]	4179	1924	2060
$A_t$ [GeV]	-5923	-4901	5449
$A_b$ [GeV]	-3994	-3817	1232
$A_{\tau }$ [GeV]	-1773	-226	2253
$M_{{\tilde{Q}}_3}$ [GeV]	2391	3539	2344
$M_{{\tilde{L}}_3}$ [GeV]	2117	1623	1163
$\tan \beta$	2.17	1.46	1.98
$\lambda$	0.53	0.55	0.49
$\kappa$	0.44	0.42	0.46
$A_{\lambda }$ [GeV]	-177	-4.78	60.12
$A_{\kappa }$ [GeV]	47	6.32	-0.81
${\mu }_{\text{eff}}$ [GeV]	-327	-307	418
$m_{H_1}$ [GeV]	125.62	125.59	125.07
$m_{H_2}$ [GeV]	504.14	428.17	618.30
$m_{H_3}$ [GeV]	612.42	465.10	767.13
$m_{A_1}$ [GeV]	220.91	106.06	37.38
$m_{A_2}$ [GeV]	602.07	444.01	629.92
$m_{H^{\pm }}$ [GeV]	597.20	432.53	620.36
${ℒ}_{\text{excl}}$ [fb-1]	1153	152	109
${\sigma }_{hh}^{\text{NLO}}$ [fb]	53.42	97.53	54.54
$K$-factor	1.96	1.97	1.96

Table 13. Table 13: NMSSM NMBP5 , 6 , 8 ( NMBP7 = NMBP5 ): The rates ( x ​ x ) ​ ( y ​ y ) H i ​ H j 𝑥 𝑥 subscript 𝑦 𝑦 subscript 𝐻 𝑖 subscript 𝐻 𝑗 (xx)(yy)_{H_{i}H_{j}} into the di-Higgs states [ H i ​ H j ] delimited-[] subscript 𝐻 𝑖 subscript 𝐻 𝑗 [H_{i}H_{j}] normalized to the SM from NLO QCD di-Higgs production through gluon fusion at s = 14 𝑠 14 \sqrt{s}=14  TeV in the final state ( x ​ x ) ​ ( y ​ y ) 𝑥 𝑥 𝑦 𝑦 (xx)(yy) : [ p ​ p → H i ​ H j → ( x ​ x ) ​ ( y ​ y ) ] / [ p ​ p → H SM ​ H SM → ( x ​ x ) ​ ( y ​ y ) ] delimited-[] → 𝑝 𝑝 subscript 𝐻 𝑖 subscript 𝐻 𝑗 → 𝑥 𝑥 𝑦 𝑦 delimited-[] → 𝑝 𝑝 superscript 𝐻 SM superscript 𝐻 SM → 𝑥 𝑥 𝑦 𝑦 [pp\to H_{i}H_{j}\to(xx)(yy)]/[pp\to H^{\text{SM}}H^{\text{SM}}\to(xx)(yy)] .

	NMBP5		NMBP6	NMBP8
$\left(b\overline{b}\right){\left(b\overline{b}\right)}_{H_i{H}_j}$ [fb]	$\left[h{A}_{\downarrow }\right]$: 2.56	$\left[A_{\downarrow }{A}_{\downarrow }\right]$: 4.57	$\left[h{A}_{\downarrow }\right]$: 44.37	$\left[A_{\downarrow }{A}_{\downarrow }\right]$: 2127
$\left(b\overline{b}\right){\left(\tau \overline{\tau }\right)}_{H_i{H}_j}$ [fb]	$\left[h{A}_{\downarrow }\right]$: 2.68	$\left[A_{\downarrow }{A}_{\downarrow }\right]$: 5.07	$\left[h{A}_{\downarrow }\right]$: 43.22	$\left[A_{\downarrow }{A}_{\downarrow }\right]$: 1710
$\left(b\overline{b}\right){\left(\gamma \gamma \right)}_{H_i{H}_j}$ [fb]	$\left[h{A}_{\downarrow }\right]$: 2.06	$\left[A_{\downarrow }{A}_{\downarrow }\right]$: 3.20	$\left[h{A}_{\downarrow }\right]$: 22.65	$\left[A_{\downarrow }{A}_{\downarrow }\right]$: 11.30

Equations

$$\begin{array}{rl}\mathcal{L}=& \frac{{\alpha }_s}{12\pi }{G}^{a\mu \nu }{G}_{\mu \nu }^a\log \left(1+\frac{h}{v}\right)\\ & =\frac{{\alpha }_s}{12\pi }{G}^{a\mu \nu }{G}_{\mu \nu }^a\left(\frac{h}{v}-\frac{h^2}{v^2}\right)+\mathcal{O}(h^3)\end{array}$$

$$V$$

$${\Phi }_1=\left(\begin{array}{c}{\varphi }_1^{+}\\ \frac{v_1+{\rho }_1+i{\eta }_1}{\sqrt{2}}\end{array}\right)\text{mbox}and{\Phi }_2=\left(\begin{array}{c}{\varphi }_2^{+}\\ \frac{v_2+{\rho }_2+i{\eta }_2}{\sqrt{2}}\end{array}\right),$$

$$\tan \beta \equiv t_{\beta }=\frac{v_2}{v_1}.$$

$$\displaystyle\left(\begin{array}[]{c}H_{1}\\ H_{2}\\ H_{3}\end{array}\right)=R\left(\begin{array}[]{c}\rho_{1}\\ \rho_{2}\\ \rho_{3}\end{array}\right)\,.$$

$$R{\mathcal{M}}^2{R}^T=\text{diag}\left(m_{H_1}^2,m_{H_2}^2,m_{H_3}^2\right),$$

$$-\frac{\pi }{2}\le {\alpha }_i<\frac{\pi }{2},$$

$$R=\left(\begin{array}[]{ccc}c_{1}c_{2}&s_{1}c_{2}&s_{2}\\ -(c_{1}s_{2}s_{3}+s_{1}c_{3})&c_{1}c_{3}-s_{1}s_{2}s_{3}&c_{2}s_{3}\\ -c_{1}s_{2}c_{3}+s_{1}s_{3}&-(c_{1}s_{3}+s_{1}s_{2}c_{3})&c_{2}c_{3}\end{array}\right)\;.$$

$$v,t_{\beta },{\alpha }_{1,2,3},m_{H_i},m_{H_j},m_{H^{\pm }},\text{mbox}Re(m_{12}^2).$$

$$\mathcal{W}$$

$$-{\mathcal{L}}_{\mathrm{mass}}$$

$$-{\mathcal{L}}_{\mathrm{tril}}$$

$$-{\mathcal{L}}_{\mathrm{gauginos}}$$

$$H_d$$

$$H_u$$

$$(H_1,H_2,H_3{)}^T={\mathcal{R}}^S(h_d,h_u,h_s{)}^T.$$

$$(A_1,A_2,G{)}^T={\mathcal{R}}^P{\mathcal{R}}^G(a_d,a_u,a_s{)}^T.$$

$$\lambda ,\kappa ,A_{\lambda },A_{\kappa },\tan \beta =v_u/v_d,{\mu }_{\mathrm{eff}}=\lambda v_s/\sqrt{2}.$$

$$10\text{mbox}GeV\le m_{H_j}<1.5\text{mbox}TeV.$$

$$\alpha (M_Z)$$

$$\displaystyle\begin{array}[]{lcllcl}m_{e}&=&0.510998928\mbox{ MeV}\;,&m_{\mu}&=&105.6583715\mbox{ MeV}\;,\\ m_{u}&=&100\mbox{ MeV}\;,&m_{d}&=&100\mbox{ MeV}\;,\\ m_{s}&=&100\mbox{ MeV}\;.\end{array}$$

$$m_{{\tilde{t}}_R}=m_{{\tilde{Q}}_3},m_{{\tilde{\tau }}_R}=m_{{\tilde{L}}_3}\text{mbox}and{m}_{{\tilde{b}}_R}=3\text{mbox}TeV.$$

$$m_{{\tilde{u}}_R,{\tilde{c}}_R}=m_{{\tilde{d}}_R,{\tilde{s}}_R}=m_{{\tilde{Q}}_{1,2}}=m_{{\tilde{L}}_{1,2}}=m_{{\tilde{e}}_R,{\tilde{\mu }}_R}=3\text{mbox}TeV.$$

$${\lambda }^2+{\kappa }^2<0.7^2.$$

$${\mu }_R=M_s=\sqrt{m_{{\tilde{Q}}_3}{m}_{{\tilde{t}}_R}}.$$

$$124\text{mbox}GeV\le m_h\le 126\text{mbox}GeV,$$

$$\displaystyle\begin{array}[]{lll}\mbox{T1:}&\quad\sigma(hh)^{\text{NLO}}=1.64\mbox{ pb},&\quad K=2.06,\\ \mbox{T2:}&\quad\sigma(hh)^{\text{NLO}}=64.96\mbox{ fb},&\quad K=1.94.\end{array}$$

$${\Sigma }_X=\sum _{i\in \text{SM}\backslash \{h\}}\text{BR}(X\to i),$$

$$\sigma (gg\to h{H}_{\downarrow }\to 4b)/\text{mbox}SM$$

$$\sigma (gg\to h{H}_{\downarrow }\to (2b)(2\tau ))/\text{mbox}SM$$

$$\sigma (gg\to h{H}_{\downarrow }\to (2b)(2\gamma ))/\text{mbox}SM$$

$$\sigma (gg\to hh\to 4b)/\text{mbox}SM=23.80,$$

Full text

I Introduction

The Higgs precision spectroscopy program that ensued after the discovery of the Higgs boson in 2012 Aad et al. (2012); Chatrchyan et al. (2012) has assumed a central role in particle physics over the past years. One reason why the Higgs couplings and properties measurements move increasingly into the focus of searches for physics beyond the Standard Model (BSM) is the lack of conclusive hints for new interactions in the plethora of BSM searches performed by the ATLAS and CMS experiments. The Higgs boson as the direct implication of electroweak symmetry breaking is typically considered as a harbinger of new physics due to its special role in the unitarization of scattering amplitudes at high energy Cornwall et al. (1973, 1974); Lee et al. (1977a, b) and its relation to naturalness of the electroweak scale Veltman (1981), to only name a couple of examples.

Although modifications of Higgs physics at the TeV scale of this size are still well within the limits set by recent 13 TeV LHC measurements, now that the SM can be considered as complete no additional ultra-violet (UV) energy scale can be predicted from the SM alone. This becomes even more pressing when coupling measurements stay consistent with the SM expectation in the future.

Additional requirements are more conveniently imposed in model-specific approaches, which try to mend apparent shortcomings of the SM such as the lack of a viable dark matter candidate or an insufficiently first-order electroweak phase transition to address criteria of baryogenesis Sakharov (1967). While model-independent approaches based on effective field theory (EFT) Grzadkowski et al. (2010) can inform UV completions that address these questions through matching calculations, the appearance of novel phenomenological signatures such as resonances or thresholds within the LHC’s kinematic coverage typically falls outside the region of reliability of these techniques.

A process that highlights the shortcomings of EFT methods in the presence of thresholds already in the SM context is di-Higgs production $pp\to hh$ (see e.g. de Florian et al. (2016); Dawson et al. (2018) for recent summaries). Therefore, EFT methods for $ggh^{n}$ that can be summarized as Ellis et al. (1976); Shifman et al. (1979); Vainshtein et al. (1980); Voloshin (1986); Kniehl and Spira (1995)

[TABLE]

reflect the destructive interference between the top triangle and box diagrams of Fig. 1. While it can be assumed that the Higgs-top coupling can be accessed at the LHC in the range of $\sim 10\%$ and the intermediate top mass scale is under theoretical control at the next-to-leading order (NLO) level***A recent summary of higher-order corrections to Higgs pair production can be found in Baglio et al. (2018). Borowka et al. (2016a, b); Heinrich et al. (2017); Baglio et al. (2018); Davies et al. (2018); Bonciani et al. (2018); Chen et al. (2015), the situation for the trilinear coupling is less clear as it probes a direction in the dimension-6 linear EFT space $\sim c_{6}(H^{\dagger}H)^{3}/\Lambda^{2}$ (we denote the SM Higgs doublet with $H$ in the following) Barger et al. (2003). As $c_{6}$ is essentially a free parameter unless a matching calculation is performed, its size is only limited by technical considerations related to perturbation theory, on which we need to rely to make phenomenological predictions. This has raised the question of how large trilinear coupling modifications can be Baglio et al. (2013); Gupta et al. (2013); Di Vita et al. (2017); Di Luzio et al. (2017) to inform di-Higgs investigations. Current CMS projections show that a sensitivity of $\lesssim 0.5~{}\sigma$ at $3~{}{\text{ab}}^{-1}$ Khachatryan et al. (2017) seems a realistic target at the LHC. Recent theoretical studies are slightly more optimistic Adhikary et al. (2018); Goncalves et al. (2018); Homiller and Meade (2018), but the Higgs trilinear coupling will only understood in great detail even at large luminosity, in particular in the light of possible top quark Yukawa coupling modifications Azatov et al. (2015). As a consequence, di-Higgs production is a main motivation for considering a major energy upgrade of the LHC or a future hadron collider Yao (2013); Barr et al. (2015); Papaefstathiou (2015); Zhao et al. (2017); Contino et al. (2017); Banerjee et al. (2018).

In concrete UV scenarios that address fundamental BSM questions, typically a number of exotic states appear in the spectrum that will not only impact the “standard” single Higgs phenomenology at the LHC but can lead to even more dramatic changes in the production of multiple Higgs bosons. For instance, the triangle and box diagrams probe different aspects of top physics in extended top sectors Grober and Muhlleitner (2011); Dolan et al. (2013); Chen et al. (2014); Grober et al. (2016); Vryonidou and Zhang (2018).

Interleaving modifications of single Higgs physics with theoretically well-motivated UV considerations can therefore turn di-Higgs production into a strong probe of new physics: new kinematic features can appear that motivate new final states and search strategies, which might currently be unconsidered; di-Higgs production can be enhanced or suppressed. Such phenomenological modifications become particularly relevant when extrapolations of standard single Higgs channels do not show a significant departure from their SM expectation in these scenarios.

We will address these questions in this work using a particular set of models that allow us to contrast precise theoretical and phenomenological requirements with concrete predictions of single- and multi-Higgs production. Imposing, e.g., a strong first-order electroweak phase transition, dark matter constraints, electric dipole measurements and consistency with current Higgs coupling measurements as well as an extrapolation thereof, we discuss the results of a comprehensive scan of the models’ parameter space with a particular emphasis on the relevance of multi-Higgs final states. We distill this scan into a number of representative benchmark points of BSM theories that highlight the importance of di-Higgs measurements in the future. In passing, we discuss how non-SM signatures are correlated with modifications of single Higgs physics BSM effects (or a lack of the latter). In turn, this also allows us to formulate an upper limit of the SM-like di-Higgs production cross section in these models when there are no conclusive hints for new physics in single Higgs phenomenology.

This work is structured as follows: In Sec. II, we outline the models that we consider for the purpose of this work. These are the CP-violating 2-Higgs-Doublet Model (C2HDM) and the Next-to-Minimal Supersymmetric Extension (NMSSM). Both models are special in the sense that they feature extended Higgs sectors that allow for Higgs-to-Higgs decays, also into final states with different Higgs bosons, or even cascade Higgs-to-Higgs decays. In Sec. III we outline the details of our scan over the model space and discuss the overall scan results. We present our benchmarks for the C2HDM and NMSSM together with their phenomenological properties in Sec. IV. We conclude in Sec. V.

II The models

II.1 The C2HDM

The 2-Higgs-Doublet Model (2HDM) Lee (1973); Gunion et al. (2000); Branco et al. (2012) is obtained from the SM by adding a second $SU(2)_{L}$ Higgs doublet. The Higgs potential of a general 2HDM with a softly broken $\mathbb{Z}_{2}$ symmetry, under which $\Phi_{1}\to\Phi_{1}$ and $\Phi_{2}\to-\Phi_{2}$, can be written as

[TABLE]

The absence of Flavour Changing Neutral Currents (FCNC) is guaranteed by extending the $\mathbb{Z}_{2}$ symmetry to the fermions. Depending on the $\mathbb{Z}_{2}$ charge assignments, there are four phenomenologically different types of 2HDMs that are summarized in Tab. 1.

Hermiticity of the potential requires all parameters to be real, except for $\lambda_{5}$ and $m_{12}^{2}$. If they have different unrelated complex phases we are in the framework of the complex or CP-violating 2HDM Ginzburg et al. (2002), which depends on ten real parameters. In the description of the C2HDM we follow the conventions of Fontes et al. (2014). The vacuum expectation values (VEVs) developed by the Higgs doublets after electroweak symmetry breaking (EWSB) could in principle be complex in the C2HDM. Since the phase can be removed by a basis change Ginzburg et al. (2002), we set it to zero without loss of generality. In terms of the complex charged fields $\phi_{i}^{+}$ and the real neutral CP-even and CP-odd fields $\rho_{i}$ and $\eta_{i}$ ($i=1,2$), respectively, the Higgs doublets are given by

[TABLE]

where $v_{1}$ and $v_{2}$ denote the VEVs of the two Higgs doublets $\Phi_{1}$ and $\Phi_{2}$, respectively, and $v_{1}^{2}+v_{2}^{2}=v^{2}$ with the SM VEV $v\approx 246$ GeV. The ratio of the VEVs is parametrised by the mixing angle $\beta$,

[TABLE]

The minimum conditions obtained from the requirement that the minimum of the potential is given by $\langle\Phi_{i}\rangle=(0,v_{i}/\sqrt{2})^{T}$, can be used to trade the parameters $m_{11}^{2}$ and $m_{22}^{2}$ for $v_{1}$ and $v_{2}$. They also yield a relation between the imaginary parts of $m_{12}^{2}$ and $\lambda_{5}$, so that one of the ten parameters is fixed. The neutral Higgs mass eigenstates $H_{i}$ ($i=1,2,3$) are obtained from the neutral components of the C2HDM basis, $\rho_{1,2}$ and $\rho_{3}\equiv(1/\sqrt{2})(-\sin\beta\,\eta_{1}+\cos\beta\,\eta_{2})$ via the rotation†††Actually, the field $\rho_{3}$ is equal to the CP-odd component of the second Higgs doublet in the Higgs basis Lavoura and Silva (1994); Botella and Silva (1995).

[TABLE]

The matrix $R$ diagonalizes the mass matrix ${\cal M}$ of the neutral states,

[TABLE]

where $m_{H_{1}}\leq m_{H_{2}}\leq m_{H_{3}}$ denote the masses of the neutral Higgs bosons. Introducing the abbreviations $s_{i}\equiv\sin\alpha_{i}$ and $c_{i}\equiv\cos\alpha_{i}$ with

[TABLE]

the mixing matrix $R$ can be parametrised as

[TABLE]

In total, the C2HDM has nine independent parameters that we choose to be El Kaffas et al. (2007)

[TABLE]

The $m_{H_{i}}$ and $m_{H_{j}}$ denote any of the three neutral Higgs boson masses. The third mass is not independent and is calculated from the other parameters El Kaffas et al. (2007). For further details, in particular also on the couplings of the C2HDM, see Fontes et al. (2018).

II.2 The NMSSM

Supersymmetric (SUSY) extensions require the introduction of at least a second Higgs doublet. In the NMSSM, the minimal field content with the doublet superfields $\hat{H}_{u}$ and $\hat{H}_{d}$ is extended by a complex superfield $\hat{S}$ (for reviews on the NMSSM, see for example Ellwanger et al. (2010); Maniatis (2010)). The NMSSM Higgs potential is derived from the superpotential, the soft SUSY breaking Lagrangian and the $D$-term contributions. The scale-invariant NMSSM superpotential reads in terms of the hatted superfields

[TABLE]

For simplicity, we only included here the third generation fermion superfields, given by the left-handed doublet quark ($\hat{Q}_{3}$), and lepton ($\hat{L}_{3}$) superfields as well as right-handed singlet quark ($\hat{t}^{c}_{R},\hat{b}^{c}_{R}$) and lepton ($\hat{\tau}^{c}_{R}$) superfields. The soft SUSY breaking Lagrangian

[TABLE]

contains the mass terms $m_{x}$ for the Higgs ($x=H_{u},H_{d},S$) and sfermion ($x=\tilde{Q}_{3},\tilde{t}_{R},\tilde{b}_{R},\tilde{L}_{3},\tilde{\tau}_{R}$) fields, obtained from the complex scalar components of the superfields. The Lagrangian with the trilinear soft SUSY breaking interactions $A_{\lambda,\kappa,t,b,\tau}$ between the sfermion and Higgs fields reads

[TABLE]

The contribution to soft SUSY breaking from the gaugino mass parameters $M_{1,2,3}$ of the bino ($\tilde{B}$), winos ($\tilde{W}$) and gluinos ($\tilde{G}$), respectively, is given by

[TABLE]

Expanding the tree-level scalar potential around the non-vanishing VEVs of the Higgs doublet and singlet fields,

[TABLE]

leads to the Higgs mass matrices for the three scalar ($h_{d},h_{u},h_{s}$), the three pseudoscalar ($a_{d},a_{u},a_{s}$) and the charged Higgs states ($h_{u}^{\pm},h_{d}^{\mp}$). We choose the VEVs $v_{u},v_{d}$ and $v_{s}$ to be real and positive. The three CP-even mass eigenstates $H_{i}$ ($i=1,2,3$) are obtained from the interaction states through rotation with the orthogonal matrix ${\cal R}^{S}$ that diagonalizes the $3\times 3$ mass matrix squared of the CP-even fields,

[TABLE]

The mass eigenstates are ordered by ascending mass, $M_{H_{1}}\leq M_{H_{2}}\leq M_{H_{3}}$. The CP-odd mass eigenstates $A_{1}$ and $A_{2}$ are obtained from a rotation ${\cal R}^{G}$ that separates the Goldstone boson, followed by a rotation ${\cal R}^{P}$ into the mass eigenstates,

[TABLE]

They are also ordered by ascending mass, $M_{A_{1}}\leq M_{A_{2}}$. Altogether, the NMSSM Higgs spectrum consists of seven physical Higgs states, three neutral CP-even, two neutral CP-odd and two charged Higgs bosons. We use the three minimisation conditions of the scalar potential to replace the soft SUSY breaking masses squared for $H_{u},H_{d}$ and $S$ in ${\cal L}_{\text{mass}}$ by the remaining parameters of the tree-level potential so that the set of six parameters parametrising the tree-level NMSSM Higgs sector is given by

[TABLE]

The sign conventions are such that $\lambda$ and $\tan\beta$ are positive, while $\kappa,A_{\lambda},A_{\kappa}$ and $\mu_{\text{eff}}$ can take both signs. Note that the Higgs boson masses are not input parameters, but dependent parameters calculated from the input values. The inclusion of higher-order corrections in the Higgs boson masses is crucial here to shift the mass of the SM-like Higgs boson to the observed value of 125 GeV.

III Details of the Scan

III.1 The C2HDM Scan

The benchmark points ‡‡‡2HDM benchmarks for double Higgs production can be found in Refs. Haber and Stål (2015); Baglio et al. (2014a). provided in this paper have to satisfy theoretical and experimental constraints. In order to find valid points, we perform a scan in the C2HDM parameter space and additionally require the mass of one of the Higgs bosons, to be identified with the SM-like one and denoted by $h$, to be $m_{h}=125.09$ GeV Aad et al. (2015a). The scan ranges are summarized in Tab. 2. For simplicity, we only consider the C2HDM Type 1 (T1) and Type 2 (T2), which cover to a large extent the phenomenological effects to be expected in the C2HDM.

Since physical parameter points with $\mbox{Re}(m_{12}^{2})<0$ are extremely rare, though possible, we neglect them in our scan. We test for compatibility with the flavour constraints on $R_{b}$ Haber and Logan (2000); Deschamps et al. (2010) and $B\to X_{s}\gamma$ Deschamps et al. (2010); Mahmoudi and Stal (2010); Hermann et al. (2012); Misiak et al. (2015); Misiak and Steinhauser (2017) as 2$\sigma$ exclusion bounds in the $m_{H^{\pm}}-\tan\beta$ plane. In accordance with Misiak and Steinhauser (2017) we therefore require $m_{H^{\pm}}$ to be above 580 GeV in the C2HDM T2, whereas in the C2HDM T1 this bound is much weaker and depends more strongly on $\tan\beta$. We verify agreement with the electroweak precision data by using the oblique parameters $S$, $T$ and $U$ - the 2HDM formulae are given in Branco et al. (2012); Dawson and Sullivan (2018) - for which we demand $2\sigma$ compatibility with the SM fit Baak et al. (2014), including the full correlation among the three parameters. We require perturbative unitarity to hold at tree level. The third neutral Higgs boson mass $m_{H_{j}\neq H_{i},h}$, which is not an independent input parameter in the C2HDM, but calculated from the other input values, is demanded to lie in the interval

[TABLE]

To avoid degenerate Higgs signals, we additionally impose $m_{H_{i}\neq h}$ to be 5 GeV away from 125 GeV. For the SM input parameters we use Tanabashi et al. (2018); Denner (2015)

[TABLE]

The remaining light quark and lepton masses have been set to Tanabashi et al. (2018); Denner (2015)

[TABLE]

In order to perform the scans and find valid parameter points we use the program ScannerS Coimbra et al. (2013); Ferreira et al. (2014). Besides the above mentioned constraints it also tests for the potential to be bounded from below and uses the tree-level discriminant of Ivanov and Silva (2015) to enforce the electroweak vacuum to be the global minimum of the tree-level Higgs potential. Agreement with the Higgs exclusion limits from LEP, Tevatron and LHC is checked by using HiggsBounds5.2.0 Bechtle et al. (2010, 2011); Bechtle et al. (2014a) and with the Higgs rates by using HiggsSignals2.2.1 Bechtle et al. (2014b). The required decay widths and branching ratios are obtained from the C2HDM implementation C2HDM_HDECAY Fontes et al. (2018) in HDECAY Djouadi et al. (1998, 2018). In the production cross sections we included the QCD corrections taken over from the SM and the Minimal Supersymmetric Extension (MSSM), where available. Electroweak corrections are consistently neglected both in production and decays. For more details on the production cross sections, see Fontes et al. (2018); Muhlleitner et al. (2017).

Working in the C2HDM, we also make sure to be in agreement with the measurements of the electric dipole moment (EDM), where the strongest constraint originates from the electron EDM Inoue et al. (2014). We take the experimental limit given by the ACME collaboration Andreev et al. (2018). Finally, we also investigate for the C2HDM if the parameters of the final data set induce a strong first order phase transition, a necessary condition for successful baryogenesis Sakharov (1967); Quiros (1994); Moore (1999), by using the C++ code BSMPT Basler and Muhlleitner (2018).

III.2 The NMSSM Scan

In order to find benchmark points§§§For NMSSM benchmarks for double Higgs production from Higgs-to-Higgs decays, see King et al. (2014). that are compatible with the recent experimental constraints we proceed as described in Costa et al. (2016); King et al. (2014); Azevedo et al. (2018), where also further details can be found. We perform a scan in the NMSSM parameter range with scan ranges summarized in Tab. 3. The remaining mass parameters of the third generation sfermions not listed in the table are chosen as

[TABLE]

The mass parameters of the first and second generation sfermions are set to

[TABLE]

The soft SUSY breaking trilinear couplings of the first two generations are set equal to the corresponding values of the third generation. In order to ensure perturbativity we applied the rough constraint

[TABLE]

In accordance with the SUSY Les Houches Accord (SLHA) format Skands et al. (2004); Allanach et al. (2009) the soft SUSY breaking masses and trilinear couplings are understood as $\overline{\text{DR}}$ parameters at the scale

[TABLE]

The SM input parameters have been chosen as in the C2HDM scan, with the exception of the top quark mass which has been set to $m_{t}=173.5$ GeV. The small difference of 1 GeV has no effect on the scan results.

The spectrum of the Higgs and SUSY particles including higher-order corrections is calculated with NMSSMTools5.2.0 Ellwanger et al. (2005); Ellwanger and Hugonie (2006, 2007); Das et al. (2012); Mühlleitner et al. (2005); Belanger et al. (2005) which also checks for the constraints from flavour and low-energy observables. It provides the input for HiggsBounds5.2.0 Bechtle et al. (2010, 2011); Bechtle et al. (2014a) to check for compatibility with the exclusion bounds from the Higgs searches. The mass of one of the neutral CP-even Higgs bosons, identified with the SM-like Higgs boson denoted by $h$, has to lie in the range

[TABLE]

and the masses of all other Higgs bosons are demanded to be separated by at least 1 GeV in order to avoid two overlapping signals. The signal strengths of this Higgs boson have to be in agreement with the signal strength fit of Aad et al. (2016) at the $2\times 1\sigma$ level. For the computation of the signal strengths we need the production cross section and branching ratios for the NMSSM Higgs bosons. To compute production through gluon fusion and $b\bar{b}$ annihilation, we take the SM cross sections and multiply them with the effective couplings obtained from NMSSMTools. The SM values are calculated with SusHi Harlander et al. (2013, 2017) and include in gluon fusion the NLO corrections with the full top quark mass dependence Spira et al. (1995) as well as the next-to-next-to-leading order (NNLO) corrections in the heavy quark effective theory Harlander and Kilgore (2002a); Anastasiou and Melnikov (2002); Harlander and Kilgore (2002b); Anastasiou and Melnikov (2003); Ravindran et al. (2003). The next-to-next-to-next-to-leading order (N3LO) corrections are taken into account in a threshold expansion Anastasiou et al. (2015a, b); Anastasiou et al. (2016); Mistlberger (2018) for Higgs masses below 300 GeV. For masses above 50 GeV, $b\bar{b}$ annihilation cross sections that match between the five- and four-flavor scheme are used, obtained in the soft-collinear effective theory Bonvini et al. (2015, 2016). They are in accordance with the results from Forte et al. (2015, 2016). For masses below 50 GeV, cross sections obtained in the Santander matching Harlander et al. (2011) are used, with the five-flavor scheme cross sections from Harlander and Kilgore (2003) and the four-flavor scheme ones from Dittmaier et al. (2004); Dawson et al. (2004); Wiesemann et al. (2015). The branching ratios are taken from NMSSMTools and cross-checked against against NMSSMCALC Baglio et al. (2014b).

The parameter points also have to satisfy the bounds from SUSY searches at the LHC. The gluino mass and the lightest squark mass of the second generation are demanded to lie above 1.85 TeV, respectively, see Aad et al. (2015b). The stop and sbottom masses are required to be above 800 GeV, respectively, Aaboud et al. (2016); Aad et al. (2015b), the slepton masses above 400 GeV Aad et al. (2015b) and the absolute value of the lightest chargino mass above 300 GeV Aad et al. (2015c).

Through an interface with micrOMEGAS Belanger et al. (2005) we obtain the relic density that we require not to exceed the value measured by the PLANCK collaboration Ade et al. (2014). The spin-independent nucleon-dark matter direct detection cross section, also provided by micrOMEGAS, is required not to violate the upper bound from the LUX experiment Akerib et al. (2017). We furthermore test for compatibility with the direct detection limits from XENON1T Aprile et al. (2018) and check the Dark Matter annihilation cross section against the results provided by FermiLat Ackermann et al. (2015).

III.3 Extrapolations

We include a range of extrapolations of single Higgs measurements to our discussion to identify an approximate “exclusion luminosity” (see below) at which single Higgs measurements will start to become sensitive to a particular scenario and spectrum. This notion will allow us to put multi-Higgs final states in direct comparison with single Higgs measurement expectations and identify interesting regions of parameter space.

In particular, we include projections for the $m_{h}\simeq 125~{}\text{GeV}$ standard single Higgs production modes, $gg\to h$ (gluon fusion), $qq\to hjj$ (weak boson fusion), $qq\to Vh$, $V=W^{\pm},Z$ (Higgs radiation), $pp\to t\bar{t}h$ (associated production with top quarks), and consider decays $h\to ZZ$, $h\to WW$, $h\to\gamma\gamma$, $h\to b\bar{b}$ and $h\to\tau^{+}\tau^{-}$. For the decays $h\to\gamma\gamma,ZZ$ we use the CMS projections provided in Ref. Khachatryan et al. (2017); these include the production modes gluon fusion, weak boson fusion and $t\bar{t}h$. We interpolate between different luminosities using a $\sqrt{\cal{L}}$ luminosity dependence at all times.

Projections for $h\to WW$ are obtained using Ref. Aad et al. (2015d) and we rescale these results taking into account the cross section differences between 13 TeV and 8 TeV using the results provided by the Higgs cross section working group Dittmaier et al. (2011). For $h\to b\bar{b}$ we consider extrapolations based on $Vh$ production Aaboud et al. (2017a), $t\bar{t}h$ production Sirunyan et al. (2018a) as well as weak boson fusion Khachatryan et al. (2015a). $h\to\tau\tau$ is based on Sirunyan et al. (2018b), which agrees with the ECFA results of Khachatryan et al. (2017) upon projecting to 3 ab*-1*.

The improved determination of the SM-like Higgs boson with $m_{h}\simeq 125~{}\text{GeV}$ needs to be contrasted with additional coverage of Higgs-like searches for masses $m_{h}\neq 125~{}\text{GeV}$. We include projections of existing resonance searches in $\gamma\gamma$ Aaboud et al. (2017b), $\tau\tau$ Aaboud et al. (2018a), $WW$ Aaboud et al. (2018b) and $ZZ$ Aaboud et al. (2018c) final states. By far the most constraining exotic searches result from $t\bar{t}$ resonance searches, and we extrapolate the results of Ref. Aaboud et al. (2018d). This analysis is performed in the context of a $Z^{\prime}$ model and can therefore be interpreted as only a rough estimate of the search potential of $t\bar{t}$ Higgs resonances. To our knowledge no comprehensive analysis of exotic heavy Higgs masses is publicly available. This is mostly due to the dedicated interference between the background and the Higgs signal that also depends on the CP character of the produced state Aaboud et al. (2017c). This has a significant impact on the sensitivity of $t\bar{t}$ final states. Including such effects is beyond the scope of this work.

III.4 Results

Both in the C2HDM and the NMSSM the enlarged Higgs sector leads to a plethora of di-Higgs production processes. In particular, they feature processes with two different Higgs bosons in the final state. Compared to the SM, the cross sections can be largely enhanced in case of resonant production of a heavy Higgs boson that subsequently decays into a pair of lighter Higgs bosons, provided the Higgs self-coupling is not too small. Moreover, the different Higgs self-couplings themselves can enhance the cross section in view of the well-known fact that in the SM the triangle and box diagrams interfere destructively. Additionally, loops with bottom quarks may play a role in scenarios with enhanced down-type Yukawa couplings for large $\tan\beta$ in the NMSSM or the C2HDM T2. In the NMSSM loops with stop and sbottom quarks also contribute to Higgs pair production, and we furthermore have the possibility to produce a di-Higgs final state with pseudoscalars. These processes can yield even larger rates as has been discussed in detail in Costa et al. (2016). However, due to supersymmetry, the Higgs self-couplings are given in terms of the gauge couplings limiting to some extent deviations from the SM. This is not the case for the C2HDM, so that effects different from the NMSSM may be expected here. In the NMSSM, on the other hand, we have the possibility of the final state Higgs bosons to decay into non-SM final states like e.g. neutralinos, inducing signatures with phenomenologically interesting features. Altogether, both models provide a large playground for possible BSM effects in Higgs pair production that can be rather different.

Furthermore, we restrict ourselves to SM final states. Most of our results show the leading order (LO) Higgs pair production cross sections. For the benchmark points we also computed the NLO QCD corrections in the limit of heavy loop particles. They typically increase the cross section by about a factor two. We have implemented the NLO QCD corrections both for the NMSSM¶¶¶They only include the corrections to the top quark loops. For NLO QCD corrections including also the squarks in the limit of vanishing external momenta, see Agostini et al. (2016). Costa et al. (2016) and C2HDM Grober et al. (2017) in the Fortran code HPAIR∥∥∥See M. Spira’s website, http://tiger.web.psi.ch/proglist.html. that was originally designed to compute the SM and MSSM Higgs pair production at NLO QCD. All Higgs pair production processes have been computed at a c.m. energy of $\sqrt{s}=14$ TeV, and we have adopted the CT14 parton densities Dulat et al. (2016) for the LO and NLO cross sections with $\alpha_{s}(M_{Z})=0.118$ at LO and NLO. The renormalisation scale has been set equal to $M_{HH}/2$, where $M_{HH}$ generically denotes the invariant mass of the final state Higgs pair. Consistent with the application of the heavy top quark limit in the NLO QCD corrections, we neglect the bottom quark loops in the LO cross section.

In view of the possibility of (largely) enhanced production of a pair of SM-like Higgs bosons, in the selection of valid scenarios, we also took into account limits set by LHC $4b$ Aaboud et al. (2018e); Collaboration (2017a); Sirunyan et al. (2018c), $(2b)(2\tau)$ Collaboration (2017b); Sirunyan et al. (2018d); Aaboud et al. (2018f) and $(2b)(2\gamma)$ Collaboration (2017c) final states from the production of a heavy scalar resonance that decays into two 125 GeV Higgs bosons.

III.4.1 C2HDM

We start by discussing the possible sizes of Higgs pair production that are compatible with all present experimental constraints. Table 4 summarizes the maximum cross section values found in the sample of valid parameter points where we additionally applied the extrapolations of Subsection III.3. Taking these into account, we only kept the points that are not excluded at 64 fb*-1*, which corresponds approximately to the present luminosity acquired by the LHC experiments. We will come back to the role of the extrapolations below. In the following, we denote the SM-like Higgs boson with a mass of 125 GeV $h$, the lighter of the non-SM like neutral Higgs bosons is called $H_{\downarrow}$, the heavier one $H_{\uparrow}$. All cross sections are calculated at LO QCD and hence still increase by approximately a factor two when QCD corrections are included.

The SM Higgs pair production cross section in gluon fusion amounts to 19.72 fb at LO and 38.19 fb at NLO******This value differs from the one given in Borowka et al. (2016a). This is because we do not include top quark mass effects here and use a different pdf set. with a $K$-factor of $K=1.95$. Inspection of the Tab. 4 shows, that both in the T1 and T2 scenarios, the maximum attained cross section for $hh$ production can exceed the SM value, in T1 by a factor of about 40 and in T2 by about a factor 1.7. This is also the case at NLO, with NLO cross sections and $K$-factors for $hh$ production in the T1 and T2 models given by

[TABLE]

The reason for the large enhancement in T1 is the resonant production of the heavier Higgs bosons $H_{\downarrow}$ and $H_{\uparrow}$ with a mass of 285 GeV and 287 GeV, respectively, that subsequently decay into a pair of SM-like Higgs bosons. This is also the reason for the enhancement in T2, where the masses of the non-SM-like Higgs bosons amount to 794 and 798 GeV. The reason for the much smaller enhancement in $hh$ production in T2 compared to T1 is the overall heavier Higgs spectrum in T2. In particular, the intermediate heavy resonances in the T2 scenarios that can produce $hh$ in their decay usually fall into the heavy mass range where the ATLAS and CMS limits on the upper cross section for $(4b)$ production drop rapidly, cf. Aaboud et al. (2018e); Collaboration (2017a). In T1, furthermore, the maximum values of di-Higgs production processes involving $H_{\downarrow}$ can compete with SM Higgs pair production or even largely exceed it. Thus the production of a SM-like Higgs boson and $H_{\downarrow}$ can be larger by a factor 2.5. This final state is interesting as it is clearly a non-SM-like signature where the experiments can use the SM-like Higgs boson to calibrate or “tag” this exotic configuration. This does not apply for T2, however, where due to the experimental constraints, the non-SM-like Higgs bosons are in general heavier than in T1, inducing small di-Higgs production processes due to a much smaller phase space.

Experimental accessibility and exclusion luminosity

In order to assess the experimental accessibility of these cross sections, however, we need to look at their decay products. We therefore applied the narrow width approximation and multiplied the produced Higgs bosons with their branching ratios in various SM final states. The most promising final states for the investigation of Higgs pair production at the LHC are the $(b\bar{b})(\gamma\gamma)$ Baur et al. (2004), the $(b\bar{b})(\tau\tau)$ Baur et al. (2003a); Dolan et al. (2012); Barr et al. (2014) and $(b\bar{b})(b\bar{b})$ Baur et al. (2003a); Ferreira de Lima et al. (2014); Wardrope et al. (2015) final states (for other final states see also Baur et al. (2002, 2003b); Papaefstathiou et al. (2013)). In Fig. 2 we show for all parameter points that pass our applied constraints for the C2HDM T1, the cross section values of SM Higgs pair production in the $(b\bar{b})(\gamma\gamma)$ final state (left) and for $H_{\downarrow}H_{\downarrow}$ production in the $(b\bar{b})(b\bar{b})$ final state (right) normalized to the corresponding SM values as a function of the exclusion luminosity. By the latter we define the luminosity at which this process would be excluded experimentally, based on the extrapolations described in Subsection III.3. First of all we notice that the figures contain parameter points at lower luminosity that should have been excluded by HiggsBounds. The reason why they are there is that HiggsBounds relies on the published experimental results and cannot check for certain signatures that become relevant in BSM Higgs sectors. Thus there exist Higgs spectra with heavy Higgs bosons that dominantly decay into top quark pairs. These would induce exotic four-top final states in heavy Higgs pair production. Such signatures compete, however, with single heavy Higgs production and subsequent decay into a top-quark pair. Applying our rough estimate on the exclusion power of the experiments for this process, based on the $Z^{\prime}$ data, such scenarios are excluded already, although they have been let through by HiggsBounds due to the lack of a dedicated experimental analysis for this. This shows the importance of experimental analyses investigating top pair final states from heavy Higgs production in order to properly assess the exclusion limits for BSM Higgs sectors - with dramatic effects on possible Higgs pair production signatures. While our rough extrapolation excludes about 0.6% of the T1 points for a luminosity of about 36 fb*-1*, the effect is much larger for the T2 sample allowed by HiggsBounds††††††HiggsBounds takes into account data at 36 fb*-1*.. Here about 22% of the points would be excluded. This is because of the overall heavy non-SM-like Higgs bosons in T2 and their prominent decays into top-quark pairs.

As can be inferred from the figures in the C2HDM T1, the production of a SM-like Higgs pair with subsequent decay into $(b\bar{b})(\gamma\gamma)$ can exceed the SM rates by up to a factor 60. This maximum enhancement factor is the same for all final states, as the branching ratios of the SM-like Higgs boson $h$ are almost the same as in the SM. In the following, we will use the quantity

[TABLE]

to classify whether a Higgs boson $X$ has a sizable non-SM branching ratio and decay phenomenology. If $\Sigma_{X}\simeq 1$ then the exotic states can be dominantly discovered in “standard” SM-Higgs-like decay channels, e.g. $X\to b\bar{b}$ or $t\bar{t}$ if the mass of $X$ permits such a decay.

In the $H_{\downarrow}H_{\downarrow}$ final state with both $H_{\downarrow}$’s decaying into bottom quarks the enhancement can even be up to a factor of about 200. The point with the maximum enhancement corresponds to the one quoted in Tab. 4 and the enhancement is due to the large di-Higgs production process of 3.2 pb and a slightly enhanced branching ratio into $b$-quarks as compared to the SM. The same factor is found for the $(b\bar{b})(\tau\bar{\tau})$ final state. Due to a smaller branching ratio into photons, however, the maximum allowed enhancement in the $(b\bar{b})(\gamma\gamma)$ final state only amounts up to a factor of 40. The $H_{\downarrow}$ in this scenario has a mass of $m_{H_{\downarrow}}=131$ GeV, and the mass of $H_{\uparrow}$ is $m_{H_{\uparrow}}=313$ GeV. Its main branching ratios are BR($H_{\uparrow}\to ZH_{\downarrow}$)= 0.53 and BR($H_{\uparrow}\to H_{\downarrow}H_{\downarrow}$) = 0.46. The maximum branching ratios of the charged Higgs boson with a mass of $m_{H^{+}}=312$ GeV are BR($H^{+}\to W^{+}H_{\downarrow}$)=0.65 and BR($H^{+}\to t\bar{b}$)=0.34. With its large di-Higgs production cross section and the large non-SM-like branching ratios, this parameter point is an interesting scenario for studying new physics effects (also beyond the Higgs pair events that we consider here).

All remaining di-Higgs production processes are less promising. Thus the enhancement factor for $hH_{\downarrow}$ production remains below 3 in the $4b$ and $2b2\tau$ final state and below 2 in the $2b2\gamma$ final state. All other final states range below the SM values.

As can already be inferred from the maximum di-Higgs production values in T2, given in Tab. 4 the situation looks much less promising in the C2HDM T2. There are very few points in $hh$ production with subsequent decay into the $(2b)(2\tau)$ and $4b$ final state that exceed the SM rate, and only by a factor of about 2.4. The maximum enhancement found in the $(2b)(2\gamma)$ final state is about 2.4. All other final states lead to smaller rates than in the SM.

From these considerations we can conclude that there are promising di-Higgs signatures with large rates in the C2HDM T1 both for SM-like Higgs pair production but also for final states with non-SM-like Higgs bosons. The exotic Higgs bosons appear in SM-like final states, however, with different kinematic correlations due to different masses. This highlights the need to conduct Higgs pair analyses in a broad range of kinematic possibilities. Furthermore, the strict constraints on T2 scenarios, would exclude the model if di-Higgs signatures much larger than in the SM are found.

III.4.2 NMSSM

In Table 5 we summarize for the NMSSM the maximum di-Higgs production cross section values found in the sample of valid parameter points that are not excluded at a luminosity of 64 fb*-1*. All cross sections are calculated at LO QCD and hence still increase by approximately a factor two when QCD corrections are included. By $A_{\downarrow}$ we denote the lighter of the two pseudoscalar Higgs bosons.

The reason for the large enhancement of $\sigma(gg\to hh)$ is the intermediate resonant production of heavy Higgs bosons $H_{\downarrow}$ and $H_{\uparrow}$ with subsequent decay into a SM-like Higgs pair. The $H_{\downarrow}H_{\downarrow}$ production cross section is so large because of the smallness of the $H_{\downarrow}$ mass, $m_{H_{\downarrow}}=39.52$ GeV. The enhancement in $hA_{\downarrow}$ production is due to the resonant $A_{2}\equiv A_{\uparrow}$ production with subsequent decay into $hA_{\downarrow}$. The huge enhancement in $A_{\downarrow}A_{\downarrow}$ production is on the one hand due to the smallness of the $A_{\downarrow}$ mass of $m_{A_{\downarrow}}=37$ GeV, on the other hand due to the resonant $H_{\uparrow}$ production with subsequent decay into $A_{\downarrow}A_{\downarrow}$ (the resonant $H_{\downarrow}$ production plays a minor role). Searches for relatively low mass states are performed in the $\gamma\gamma$ Sirunyan et al. (2018e) and $\tau\tau$ Sirunyan et al. (2018f) channels, however, with rather limited sensitivity.

Experimental accessibility and exclusion luminosity

In Fig. 3 we show for all parameter points that pass our applied constraints, the NMSSM cross section values of SM Higgs pair production in the $(b\bar{b})(b\bar{b})$ final state (left) and for $A_{\downarrow}A_{\downarrow}$ production in the $(b\bar{b})(b\bar{b})$ final state (right) normalized to the corresponding SM values as a function of the exclusion luminosity.

As can be inferred from Fig. 3 (left), the $4b$ final state rates from SM-like Higgs pair production exceed the SM reference value by less than a factor 10 and only for lower exclusion luminosities. Large enhancement factors are basically limited by the LHC upper limits on heavy resonant scalar production with subsequent decay into a SM-like Higgs pair. The situation looks even less promising in the production of an SM-like Higgs boson together with the lighter of the CP-even non-SM-like Higgs bosons, where only an enhancement factor slightly above 2.3 at most is found. This is the case for high exclusion luminosities beyond 1 ab*-1* so that nevertheless this process might be accessible at high luminosities. The situation is different in the production of $hA_{\downarrow}$. Because the lighter pseudoscalar can be relatively light and decays dominantly into $(b\bar{b})$‡‡‡‡‡‡Note that typical trigger criteria are too selective to directly observe $pp\to A_{\downarrow}\to b\bar{b}$. we can have enhancement factors above 10 up to about 45 in the $4b$ state. This makes it particularly interesting, moreover in view of the exotic final state with two different Higgs masses in di-Higgs production. The enhancement factors can become huge in $A_{\downarrow}A_{\downarrow}$ production, which is mainly due to the lightness of $A_{\downarrow}$. In $4b$ production it can be up to 1000. For larger exclusion luminosities the enhancement factor can still be a factor up to 10, cf. Fig. 3 (right). In the $(b\bar{b})(\gamma\gamma)$ final state the enhancement can be larger than 100 up to about 360.

Figure 4 nicely summarizes the situation of the enhanced di-Higgs cross sections involving very light Higgs bosons. It shows the production of $A_{\downarrow}A_{\downarrow}$ with subsequent decay in the $4b$ final state normalized to the value of the corresponding process for the SM-like di-Higgs production, as a function of the mass of the light pseudoscalar. The color code denotes the exclusion luminosity. For very light masses below 125 GeV (note the gap around 125 GeV is due to our scan procedure) the rates are largely enhanced because of the large di-Higgs production cross sections. With increasing mass the rates decrease. The exclusion luminosities are high for exotic Higgs masses above 125 GeV and below the top-pair threshold. Above the top-pair threshold the exclusion luminosities are much lower due to the exclusion limits in the top-pair final state. For masses below the SM-like Higgs mass, however, there are parameter points where the exclusion luminosities can exceed the 100 fb*-1* and even 1 ab*-1* while still featuring large rates. The reason is that these points are not excluded from single Higgs searches as light Higgs states with dominant decays into $b\bar{b}$ final states are difficult to probe. On the other hand this enhancement combined with the large di-Higgs production cross section implies huge $4b$ final state rates, that may be tested at the high luminosities, but with associated experimental difficulties. This is a nice example of the interplay between difficult single-Higgs searches and large exotic di-Higgs rates, where new physics may be found.

IV Benchmarks and Phenomenology

IV.1 Type 1 Benchmarks

We describe a representative set of benchmarks of the C2HDM T1 model and their associated (exotic) multi-Higgs phenomenology. The input parameters, the derived third neutral Higgs boson mass, the CP-odd admixtures in terms of the squared mixing matrix elements $R_{i3}^{2}$ and the exclusion luminosity ${\cal L}_{\text{excl}}$ are given in Tab. 6. We also give the NLO QCD gluon fusion $hh$ production cross section at $\sqrt{s}=14$ TeV together with its $K$-factor, given by the ratio of the NLO cross section to the LO one. In Table 7 we present the $4b$, $(2b)(2\tau)$ and $(2b)(2\gamma)$ rates from Higgs pair production normalized to the rate expected in the SM from Higgs pairs relevant for the discussion of the various benchmark points.

• •

T1BP1 - Highest_exclusion_lumi: The point with the highest exclusion luminosity in the complete sample.

The exclusion luminosity for this point is found to be 11.5 ab*-1*, i.e. well above the LHC design luminosity even after the high luminosity phase. All di-Higgs cross sections involving non-SM-like Higgs bosons have values below the SM reference point. Altogether this benchmark point behaves very SM-like as expected for such a high exclusion luminosity. The neutral Higgs mass spectrum is relatively degenerate with all masses in the vicinity of the SM-like Higgs boson at approximately 125 GeV. The SM-like Higgs pair production cross section for this point is SM-like; $b\bar{b}\gamma\gamma$ is about 10% below the SM expectation while $b\bar{b}\tau\tau$ and $b\bar{b}b\bar{b}$ are slightly enhanced by 5%, cf. Tab. 7. The cases when a SM Higgs is accompanied by an additional exotic Higgs are around 30% of the SM expectation in $b\bar{b}\tau\tau$. With the latest improvements in hadronic tau tagging Cadamuro (2015); Khachatryan et al. (2015b); Mastrolorenzo (2016) such a signature might be phenomenologically accessible at the LHC with large luminosity.

• •

T1BP2 - HighLMaxHsmHl: Among the points with exclusion luminosities $\geq 1~{}\text{ab}^{-1}$ the point with the maximum cross section $gg\to hH_{\downarrow}$.

We have $\sigma(gg\to hH_{\downarrow})=28.47$ fb at LO. As summarized in Tab. 7, for the final states involving $b$’s, $\tau$’s and $\gamma$’s we find for this di-Higgs final state relative to the production of a SM Higgs pair with subsequent decay into the same final state:

[TABLE]

The neutral Higgs mass spectrum in this case is slightly split while the pair production of the SM-like Higgs bosons follows largely the SM paradigm. The mass of the lighter neutral state is rather close to the SM boson, which allows us to compare the rates with the SM itself. $H_{\downarrow}$ has a decay phenomenology that is SM-like. At $3~{}\text{ab}^{-1}$ we can therefore expect around 26k exotic $b\bar{b}\tau\tau$ events while the $b\bar{b}b\bar{b}$ rate is enhanced by 70% over the SM expectation.

• •

T1BP3 - NonSMHsmHL: Among the points with exclusion luminosities $\geq 1~{}\text{ab}^{-1}$ and $\Sigma_{H_{\downarrow}}\leq 0.5$ (i.e. dominantly non-SM-like decays for $H_{\downarrow}$) the point with the maximum cross section $gg\to hH_{\downarrow}$.

This point will have a highly non-SM decay phenomenology and cascade decays are an interesting avenue to look for such a scenario. All di-Higgs production cross sections involving non-SM-like Higgs bosons lie below 5 fb. While the point has interesting signatures for non-SM-like single Higgs production, exotic di-Higgs production is not a very promising avenue. The cross sections are far below the SM value, as no resonance-enhancement is possible as the Higgs mass values are too close. Additionally the branching ratios in SM-like Higgs states are very small, as decays into non-SM-like final states dominate. We have $\Sigma_{H_{\downarrow}}=0.180$, $\Sigma_{H_{\uparrow}}=0.149$. The SM-like Higgs pair production is consistent with the SM expectation within 10% for this point. Exotic production of $H_{\downarrow}$ together with the SM Higgs leads to around 220 exotic $(b\bar{b})(\tau\tau)$ events and around 1900 $(b\bar{b})(b\bar{b})$ at $3~{}\text{ab}^{-1}$.

• •

T1BP4 - MaxLEnhancedHsmHl : Among the points with $gg\to hH_{\downarrow}\geq 10$ fb the point with the maximum exclusion luminosity. The SM-like Higgs pair production modes fall again within $\sim 10\%$ of the SM expectation. The mass of the lightest neutral Higgs boson of about $120~{}\text{GeV}$ is again reasonably close to the SM-like Higgs to allow a direct comparison of expected rates, which are slightly smaller than the SM. Specifically, the exotic $b\bar{b}\tau\tau$ and $b\bar{b}b\bar{b}$ modes are 35% smaller than what we would expect for the SM mode with $m_{h}\simeq 125~{}\text{GeV}.$

• •

T1BP5 - EnhancedHsmHsm: Among the points with exclusion luminosities $\geq 1$ ab*-1* the point with the maximum cross section $gg\to hh$.

The exotic Higgs bosons lie well above the SM-like Higgs state which opens the possibility of enhancing the SM-like di-Higgs production due to resonant enhancement of intermediate $H_{\downarrow/\uparrow}$ production. We find explicitly

[TABLE]

• •

T1BP6 - EnhancedHlHl: Among the points with exclusion luminosities $\geq 1$ ab*-1* the point with the maximum cross section $gg\to H_{\downarrow}H_{\downarrow}$.

The exclusion luminosity is 2.58 $\text{ab}^{-1}$ and we have a light $H_{\downarrow}$ with a mass just above half the SM-like mass so that the branching ratios of the latter remain in accordance with the LHC data. The di-Higgs production of $H_{\downarrow}H_{\downarrow}$ amounts to 1.249 pb at LO QCD. Comparing to the SM expectation, we have

[TABLE]

due to the suppressed decay $H_{\downarrow}\to\gamma\gamma$. The light Higgs almost exclusively decays into $b$ pairs at a branching fraction of $85\%$. This means that such a state is difficult to observe in single Higgs production as trigger criteria typically remove such events from the busy hadronic LHC environment. There is a possibility to observe this state in its $\tau$ modes ($\simeq 8\%$). This point is a nice example how the $4b$ mode can be an important BSM discovery tool when SM mass correlations are relaxed.

We note that none of these points features a strong first order electroweak phase transition. In general we observe that points with a strong first order phase transition do not lead to enhanced rates in the four-particle final states; the cross section for SM-like di-Higgs production is close to the SM value while other di-Higgs production cross sections are smaller than the SM expectation.

In addition to these points which are all characterized by relatively light exotic states we include benchmarks with heavy neutral exotics. This is achieved by adding the additional requirement $m_{H_{\uparrow}}\geq 400$ GeV to the defining criteria of the benchmark points quoted above. We restrict ourselves to three benchmark points highlighting special features. They are called T1BP1_H, T1BP2_H and T1BP5_H in analogy to the their lighter mass spectrum counterparts T1BP1,2,5. Detailed information on the benchmark points is summarized in Tab. 8.

• •

T1BP1_H - Highest_exclusion_lumi: The exclusion luminosity for this point is found to be 2.46 ab*-1*. All di-Higgs cross sections involving non-SM-like Higgs bosons have values below the SM reference as the Higgs pair production cross section steeply falls for heavy Higgs production. The additional resonant structures that are sourced in the SM-like Higgs pair production amount to an increase above the SM expectation by a factor of $\sim 5.8$ across the standard search channels $4b,(2b)(2\tau),(2b)(2\gamma)$.

• •

T1BP2_H - HighLMaxHsmHl : We have $\sigma(gg\to hH_{\downarrow})=2.34$ fb at LO, which is rather large given the mass of the exotic Higgs. The decay phenomenology of the additional Higgs is completely dominated by decays into top-final states. In this sense the single Higgs production and exotic $hH_{\downarrow}$ production are fully correlated. The exclusion luminosity is $1.8~{}\text{ab}^{-1}$ and results from the extrapolation of the $t\bar{t}$ resonance search. This point, although not relevant for di-Higgs analyses shows how single Higgs measurements in the $t\bar{t}$ channel influence multi-Higgs final states. Such a benchmark could be adopted to further clarify the role of single Higgs measurements for exotic multi-Higgs final states.

• •

T1BP5_H - EnhancedHsmHsm: We have enhanced SM-like di-Higgs production due to resonant enhancement of intermediate $H_{\downarrow}$ and $H_{\uparrow}$ production. We find

[TABLE]

This point hence gives access to SM-like Higgs pair production even for a heavy Higgs spectrum, and is an immediate sign of BSM physics as the di-Higgs cross section is enhanced.

IV.2 Type 2 Benchmarks

As already visible from Tab. 4, the C2HDM T2 model gives rise to less spectacular signatures than the C2HDM T1. We give two representative scenarios below. In general, the spectrum is much heavier than for T1. There is no scenario with the SM-like Higgs boson being $H_{2}$; the SM Higgs is always the lightest state $H_{1}$. The input parameters for these points as well as further relevant information are summarized in Tab. 9.

• •

T2BP1 - Highest_exclusion_lumi: The exclusion luminosity for this point is found to be 2.66 ab*-1*. All di-Higgs cross sections involving non-SM-like Higgs bosons have values below the SM reference value. Altogether this benchmark point behaves very SM-like as expected for such a high exclusion luminosity where $t\bar{t}$ resonance searches become sensitive to this scenario.

• •

T2BP2 - EnhancedHSMHSM: The di-Higgs production into a SM-like Higgs pair is enhanced so that we get

[TABLE]

Despite a lower exclusion luminosity this point behaves also very SM like and barely exceeds the rates into SM-like final states of T2BP1 which has a higher exclusion luminosity.

Overall, it will be difficult to probe the C2HDM T2, which features a heavy Higgs spectrum, in di-Higgs production. This is partly also due to the fact that enhanced SM-like Higgs pair production cross sections are already excluded by the LHC limits on resonant heavy scalar production with decay in a SM Higgs boson pair.

IV.3 NMSSM Benchmarks

Let us finally turn to the NMSSM. The criteria for selecting the benchmark points are as follows.

• •

NMBP1: The point with the largest $4b$ rate from SM-like Higgs boson pair production with an exclusion luminosity above 1 ab*-1*.

• •

NMBP2: The same as BP1 but with an exclusion luminosity beyond 100 fb*-1*.

• •

NMBP3: The point with the largest $4b$ rate from the production of an SM-like Higgs boson and the lighter of the CP-even non-SM-like Higgs bosons, $H_{\downarrow}$, with an exclusion luminosity above 1 ab*-1*.

• •

NMBP4: The point with the largest $4b$ rate from $H_{\downarrow}H_{\downarrow}$ production with an exclusion luminosity above 100 fb*-1*.

We also provide benchmark points for di-Higgs final states involving a light pseudoscalar $A_{\downarrow}$ in the final state:

• •

NMBP5: The point with the largest $4b$ rate from $hA_{\downarrow}$ production with an exclusion luminosity above 1 ab*-1*.

• •

NMBP6: The same as BP5 but with an exclusion luminosity beyond 100 fb*-1*.

• •

NMBP7: The point with the largest $4b$ rate from $A_{\downarrow}A_{\downarrow}$ production and an exclusion luminosity above 1 ab*-1*. It turns out that NMBP7 is identical to NMBP5.

• •

NMBP8: The same definition as for NMBP7 but with an exclusion luminosity beyond 100 fb*-1*.

The input values, the derived Higgs boson masses, the exclusion luminosity and the NLO QCD cross section for $hh$ production with its $K$-factor for the various benchmark points are listed in Tabs. 10 and 12. Note, that we include benchmark points with an exclusion luminosity around 100 fb*-1* in case the rates are much enhanced compared to the SM as in this case a luminosity of 100 fb*-1* might be enough to test this parameter point. As the rates for di-Higgs production involving a heavy scalar or pseudoscalar are low we do not present benchmarks for these cases.

From Tab. 11 we can read off that the $4b$ rates from $hh$ production for an exclusion luminosity above 100 fb*-1* (NMBP2) can be up almost a factor of about 4.8, so that this process will be difficult to be accessed at the lower luminosity. Assuming an exclusion luminosity above 1 ab*-1* (NMBP1) the enhancement compared to the SM rate is around 2.4. Again, larger enhancements in the final state with an SM-like Higgs boson pair are excluded by the limits provided from ATLAS and CMS Collaboration (2017a); Sirunyan et al. (2018c) . In the $hH_{\downarrow}$ final state the enhancement factor is only 2.3 as the di-Higgs production cross section

[TABLE]

is not much larger than in the SM. For $H_{\downarrow}H_{\downarrow}$ production we have, however,

[TABLE]

leading to an enhancement factor of up to 540. The $H_{\downarrow}$ mass is much lower here than in the former case inducing dominant branching ratios into $b\bar{b}$. The large cross section is mainly due to the small $H_{\downarrow}$ mass. The resonant enhancement plays a minor role here. Both benchmark points are special in the sense that the SM-like Higgs boson is not the lightest but the second lightest CP-even Higgs boson in the spectrum. We hence have a light CP-even Higgs boson in these scenarios. The $H_{\downarrow}$ is very singlet-like in both cases and decays with a branching ratio of about 0.9 into $b\bar{b}$.

From Tab. 13 we can read off that the $4b$, $(2b)(2\tau)$ and $(2b)(2\gamma)$ final states from both $hA_{\downarrow}$ and $A_{\downarrow}A_{\downarrow}$ production can be enhanced above the SM rate even for an exclusion luminosity of 1 ab*-1*. For the lower exclusion luminosity the enhancement can be huge, in particular in the 4 fermion final state from $A_{\downarrow}A_{\downarrow}$ production (NMBP8). These enhancements are due to large di-Higgs production cross sections which at NLO QCD amount to

[TABLE]

The enhanced di-Higgs cross section values are on the one hand due to the light pseudoscalar masses and on the other hand due to resonant scalar production for $A_{\downarrow}A_{\downarrow}$ or pseudoscalar production for $hA_{\downarrow}$ production. As already noted in the discussion of Fig. 4 this is an example where new physics may lead to huge measurable effects in Higgs pair production while the single Higgs process, here $A_{\downarrow}$ production in gluon fusion, is difficult to access. This is a prime example that demonstrates that despite the very SM-like nature of the 125 GeV Higgs boson Higgs pair production can be far from being SM-like.

We finally remark that in all NMSSM scenarios the stop masses are quite large, of the order of 1 TeV and larger.

V Conclusions

Multi-Higgs final states are statistically limited at the LHC, but are key processes to gain a precise understanding of the mechanism of electroweak symmetry breaking. Phenomenologically, they are highly correlated with measurements in single Higgs final states. The question of how much additional information can be gained from the investigation of multi-Higgs final states is therefore best addressed using concrete BSM extensions.

Another particularly relevant question when considering di-Higgs final states is whether they could play a key discovery tool for BSM interactions. This could happen at the LHC in situations when single Higgs analyses are simply not competitive due to, e.g., trigger thresholds that might be mitigated in more complex multi-Higgs final states. Large branching ratios of additional scalars into SM-like Higgs bosons serve as an additional avenue to observe resonantly enhanced SM-like Higgs production. In these scenarios, the kinematic correlations are often significantly modified compared to the SM.

In this work we have performed a comprehensive scan over the complex 2HDM and the NMSSM, with a particular emphasis on the expected di-Higgs phenomenology in these models, taking into account a variety of current constraints and future projections. We find that in particular in the C2HDM Type 1 models, the di-Higgs phenomenology can significantly differ from the SM expectation. The differences range from new signatures in SM-like search channels for light Higgs bosons all the way to new resonant structures in di-Higgs final states for larger exotic Higgs masses. In particular, in final states involving light Higgs bosons the multi-fermion final states can be significantly enhanced compared to the SM case. This is also the case in the NMSSM where we can have light scalars or pseudoscalars in the spectrum. In the C2HDM Type 2 models with its heavy Higgs spectrum, the multi-Higgs final states play a more subdominant role as new physics discoveries typically occur in $t\bar{t}$ resonance searches before (exotic) di-Higgs production becomes relevant.

We have distilled our scans into a representative number of benchmark points that not only reflect the phenomenological possibilities that present themselves in non-standard Higgs sectors, but also point to a particular range of phenomenological situations:

Firstly, the expected sensitivity of the $t\bar{t}$ resonance search is a crucial factor in deciding on the relevance of di-Higgs searches. The multi-Higgs signal is typically driven by top-mediated gluon fusion. Therefore, the decay to top final states is directly correlated with a large di-Higgs cross section for resonant production as well as enhancements in the decay to photons for non-resonant production.

Secondly, when exotic Higgs masses fall below the $t\bar{t}$ threshold, di-Higgs final states typically follow the SM decay rates with compressed neutral Higgs masses. This highlights the necessity to achieve a high mass resolution in the standard search channels $(2b)(2b),(2b)(2\tau),(2b)(2\gamma)$, even when SM-mass correlations are abandoned. While di-Higgs production would be enhanced in this instance, providing clear evidence of the presence of BSM interactions, their precise nature would remain elusive to some extent.

Alternatively, additional Higgs exotics can create multiple resonant features leading to a large enhancement of the total SM Higgs pair production rate. The extrapolated signal-strength constraints locate viable candidates for enhanced di-Higgs production in $b\bar{b}$ final states, which are difficult to access experimentally in single Higgs production. Here di-Higgs production can play a significant role as a discovery tool for BSM interactions due to smaller backgrounds and better kinematical handles.

Thirdly, relatively light Higgs bosons with significant branching ratios can lead to a strong enhancement in multi-fermion final states. Such signatures are already studied by the experimental collaborations. Our results indicate the importance of these analyses in the future.

While we have specifically focused on di-Higgs production, it is clear that these scenarios can have interesting non-SM signatures that can be exploited to observe or constrain a certain parameter point in a more targeted, yet parameter point-dependent way. We will leave this for future work.

Acknowledgements.

We want to thank Jonas Wittbrodt for providing us with the C2HDM data sample. P.B. acknowledges the support by the DFG-funded Doctoral School Karlsruhe School of Elementary and Astroparticle Physics: Science and Technology (KSETA). S.D. is supported by the U.S. Department of Energy under Grant Contract de-sc0012704. C.E. is supported by the IPPP Associateship scheme and by the UK Science and Technology Facilities Council (STFC) under grant ST/P000746/1.