Probing a Scalar Singlet-Catalyzed Electroweak Phase Transition with Resonant Di-Higgs Production in the $4b$ Channel

TL;DR

This paper assesses the potential of the 14 TeV HL-LHC to detect a heavy Higgs boson decaying into two SM-like Higgs bosons in the 4b final state within a scalar singlet extended Standard Model, focusing on electroweak phase transition implications.

Contribution

It introduces an optimized analysis using boosted decision trees for resonant di-Higgs production in the 4b channel, enhancing detection prospects for a strong first order electroweak phase transition.

Findings

Potential discovery of heavy Higgs boson below 500 GeV with 3 ab$^{-1}$ luminosity.

Enhanced sensitivity of the 4b channel over other decay channels for masses above 500 GeV.

Compatibility with previous projections for other di-Higgs final states.

Abstract

We investigate the prospective reach of the 14 TeV HL-LHC for resonant production of a heavy Higgs boson that decays to two SM-like Higgs bosons in the $4b$ final state in the scalar singlet extended Standard Model. We focus on the reach for choices of parameters yielding a strong first order electroweak phase transition. The event selection follows the $4b$ analysis by the ATLAS Collaboration, enhanced with the use of a boosted decision tree method to optimize the discrimination between signal and background events. The output of the multivariate discriminant is used directly in the statistical analysis. The prospective reach of the $4b$ channel is compatible with previous projections for the $bb\gamma\gamma$ and $4\tau$ channels for heavy Higgs boson mass $m_2$ below 500 GeV and superior to these channels for $m_2 > 500$ GeV. With 3 ab$^{-1}$ of integrated luminosity, it is possible to discover the heavy Higgs boson in the $4b$ channel for $m_2< 500$ GeV in regions of parameter space yielding a strong first order electroweak phase transition and satisfying all other phenomenological constraints.

Tables

Table 1. Table 1: Values of the various xSM independent and dependent parameters for each of the benchmark values consistent with a SFOEWPT chosen to maximize the σ ​ ( p ​ p → h 2 ) × BR ​ ( h 2 → h 1 ​ h 1 ) 𝜎 → 𝑝 𝑝 subscript ℎ 2 BR → subscript ℎ 2 subscript ℎ 1 subscript ℎ 1 \sigma(pp\to h_{2})\times\mathrm{BR}(h_{2}\to h_{1}h_{1}) value at the 14 TeV LHC.

Benchmark	$\cos \theta$	$m_2$	${\Gamma }_{h_2}$	$x_0$	$\lambda$	$a_1$	$a_2$	$b_3$	$b_4$	${\lambda }_{111}$	${\lambda }_{211}$	$\sigma$	BR
		(GeV)	(GeV)	(GeV)		(GeV)		(GeV)		(GeV)	(GeV)	(pb)
B1	0.974	327	0.929	60.9	0.17	-490	2.65	-361	0.52	45	62.2	0.56	0.33
B2	0.980	362	1.15	59.6	0.17	-568	3.26	-397	0.78	44.4	76.4	0.48	0.40
B3	0.983	415	1.59	54.6	0.17	-642	3.80	-214	0.16	44.9	82.5	0.36	0.33
B4	0.984	455	2.08	47.4	0.18	-707	4.63	-607	0.85	46.7	93.5	0.26	0.31
B5	0.986	511	2.44	40.7	0.18	-744	5.17	-618	0.82	46.6	91.9	0.15	0.24
B6	0.988	563	2.92	40.5	0.19	-844	5.85	-151	0.08	47.1	104	0.087	0.23
B7	0.992	604	2.82	36.4	0.18	-898	7.36	-424	0.28	45.6	119	0.045	0.30
B8	0.994	662	2.97	32.9	0.17	-976	8.98	-542	0.53	44.9	132	0.023	0.33
B9	0.993	714	3.27	29.2	0.18	-941	8.28	497	0.38	44.7	112	0.017	0.20
B10	0.996	767	2.83	24.5	0.17	-920	9.87	575	0.41	42.2	114	0.0082	0.22
B11	0.994	840	4.03	21.7	0.19	-988	9.22	356	0.83	43.9	83.8	0.0068	0.079

Table 2. Table 2: Values of the various xSM independent and dependent parameters for each of the benchmark values consistent with a SFOEWPT chosen to minimize the σ ​ ( p ​ p → h 2 ) × BR ​ ( h 2 → h 1 ​ h 1 ) 𝜎 → 𝑝 𝑝 subscript ℎ 2 BR → subscript ℎ 2 subscript ℎ 1 subscript ℎ 1 \sigma(pp\to h_{2})\times\mathrm{BR}(h_{2}\to h_{1}h_{1}) at the 14 TeV LHC.

Benchmark	$\cos \theta$	$m_2$	${\Gamma }_{h_2}$	$x_0$	$\lambda$	$a_1$	$a_2$	$b_3$	$b_4$	${\lambda }_{111}$	${\lambda }_{211}$	$\sigma$	BR
			(GeV)	(GeV)		(GeV)		(GeV)		(GeV)	(GeV)	(fb)
BM1	0.9999	329	0.00593	111	0.13	-812	3.61	-99.8	0.35	31.8	7.30	1.1	0.71
BM2	0.9995	363	0.0549	80.6	0.13	-699	4.16	-91.5	0.57	32.2	21.6	8.2	0.68
BM3	0.9803	419	1.32	234	0.18	-981	1.56	0.417	0.96	39.0	17.5	6.9	0.018
BM4	0.9997	463	0.0864	56.9	0.13	-763	6.35	113	0.73	32.2	27.4	3.0	0.63
BM5	0.9994	545	0.278	50.3	0.13	-949	8.64	152	0.57	33.0	51.6	2.9	0.62
BM6	0.9991	563	0.459	33.0	0.13	-716	9.25	-448	0.96	33.7	66.8	3.7	0.62
BM7	0.9836	609	4.03	34.2	0.22	-822	4.53	-183	0.57	47.8	45.2	2.2	0.030
BM8	0.9870	676	4.48	30.3	0.22	-931	5.96	-680	0.43	48.4	55.2	1.3	0.037
BM9	0.9904	729	4.22	27.3	0.21	-909	6.15	603	0.94	45.7	61.0	0.78	0.045
BM10	0.9945	792	3.36	22.2	0.18	-936	9.47	-848	0.66	43.5	92.4	0.77	0.12
BM11	0.9944	841	3.95	21.2	0.19	-955	8.69	684	0.53	43.3	73.4	0.28	0.062

Table 3. Table 3: Cross section, K factor, and acceptance times efficiency for the different sources of background at the 13 TeV LHC. The expected event yields predicted by the simulation used in this analysis are compared with the expected yields from the ATLAS analysis for 3.2 fb -1 at s = 13 𝑠 13 \sqrt{s}=13  TeV. The cross-section values in the second column include the K factors listed in the third column.

Backgrounds	${\sigma }_{\mathrm{parton}}\left(\mathrm{pb}\right)$	K factor	Efficiency	Expected yield	ATLAS Ref. Aaboud:2016xco
$4b$	287.24	1.72 (NLO QCD) Alwall:2014hca	4.02$\times {10}^{-5}$	37	43
$t\overline{t}$	72	1.60 (N3LO QCD) Muselli:2015kba	1.87$\times {10}^{-5}$	4.0	4.3

Table 4. Table 4: Cross section, K factor, acceptance times efficiency, and estimated event yields for the different sources of background at the 14 TeV LHC before BDT selection. The cross-section values in the second column include the K factors listed in the third column. Uncertainties in the cross-section values are discussed in the text.

Backgrounds	${\sigma }_{\mathrm{parton}}^{\mathrm{NLO}}\left(\mathrm{pb}\right)$	K factor	Efficiency	Expected yield
$4b$	${130}_{-24\%}^{+28\%}$	1.4 (NLO QCD) Czakon:2015qwa	2.99$\times {10}^{-2}$	$1.17\times {10}^7$
$t\overline{t}$	${110}_{-5.8\%}^{+3.8\%}$	2.03 (N2LO +N2LL QCD) Czakon:2013goa	5.58$\times {10}^{-3}$	$1.84\times {10}^6$

Table 5. Table 5: Measurements of the single Higgs boson cross section by the ATLAS Collaboration relative to the SM prediction for different production mechanisms as used in the global fit.

	$\gamma \gamma$	$\tau \tau$	$W{W}^{\ast }$	$Z{Z}^{\ast }$	$bb$
$ggH$	${0.81}_{-0.18}^{+0.19}$Aaboud:2018xdt	${1.02}_{-0.55}^{+0.63}$Aaboud:2018pen	${1.10}_{-0.20}^{+0.21}$Aaboud:2018jqu	${1.11}_{-0.23}^{+0.25}$ Aaboud:2017vzb	N.A.
$VBF$	${2.0}_{-0.5}^{+0.6}$Aaboud:2018xdt	${1.18}_{-0.54}^{+0.60}$Aaboud:2018pen	${0.62}_{-0.35}^{+0.36}$Aaboud:2018jqu	${4.0}_{-1.5}^{+1.7}$ Aaboud:2017vzb	N.A.
$VH$	${0.7}_{-0.8}^{+0.9}$Aaboud:2018xdt	N.A.	N.A.	N.A.	${1.08}_{-0.43}^{+0.47}\left(WH\right)$ ${1.2}_{-0.31}^{+0.33}\left(ZH\right)$Aaboud:2018zhk
$ttH$	${1.39}_{-0.42}^{+0.48}$Aaboud:2018urx	N.A.	N.A.	N.A.	${0.79}_{-0.60}^{+0.61}$Aaboud:2018urx

Equations

$$Y_B=(8.59\pm 0.11)\times 10^{-11}.$$

$$V(H,S)=\text{texttimes}-{\mu }^2\left(H^{†}H\right)+\lambda {\left(H^{†}H\right)}^2+\frac{a_1}{2}\left(H^{†}H\right)S$$

$$\text{texttimes}+\frac{a_2}{2}\left(H^{†}H\right)S^2+\frac{b_2}{2}{S}^2+\frac{b_3}{3}{S}^3+\frac{b_4}{4}{S}^4,$$

$$\displaystyle{\rm det}\left(\begin{array}[]{c c}\partial^{2}V/(\partial s^{2})&\partial^{2}V/(\partial s\partial h)\\ \partial^{2}V/(\partial h\partial s)&\partial^{2}V/(\partial h^{2})\end{array}\right)>0.$$

$$\frac{1}{4}\left((b_4+\lambda -a_2){\cos }^4\alpha +(a_2-2b_4){\cos }^2\alpha +b_4\right).$$

$${\frac{dV}{dh}|}_{h=0,s=0}=0,{\frac{dV}{ds}|}_{h=0,s=0}=0,$$

$$\text{texttimes}{\mu }^2=\lambda v_0^2+\left(a_1+a_2{x}_0\right)\frac{x_0}{2},$$

$$\text{texttimes}{b}_2=\text{texttimes}-b_3{x}_0-b_4{x}_0^2-\frac{a_1{v}_0^2}{4x_0}-\frac{a_2{v}_0^2}{2}.$$

$$\text{texttimes}{b}_3{x}_0+2b_4{x}_0^2-\frac{a_1{v}_0^2}{4x_0}-\frac{(a_1+2a_2{x}_0{)}^2}{8\lambda }>0.$$

$$\left|\frac{a_1}{2}\right|,\left|\frac{a_2}{2}\right|,\left|\frac{b_4}{4}\right|<4\pi .$$

$$m_h^2\equiv \text{texttimes}\frac{d^{2}V}{d{h}^2}=2\lambda v_0^2,$$

$$m_s^2\equiv \text{texttimes}\frac{d^{2}V}{d{s}^2}=b_3{x}_0+2b_4{x}_0^2-\frac{a_1{v}_0^2}{4x_0},$$

$$m_{hs}^2\equiv \text{texttimes}\frac{d^{2}V}{dhds}=\left(a_1+2a_2{x}_0\right)\frac{v_0}{2}.$$

$$m_{2,1}^2=\text{texttimes}\frac{m_h^2+m_s^2\pm \left|m_h^2-m_s^2\right|\sqrt{1+\text{texttimes}{\left(\frac{4m_{hs}^2}{m_h^2-m_s^2}\right)}^2}}{2},$$

$$\left(\begin{array}{c}{h}_1\\ h_2\end{array}\right)=\left(\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right)\left(\begin{array}{c}h\\ s\end{array}\right),$$

$$\sin 2\theta =\text{texttimes}\frac{2m_{hs}^2}{\left(m_1^2-m_2^2\right)}=\text{texttimes}\frac{\left(a_1+2a_2{x}_0\right)v_0}{\left(m_1^2-m_2^2\right)}.$$

$$g_{h_{1}XX}=\cos \theta g_{hXX}^{\mathrm{SM}},g_{h_{2}XX}=\sin \theta g_{hXX}^{\mathrm{SM}},$$

$${\mu }_{h_1\to XX}=\frac{{\sigma }_{h_1}\cdot \text{BR}}{{\sigma }_{h_1}^{\mathrm{SM}}\cdot {\text{BR}}^{\mathrm{SM}}}={\cos }^2\theta ,$$

$${\lambda }_{211}=\frac{1}{4}[(a_1+2a_2{x}_0){\cos }^3\theta +4v_0(a_2-3\lambda ){\cos }^2\theta \sin \theta$$

$$+(a_1+2a_2{x}_0-2b_3-6b_4{x}_0)\cos \theta {\sin }^2\theta -2a_2{v}_0{\sin }^3\theta ].$$

$${\Gamma }_{h_2\to h_1{h}_1}=\frac{{\lambda }_{211}^2\sqrt{1-\text{texttimes}4m_1^2/m_2^2}}{8\pi m_2},$$

$${\Gamma }_{h_2}={\sin }^2\theta {\Gamma }^{\mathrm{SM}}(m_2)+{\Gamma }_{h_2\to h_1{h}_1},$$

$${\mu }_{h_2\to XX}={\sin }^2\theta \left(\frac{{\sin }^2\theta {\Gamma }^{\mathrm{SM}}(m_2)}{{\Gamma }_{h_2}}\right),$$

$${\sigma }_{h_1{h}_1}={\sigma }^{\mathrm{SM}}(m_2)\times s_{\theta }^2\frac{{\Gamma }_{h_2\to h_1{h}_1}}{s_{\theta }^2{\Gamma }^{\mathrm{SM}}(m_2)+{\Gamma }_{h_2\to h_1{h}_1}},$$

$$\Delta \mathcal{O}=(c_{\theta }^2-1){\mathcal{O}}^{\mathrm{SM}}(m_1)+s_{\theta }^2{\mathcal{O}}^{\mathrm{SM}}(m_2)=s_{\theta }^2\left[{\mathcal{O}}^{\mathrm{SM}}(m_2)-{\mathcal{O}}^{\mathrm{SM}}(m_1)\right].$$

$$\begin{array}[]{c}\Delta S\equiv S-S_{\mathrm{SM}}=0.06\pm 0.09\\ \Delta T\equiv T-T_{\mathrm{SM}}=0.10\pm 0.07\end{array}\quad\quad\quad\quad\rho_{ij}=\left(\begin{array}[]{cc}1&0.91\\ 0.91&1\end{array}\right),$$

$$\Delta {\chi }^2(m_2,c_{\theta })=\sum _{i,j}\left[\Delta {\mathcal{O}}_i(m_2,c_{\theta })-\Delta {\mathcal{O}}_i^0\right]({\sigma }^2{)}_{ij}^{-1}\left(\Delta {\mathcal{O}}_j(m_2,c_{\theta })-\Delta {\mathcal{O}}_j^0\right),$$

$$\bar{v}(T)/\sqrt{2}=\bar{\varphi }(T)\cos \alpha (T),\bar{x}(T)=\bar{\varphi }(T)\sin \alpha (T),$$

$$\cos \alpha (T_c)\frac{\bar{\varphi }(T_c)}{T_c}\gtrsim 1.$$

$$a_1/\text{TeV},b_3/\text{TeV}\in [-1,1],x_0/\text{TeV}\in [0,1],b_4,\lambda \in [0,1],$$

$$p_T^{\mathrm{lead}}>\{\begin{array}{ll}400\mathrm{GeV}& \mathrm{if}{m}_{4j}>910\mathrm{GeV},\\ 200\mathrm{GeV}& \mathrm{if}{m}_{4j}<600\mathrm{GeV},\\ 0.65m_{4j}-190\mathrm{GeV}& \mathrm{otherwise},\end{array}$$

$$p_T^{\mathrm{subl}}>\{\begin{array}{ll}260\mathrm{GeV}& \mathrm{if}{m}_{4j}>990\mathrm{GeV},\\ 150\mathrm{GeV}& \mathrm{if}{m}_{4j}<520\mathrm{GeV},\\ 0.23m_{4j}+30\mathrm{GeV}& \mathrm{otherwise},\end{array}$$

$$|\Delta {\eta }_{\mathrm{dijets}}|<\{\begin{array}{ll}1.0& \mathrm{if}{m}_{4j}<820\mathrm{GeV},\\ 1.6\times 10^{-3}{m}_{4j}-0.28& \mathrm{otherwise}.\end{array}$$

Full text

I Introduction

After the discovery of the Higgs boson at the Large Hadron Collider (LHC) Aad:2012tfa ; Chatrchyan:2012xdj , understanding the details of electroweak symmetry-breaking (EWSB) in the context of the thermal history of the universe remains an important challenge for particle physics. In particular, it is possible that EWSB was accompanied by generation of the cosmic baryon asymmetry if new physics beyond the Standard Model (BSM) was active during that era. The Planck measurement of this asymmetry, characterized by the baryon-to-entropy density ratio $Y_{B}=n_{b}/s$, gives Ade:2013zuv :

[TABLE]

Explaining the origin and magnitude of $Y_{B}$ is a key problem for BSM scenarios. Electroweak baryogenesis (EWBG) is one of the appealing possibilities, in part due to its linking of $Y_{B}$ to EWSB and in part due to its testability in current and near future experiments. Three “Sakharov conditions” Sakharov:1967dj need to be satisfied for a successful EWBG: baryon number (B) violation, C and CP violation, and departure from thermal equilibrium (through a strong first order electroweak phase transition) or a breakdown of CPT symmetry. In the Standard Model (SM), the first condition – baryon number – violation can be induced by the process of electroweak sphalerons. However, the CP violation in the SM is too feeble, and the EWSB transition is a crossover transition given the observed SM Higgs mass $m_{h}\sim 125$ GeV Aoki:1999fi ; Csikor:1998eu ; Laine:1998jb ; Gurtler:1997hr ; Kajantie:1996mn . Therefore, the minimal SM cannot generate a successful strong first order electroweak phase transition (SFOWEPT). On the other hand, if new scalars exist in addition to the SM Higgs doublet, their interactions with the SM Higgs doublet may catalyze a SFOEWPT, thereby providing the necessary conditions for successful EWBG.111 New CP-violating interactions would also be required, a topic we do not treat further here.

In this paper, we focus on the singlet extension to the SM, the xSM, which is proven to be able to give a SFOEWPT Profumo:2007wc ; Espinosa:2011ax . In the xSM, after EWSB, the gauge eigenstates of the singlet scalar and the SM Higgs doublet mix with each other to form the mass eigenstates $h_{1}$ (SM-like) and $h_{2}$ (singlet-like). Further, we restrict our study to searching for a signal of the on-shell production of the heavy singlet-like Higgs $h_{2}$ decaying into two SM-like Higgs $h_{1}$ (i.e. $m_{2}>2m_{1}$), because the regions of parameter space that can generate SFOEWPT simultaneously tend to enhance the $h_{2}h_{1}h_{1}$ tri-linear couplings Profumo:2007wc ; Espinosa:2011ax ; No:2013wsa . Currently, the ATLAS and CMS experiments are searching for a resonant di-Higgs signal through different Higgs decay final states: $4b$ Aaboud:2018knk ; Sirunyan:2018zkk , $bbWW^{*}$ or $bbZZ^{*}$ Aaboud:2018zhh ; Sirunyan:2017guj , $bb\tau\tau$ Aaboud:2018sfw ; Sirunyan:2017djm , $bb\gamma\gamma$ Aaboud:2018ftw ; Sirunyan:2018iwt , $WW^{*}WW^{*}$ Aaboud:2018ksn , and $\gamma\gamma WW^{*}$ Aaboud:2018ewm . Thus far, no significant excess over SM backgrounds has been observed. On the theoretical side, several studies have been performed in the parameter regions that are viable for SFOEWPT. The singlet-like $h_{2}$ with a relatively light mass ($\sim$270 GeV) can be discovered in the $bb\tau\tau$ final state at the 14 TeV LHC with a luminosity of 100 fb*-1* No:2013wsa . In the $bb\gamma\gamma$ and $4\tau$ final states, a discovery is possible for $m_{2}$ up to 500 GeV at the 14 TeV high-luminosity LHC (HL-LHC) with a luminosity of 3 ab*-1* Kotwal:2016tex . In the $bbWW^{*}$ final state, a resonant signal can be discovered for $m_{2}$ in the range between 350 GeV and 600 GeV at the 13 TeV LHC with a luminosity of 3 ab*-1* Huang:2017jws .

In this paper, we study the prospective discovery/exclusion in the $4b$ final state at the 14 TeV HL-LHC with a luminosity of 3 ab*-1*. To that end, we first identify 22 benchmark points with $m_{2}\in[300,850]$ GeV that produce the maximal and minimal di-Higgs signal rate $\sigma_{h_{2}}\times{\rm BR}(h_{2}\to h_{1}h_{1})$ in consecutive 50 GeV intervals. The selected benchmark points satisfy all the current phenomenological constraints from the Higgs signal rate and electroweak precision data, and also satisfy the theoretical constraints from vacuum stability, perturbativity, and a SFOEWPT. We perform a full simulation of signal and background processes with the MadGraph5 parton level event generator Alwall:2014hca using PYTHIA6 Sjostrand:2006za to simulate the parton shower and the DELPHES3 fast detector simulation deFavereau:2013fsa . Further, we use the TMVA package TMVA to implement the Boosted Decision Tree (BDT) algorithm to optimize the event selection, finally obtaining the signal significance from the BDT score distributions of background and signal events.

Based on this analysis and the results shown in Fig. 4 below, we arrive at the following conclusions:

• •

For singlet-like Higgs masses below 500 GeV, the significance of the $4b$ final state is competitive with the $bb\gamma\gamma$ and $4\tau$ final states, and it is possible to make a discovery at the 14 TeV HL-LHC with a luminosity of 3 ab*-1* for some portions of the SFOEWPT-viable parameter space.

• •

For singlet-like Higgs masses above 500 GeV, the significance of the $4b$ final state is higher than in the $bb\gamma\gamma$ and $4\tau$ final states but somewhat below recent projections for the $bbWW^{\ast}$ final state.

• •

With the results of the benchmark models that produce minimal di-Higgs signal rate, we found that it is impossible to exclude (at the 95% confidence level) all portions of parameter space consistent with a SFOEWPT and present phenomenological constraints at the HL-LHC.

The discussion of our analysis leading to these conclusions is organized as follows: Sec. II introduces the xSM framework and describes both theoretical and phenomenological constraints. In Sec. III, we describe the requirements for a SFOEWPT and the parameter scan. In Sec. IV, we discuss the simulation and analysis of the $4b$ signal and background in detail and also present prospects for the 14 TeV HL-LHC. Section V is dedicated to the conclusions. In the Appendix, we perform a global analysis of ATLAS Run 2 single Higgs measurements and present the distributions of the kinematic variables used in the BDT analysis.

II The xSM

II.1 The Model

The most general, renormalizable scalar potential in the xSM model is given by:

[TABLE]

where $S$ is the real singlet and $H$ is the SM Higgs doublet. When $S$ obtains a vacuum expectation value (vev, see below), the $a_{1}$ and $a_{2}$ parameters induce mixing between the singlet scalar and the SM Higgs doublet, thereby providing a portal for the singlet scalar to interact with other SM particles. A $\mathbb{Z}_{2}$ symmetry is present in the absence of $a_{1}$ and $b_{3}$ terms, a necessary condition for $S$ to be a viable dark matter candidate. In what follows, however, we retain both parameters in our study as they play an important role in the strength of the electroweak phase transition (EWPT) and also in the di-Higgs signal rate at collider experiments.

After EWSB, $H\to(v_{0}+h)/\sqrt{2}$ with $v_{0}=246$ GeV, and $S\to x_{0}+s$ where $x_{0}$ is the vev for $S$ without loss of generality. The stability of the scalar potential requires the quartic coefficients along all the directions in the field space to be positive. This translates into a requirement of a positive Hessian determinant of the potential with respect to fields $s$ and $h$:

[TABLE]

This leads the bounds $\lambda>0$, $b_{4}>0$ and $a_{2}>-2\sqrt{\lambda b_{4}}$. Another way to obtain these bounds is by parameterizing ($h$,$s$) as ($r\cos\alpha$, $r\sin\alpha$) in the field space, and we are able to extract the quartic coefficients of $r$ along the $\alpha$ direction in the field space:

[TABLE]

Requiring the above expression be larger than zero for any value of $\cos\alpha$ also leads to the same conditions.

Utilizing the minimization conditions,

[TABLE]

one can express two potential parameters in Eq. (II.1) in terms of the vevs and other parameters:

[TABLE]

Two additional conditions need to be satisfied for ($v_{0}$, $x_{0}$) to be a stable minimum. One of them is that ($v_{0}$, $x_{0}$) minimizes the potential locally, implying that:

[TABLE]

Also, this minimum point should be a global minimum, a requirement that we impose numerically.

As for the perturbativity consideration, we have the following naïve requirements on the quartic couplings:

[TABLE]

However, as discussed in Sec. III when scanning over the parameter space for benchmark points we implement more stringent bounds on those parameters compared with the above requirements. One may refer to Refs. Robens:2015gla ; Robens:2016xkb ; Gonderinger:2009jp for more details about the perturbativity bound in the xSM.

Now we obtain the elements of the mass-squared matrix by:

[TABLE]

After the diagonalization of the above mass matrix, the physical masses of two neutral scalars can be expressed as:

[TABLE]

with $m_{2}>m_{1}$ by construction. The mass eigenstates and gauge eigenstates are related by a rotation matrix:

[TABLE]

where $h_{1}$ is the SM-like Higgs boson with $m_{1}=125$ GeV, and $h_{2}$ is identified as the singlet-like mass eigenstate. The mixing angle $\theta$ can be expressed in terms of the vevs, physical masses and potential parameters:

[TABLE]

From Eq. (13), one can observe that the couplings of $h_{1}$ and $h_{2}$ to the SM vector bosons and fermions are rescaled with respect to their SM Higgs couplings:

[TABLE]

where $XX$ represents final states consisting of pairs of SM vector bosons or fermions. In this case, all the signal rates associated with the single Higgs measurements are rescaled by the mixing angle only:

[TABLE]

where $\sigma_{h_{1}}$ and BR are the production cross section and branching ratio in the xSM, and the quantities with the superscript SM are the corresponding values in the SM. In the xSM for $m_{2}>m_{1}$, we have BR = $\text{BR}^{\mathrm{SM}}$ because the partial width of each decay mode is rescaled by $\cos^{2}\theta$ and there is no new decay channel appearing.

In order to investigate the di-Higgs production, we also require the tri-Higgs couplings. The one relevant for the resonant di-Higgs production is $\lambda_{211}$:

[TABLE]

In this work, we focus on the situation where $m_{2}>2m_{1}$ such that a resonant production of $h_{2}$ and a subsequent decay to $h_{1}h_{1}$ is allowed. Therefore, we are able to calculate the partial width $\Gamma_{h_{2}\to h_{1}h_{1}}$:

[TABLE]

and the total width of $h_{2}$:

[TABLE]

where $\Gamma^{\mathrm{SM}}(m_{2})$ represents the total width of the SM Higgs boson with a mass of $m_{2}$, which is taken from Ref. Heinemeyer:2013tqa . The signal rate for $pp\to h_{2}\to XX$ normalized to the SM value is given by:

[TABLE]

which will be used to constrain the parameter space in the next section. The production cross section for the process $pp\to h_{2}\to h_{1}h_{1}$ can also be calculated:

[TABLE]

where $s_{\theta}\equiv\sin\theta$ and, for future reference, $c_{\theta}\equiv\cos\theta$.

II.2 Phenomenological Constraints on the Model Parameters

The mixing angle $\theta$ between the singlet and the SM Higgs doublet in the xSM is constrained by measurements of the single SM-like Higgs signal strengths. We obtain a 95% C.L. upper limit on $\sin^{2}\theta$ of 0.131 by performing a global fit with current ATLAS Run 2 single Higgs measurements as discussed in Appendix A.1.

The LHC searches for the heavy neutral Higgs boson also provide constraints on the parameter space. Here, we take into account the existing limits on both the $h_{2}\to VV$ Sirunyan:2018qlb ; Aaboud:2018bun ; Aad:2015kna ; Aad:2015ipg ; Chatrchyan:2013yoa ; Khachatryan:2015cwa and the $h_{2}\to h_{1}h_{1}$ decays, where $h_{1}h_{1}$ decay into 4$b$ Aaboud:2018knk ; Sirunyan:2018zkk , $b\bar{b}\gamma\gamma$ Aaboud:2018ftw ; Sirunyan:2018iwt , or $b\bar{b}\tau\tau$ Aaboud:2018sfw ; Sirunyan:2017djm . The constraints on the ($m_{2}$, $c_{\theta}$) plane can be found in our previous work Huang:2017jws . We will also guarantee each benchmark point in the parameter scan in the next section satisfies all the limits mentioned above.

Finally, we discuss the constraints from electroweak precision observables (EWPO). The mixing between the singlet scalar and the SM Higgs doublet induces modifications of the oblique parameters $S$, $T$, and $U$ with respect to their SM values. From Eq. (13), the deviation in oblique parameters ${\cal O}$, denoted by $\Delta{\cal O}$, can be expressed in terms of the SM Higgs contribution to that parameter, ${\cal O}^{\mathrm{SM}}(m)$ Peskin:1991sw ; Hagiwara:1994pw and the mixing angle $\theta$, where $m$ is either $m_{1}$ or $m_{2}$:

[TABLE]

In the xSM, the parameter $U=0$ is a good approximation; we therefore focus only on the deviations in the $S$ and $T$ parameters, which we take from the Gfitter group Baak:2014ora :

[TABLE]

where $\rho_{ij}$ is the covariance matrix in the ($S$,$T$) plane. Again, we will impose the criteria in the parameter scan in the next section such that for each benchmark point, $\Delta\chi^{2}(m_{2},c_{\theta})$ defined below is less than 5.99, which corresponds to deviations of $S$ and $T$ parameters within 95% C.L.:

[TABLE]

where the $\Delta\mathcal{O}^{0}_{i}$ denote the central values in Eq. (23) and $(\sigma^{2})_{ij}\equiv\sigma_{i}\rho_{ij}\sigma_{j}$, with $\sigma_{i}$ being the error in $S$ or $T$ as indicated in Eq. (23). One can observe from Fig. 1 in Ref Huang:2017jws that in general the upper limit for $\sin^{2}\theta$ extracted from EWPO is more stringent than the bound obtained from the Higgs global fit, with a limit changing from 0.12 for $m_{2}=250$ GeV to 0.04 for $m_{2}=950$ GeV.

III Electroweak Phase Transition and Benchmarks for Di-Higgs Production

The character of EWPT is understood in terms of the finite-temperature effective potential, $V^{T\neq 0}_{eff}$ However, the fact that the standard derivation of $V^{T\neq 0}_{eff}$ suffers from gauge dependence is well known which is discussed in depth in Ref. Quiros:1994dr . Here we employ a high-temperature expansion to restore the gauge independence in our analysis (see Ref. Profumo:2014opa for details). In such a case, we include in our finite temperature effective potential the $T=0$ tree level potential and the gauge-independent thermal mass corrections to $V^{T\neq 0}_{eff}$, which are crucial to restore electroweak symmetry at high temperature. In this limit, the $a_{1}$ and $b_{3}$ parameters will generate a tree-level barrier between the broken and unbroken electroweak phases, thereby allowing for a first-order EWPT. We also note that the presence of the $a_{2}$ term may also strengthen the first order transition, as discussed in Ref. Profumo:2007wc . In the high-temperature limit, we follow Ref. Pietroni:1992in ; Profumo:2007wc and write the $T$-dependent, gauge-independent (indicated by the presence of a bar) vevs in a cylindrical coordinate representation as

[TABLE]

with $\bar{v}(T=0)=v_{0}$ and $\bar{x}(T=0)=x_{0}$. The critical temperature $T_{c}$ is defined as the temperature at which the broken and unbroken phases are degenerate, i.e. $V^{T\neq 0}_{eff}(\phi,\alpha\neq\pi/2,T_{c})=V^{T\neq 0}_{eff}(\phi,\alpha=\pi/2,T_{c})$. Once the critical temperature is found, one is able to evaluate the quenching effect of the sphaleron transitions in the broken electroweak phase (see, e.g., Ref. Morrissey:2012db ), which is related to the energy of the electroweak sphaleron that is proportional to the vev of SU(2)L doublet $\bar{v}(T)$. A first-order EWPT is strong when the quenching effect is sufficiently large, and the criterion is approximately given by:

[TABLE]

To select the benchmarks parameter points for the collider simulation, we perform a scan over the parameters $a_{1}$, $b_{3}$, $x_{0}$, $b_{4}$, and $\lambda$ within the following ranges:

[TABLE]

while the remaining parameters are fixed from the input values of $v_{0}=246$ GeV and $m_{h}=125$ GeV. Our lower bound on quartic couplings $b_{4}$ and $\lambda$ guarantees tree-level vacuum stability. We also require a naïve perturbativity bound on the Higgs portal coupling $a_{2}/2\lesssim 5$. For each set of randomly chosen parameters, we calculate $c_{\theta}$, $m_{2}$, and $\lambda_{211}$, and only keep the points that satisfy all the phenomenological constraints mentioned in the previous section (Higgs signal rate, LHC search for heavy Higgs $h_{2}$, and EWPO). We then pass these sets of parameters into the CosmoTransitions package Wainwright:2011kj and numerically evaluate all the quantities related to the EWPT, such as critical temperature, sphaleron energy, tunneling rate into the electroweak symmetry-broken phase, using as an input the xSM finite temperature effective potential in the high-temperature limit. Finally, we only keep the sets of parameters that satisfy the strong first-order EWPT criterion defined above and also have a sufficient tunneling rate to prevent the universe from remaining in a metastable vacuum.

From the randomly chosen parameters satisfying the foregoing requirements, we identify benchmark points with maximum and minimum signal rate $\sigma(pp\to h_{2})\times{\rm BR}(h_{2}\to h_{1}h_{1})$ from 11 consecutive $h_{2}$ mass windows of width 50 GeV ranging from 300 to 850 GeV. The upper bound of $m_{2}=850$ GeV is obtained by the observation that we did not find a choice of parameters for $m_{2}$ larger than 850 GeV that give a SFOEWPT, even though our scan in $m_{2}$ reaches one TeV. We list all the benchmark points in Tables 1 and  2. We would like to mention that the benchmark points B3 and B4 for maximum signal rate in Table 1 has already been excluded by the CMS $h_{2}\to ZZ$ search Sirunyan:2018qlb , but we retain them here to make contact with the results of previous studies for comparison. In contrast, the new ATLAS and CMS limits on resonant di-Higgs production in the $b{\bar{b}}\tau\tau$ channelAaboud:2018sfw ; Sirunyan:2017djm do not yet appear to constrain the SFOEWPT-viable parameter space.

IV 4b Final State Analysis

IV.1 Reproduction of 13 TeV LHC results

For the signal process, the $h_{2}$ mass is varied from 300 GeV to 1500 GeV in steps of 100 GeV. For the background processes, we generate $pp\to 4b$ and $pp\to t\bar{t}$ with top quarks decaying hadronically. We follow the ATLAS resolved analysis in Ref. Aaboud:2016xco , and reproduce the signal efficiency and background distributions in Figs. 1 and 2, respectively. Parton level signal and background events are generated with MG5_AMC@NLOv2.4.3 Alwall:2014hca and the NNPDF2.3QED LO set of parton distribution functions Ball:2013hta . For the $4b$ QCD background, we generate events with the process p p > b b b$\sim$ b$\sim$, while all other parton level cuts are set to the Madgraph default values. For the $t\bar{t}$ background, we generate events with the process p p > t t$\sim$,(t > b c s$\sim$),(t$\sim$ > b$\sim$ c$\sim$ s) plus one additional jet with jet matching. The xqcut in the run card is set to 20 GeV, and other cuts are kept at the default settings. The events are interfaced with PYTHIA6 Sjostrand:2006za for parton showering, fragmentation and hadronization. DELPHES3 deFavereau:2013fsa is used to simulate the detector response. The default CMS DELPHES card is used rather than the ATLAS DELPHES card as it better approximates the $b$-tagging and jet reconstruction performance. Jets are constructed using the anti-$k_{t}$ clustering algorithm with a radius parameter $R$ set to 0.4, and the efficiency for a $b$-quark-initiated jet to pass the $b$-tagging requirements is parameterized as a function of the jet transverse momentum $p_{T}$ in a manner corresponding to an average 70% efficiency working point described in Ref. Chatrchyan:2012jua . (This is the default setting in the DELPHES CMS card).

The selection criteria for the ATLAS analysis are as follows:

• •

Events must have at least four $b$-tagged jets with $p_{T}>$ 40 GeV and $|\eta|<2.5$. If the number of $b$-tagged jets is greater than four, the four jets with the highest $p_{T}$ are selected to reconstruct two dijet systems in each event.

• •

Two dijet systems are formed using the selected $b$-tagged jets. The two jets in each dijet system are required to have $\Delta R<1.5$ and the transverse momentum of the leading (subleading) dijet system must be greater than 200 (150) GeV.

• •

The leading and subleading dijet systems must satisfy the following set of requirements depending on the reconstructed invariant mass ($m_{\mathrm{4j}}$) of the four selected $b$-tagged jets:

[TABLE]

[TABLE]

[TABLE]

• •

To reduce the $t\bar{t}$ background, we impose a “$t\bar{t}$ veto” as follows. A set of $W$-boson candidates is formed by combining one of the $b$-tagged jets in the dijet system with any extra jet in the event that satisfies $p_{T}>30$ GeV and $|\eta|<2.5$ as well as $\Delta R<1.5$ relative to the dijet system. Top-quark candidates are then formed by combining the dijet system with each of the extra jets that are selected. An event is vetoed if the invariant masses of the $W$-boson ($m_{W}$) and top-quark ($m_{t}$) candidates satisfy the following condition for any possible choice of extra jet and $b$-tagged jet from either of the dijet systems in the event:

[TABLE]

• •

Finally, the signal region is defined by the following requirement on the invariant masses of the leading and subleading dijet systems forming the two Higgs boson candidates:

[TABLE]

The central values for $m_{\mathrm{2j}}^{\mathrm{lead}}$ and $m_{\mathrm{2j}}^{\mathrm{subl}}$ in the above equation are somewhat lower than in the ATLAS analysis Aaboud:2016xco to account for differences in the treatment of jets in DELPHES compared to the ATLAS simulation.

The acceptance times efficiency values for signal events with $m_{2}$ ranging from 500 to 1000 GeV are compared with the ATLAS results in Fig. 1. Overall, the signal region efficiencies obtained in this analysis agree well with those from Ref. Aaboud:2016xco .

The background event yields in the signal region are summarized in Table 3. In addition to the yields from $4b$ and $t\bar{t}$ production, the contribution from $bbcc$ production with the $c$-quark jets passing the $b$-tagging requirements is estimated assuming that the kinematic distributions of jets in $bbcc$ events are similar to those of $4b$ events:

[TABLE]

where $N_{4b}$ is the estimated number of QCD $4b$ events, $\sigma_{bbcc}$ and $\sigma_{4b}$ are parton level cross sections for $bbcc$ and $4b$ processes, $\epsilon^{\mathrm{tag}}_{c}$ and $\epsilon^{\mathrm{tag}}_{b}$ are the $b$-tagging efficiencies for $c$-quark and $b$-quark jets taken to be 0.2 Aad:2015ydr and 0.7, respectively. The expected number of $bbcc$ background events is two with a luminosity of 3.2 fb*-1*, i.e. about 5% of the total background. The distribution of the reconstructed $4b$ invariant mass for background events is shown in Fig. 2. This distribution is obtained after rescaling the momenta of the dijet systems such that their invariant masses are equal to 125 GeV. Good agreement is observed between the background estimate from this analysis and the ATLAS results from Ref. Aaboud:2016xco .

IV.2 Predictions for 14 TeV HL-LHC

After demonstrating that we are able to reproduce the ATLAS results obtained at the 13 TeV LHC, we evaluate the prospects at the 14 TeV HL-LHC with a modified analysis. The event selection is modified according to Ref. ATLAS:2016ixk :

• •

Events are required to have at least four $b$-tagged jets with $p_{T}>30$ GeV and $|\eta|<2.5$.

• •

Dijet systems are formed such that the separation $\Delta R_{\mathrm{jj}}$ between the two jets satisfies the following requirements:

[TABLE]

[TABLE]

• •

If more than one pair of dijet systems satisfies this constraint, the pair with the smallest variable $D_{h_{1}h_{1}}$ is selected with

[TABLE]

In order to optimize the separation between signal and background events, the analysis in this paper relies on a BDT trained on half of the simulated signal and background events and validated with the other half. Separate training is performed for each benchmark point studied. The kinematic quantities included in the training of the BDT are $p_{T}^{\rm lead},\ p_{T}^{\rm subl},\Delta R_{\mathrm{jj}}^{\mathrm{lead}},\Delta R_{\mathrm{jj}}^{\mathrm{subl}},\Delta R_{h_{1}h_{1}},\Delta\phi_{h_{1}h_{1}},\Delta\eta_{h_{1}h_{1}},m^{\rm lead}_{\rm 2j},m^{\rm subl}_{\rm 2j},X_{h_{1}h_{1}},m_{\mathrm{4j}}$, where the variable $X_{h_{1}h_{1}}$ is defined as

[TABLE]

Among those variables, $\Delta R_{\mathrm{jj}}^{\mathrm{lead}}$, $\Delta R_{\mathrm{jj}}^{\mathrm{subl}}$, and $m_{\mathrm{4j}}$ are consistently ranked high in terms of discrimination power for all benchmark points. To derive the optimal sensitivity, BDT score distributions for signal and background events are rebinned such that each bin contributes the maximum $S/\sqrt{B}$ ($S$ and $B$ are the numbers of signal and background events in that bin), starting from the bin with the highest BDT score where the signal contributes the most. This rebinning also requires a minimum of ten background events per bin to minimize the impact of statistical fluctuations. As an illustration, the rebinned BDT score distributions for two benchmark points are shown in Fig. 3.

The production cross sections and the efficiencies of backgrounds before the BDT selection are summarized in Table 4. The uncertainties for these backgrounds correspond to the theoretical uncertainties associated with variations in the renormalization and factorization scales and uncertainties from the parton distribution functions. In the case of the $pp\to 4b$ background process, a parton-level requirement of $\Delta R_{bb}>0.3$ is imposed in the generation of events with MadGraph5 to allow the use of the NLO cross-section calculation at $\sqrt{s}=14$ TeV from Ref. Czakon:2015qwa . Such a requirement is consistent with the radius parameter $R=0.4$ used in the anti-$k_{t}$ algorithm as this sets an effective lower bound of $\Delta R=0.4$ between two jets.

To evaluate the sensitivity to di-Higgs scalar resonances, we calculate the ${\rm CL}_{b}$ value from the rebinned BDT score distributions with the profile likelihood method using the asymptotic formula described in Refs. Cowan:2013pha ; Cowan:2010js . The quantity $1-{\rm CL}_{b}$, which represents the probability that the background-only model yields an observed number of events at least as large as the expectation for the signal plus background model, is then translated into the corresponding $N_{\sigma}$ Gaussian significance. As a test of the statistical analysis, it was verified that the 95% upper limit on the cross section as a function of resonance mass derived from our emulation of the 13 TeV ATLAS analysis (discussed in Sec. IV.1) agrees with the results from Ref. Aaboud:2016xco within 10% for $h_{2}$ masses up to 750 GeV and within 20% up to 850 GeV. The slight deviation at higher mass may be due to the use of the asymptotic formula which is known to produce upper limits that are too aggressive for the low number of expected events at high mass with only 3.2 fb*-1* of luminosity.

The significance $N_{\sigma}$ as a function of resonance mass is shown in Fig. 4, where the upper and lower boundaries of the band correspond to the influence of uncertainties in the production cross sections for the $4b$ and $t\bar{t}$ backgrounds as given in Table 4. The two boundaries are obtained by coherently changing the number of events for the two backgrounds by the 1$\sigma$ uncertainties listed in Table 4, computing the CLb according to the method mentioned above, and then converting CLb to $N_{\sigma}$. One can observe that, with 3 ab*-1* of integrated luminosity at the 14 TeV HL-LHC, the benchmark points with maximum signal rate up to $m_{2}=500$ GeV can be discovered with $N_{\sigma}>5$. If the future HL-LHC experiments do not observe a signal, then one can exclude the maximum signal rate benchmark points up to $m_{2}=680$ GeV at 95% C.L.

The significance is compared to that obtained with the same method for the $bb\gamma\gamma$ and $4\tau$ channels at the 14 TeV HL-LHC Kotwal:2016tex and for the $bbWW$ channel at the 13 TeV LHC Huang:2017jws in Fig. 5. We only compare the benchmark points from BM3 to BM11 because the first two BM points are different from those in Ref. Kotwal:2016tex . We find that for a heavy Higgs mass $m_{2}$ less than 500 GeV, the $bb\gamma\gamma$ channel is the most sensitive channel in the search for a resonant di-Higgs signal. Moreover, the $4b$ channel is competitive with the $bb\gamma\gamma$ channel, which could serve as a complementary check if a signal is observed in the $bb\gamma\gamma$ channel. However, for $m_{2}$ larger than 500 GeV the $4b$ channel provides better sensitivity than the $bb\gamma\gamma$ or $4\tau$ channels but not as good as the $bbWW^{*}$ channel Huang:2017jws . We note, however, that the analysis given in Ref. Huang:2017jws employs a novel Heavy Mass Estimator (HME) and assumptions that the systematic uncertainties will be improved compared to those quoted in the recent CMS $bbWW^{*}$ analysis Sirunyan:2017guj that did not implement the HME. These differences may account for the stronger projected limits given in Ref. Huang:2017jws than one would infer by rescaling the results in Ref. Sirunyan:2017guj by the improved statistics expected for the HL-LHC Luca . We also note that Ref. Huang:2017jws assumes an ATLAS-CMS combination, thereby doubling the number of events. We do not make such an assumption in the present study.

V Conclusion

Investigating the thermal history of EWSB is important for determining whether or not the cosmic matter-antimatter asymmetry was generated through EWBG. Monte Carlo simulations indicate that the EWSB transition is cross over in the minimal SM, given the the observed Higgs mass. In this context of a SM-only universe, EWBG cannot occur. However, introducing new scalar degrees of freedom can change the behavior of thermal effective potential and make the SFOEWPT possible during the EWSB era. Adding a real scalar singlet is one of the simplest ways to extend the SM – yielding the xSM – and realize this possibility. Previous studies have demonstrated the existence of a strong correlation between an enhanced coupling of a heavy singlet-like scalar to a SM-like di-Higgs pair and the occurrence of a SFOEWPT in the xSM parameter space. Therefore, there exists strong motivation to search for resonant production of heavy singlet-like scalar that decays to SM-like di-Higgs state as a probe of SFOEWPT in the xSM.

In this paper, we focused on the possibility of discovering at the HL-LHC a resonant gluon fusion production of the heavy singlet-like scalar in the xSM that decays into a pair of SM-like Higgs with a four $b$-quark final state. The four $b$-quark final state is a promising channel, given its large branching ratio, but it also suffers from a significant QCD background. In analyzing this process, we first validated our simulation against the ATLAS 13 TeV cut-based analysis, then implemented the BDT, a multi-variable analysis method to help to classify signal and background events for the HL-LHC. We selected 11 benchmark points for both maximum and minimum di-Higgs signal rates that yield a SFOEWPT and that satisfy all the theoretical and phenomenological bounds for a heavy singlet-like scalar mass in successive 50 GeV energy bins ranging from 300 to 850 GeV. We then analyzed the signal significance for the 14 TeV HL-LHC with a luminosity of 3 ab*-1*. We also compared the results with earlier projections for the $bb\gamma\gamma$ and $4\tau$ channels and find that for the mass of the singlet-like scalar larger than 500 GeV, the significance for the $4b$ channel is superior to both of these other channels. For heavy singlet-like scalar mass less than 500 GeV, the significance for the $4b$ state with maximum signal rate can be larger than 5. This significance is comparable to that of the $bb\gamma\gamma$ final state, and is somewhat better than that projected for the $4\tau$ final state. While our projection for the reach using the $4b$ channel is somewhat below that for the $bbWW^{\ast}$ channel as analyzed in Ref. Huang:2017jws , the latter work utilized a new Heavy Mass Estimator and assumptions about future reductions in systematic uncertainties that await validation with new data. Thus, inclusion of the $4b$ channel in a comprehensive search strategy that also includes the $bb\gamma\gamma$, $4\tau$, and $bbWW^{\ast}$ channels is strongly motivated. In terms of exclusion, we find that for the future 14 TeV HL-LHC, one can exclude the mass of a heavy singlet-like scalar up to around 680 GeV for the benchmark points with maximum signal rate. However, a signal in the case of the minimum signal rate benchmark points is far from being excluded. Therefore in this sense, a future 100 TeV $pp$ collider may be required to fully exclude the possibility of generating SFOEWPT in the xSM.

Acknowledgement

HL and MJRM were supported in part under U.S. Department of Energy contract DE-SC0011095. HL is supported by the National Science Foundation of China under Grants No. 11875003. SW is supported in part under U.S. Department of Energy contract DE-SC0010004.

Appendix A Appendix

A.1 Single Higgs global fit with ATLAS Run 2 Results

We use the public code Lilith Bernon:2015hsa to implement a global fit with following observables:

[TABLE]

Here, $X$ represents the production mode (e.g. gluon fusion, vector boson fusion etc.) and $Y$ represents the final states into which the SM-like Higgs decays. We list the data and the corresponding references we use from ATLAS Run 2 results in Table 5. The statistical $\chi^{2}$ is:

[TABLE]

where $C^{-1}$ is the inverse of the covariance matrix $cov[\mu^{obs}_{i},\mu^{obs}_{j}]$. In principle we need to know the whole $n\times n$ covariance matrix ($n$ is the number of observables we use in the global fit) to compute the $\chi^{2}$, but doing so is impossible, as the full matrix is not provided by the experimental collaboration. Instead, we ignore the off-diagonal part in the covariance matrix and approximate the $\chi^{2}$ as:

[TABLE]

where $\sigma_{X,Y}$ denotes the $1\sigma$ uncertainty for the given observable. The treatment of asymmetric uncertainties is discussed in the Lilith documentation Bernon:2015hsa .

Figure 6 shows the result of our global fit in the $\Delta\chi^{2}$ vs $\sin^{2}\theta$ plane, the $\Delta\chi^{2}$ is defined by:

[TABLE]

where $\chi^{2}_{min}$ is the minimum value of $\chi^{2}$ in the scan. This translates into a 95% C.L. upper bound on $\sin^{2}\theta<$0.131, given by the requirement that $\Delta\chi^{2}<$3.841.

Appendix B Distributions of BDT variables

We plot the signal and background distributions of kinematic variables used in the BDT analysis here. The signal is taken to be the benchmark point B7 in Table 1.