Resonant Di-Higgs Production in the b{\bar{b}}WW Channel: Probing the Electroweak Phase Transition at the LHC


We analyze the prospects for resonant di-Higgs production searches at the LHC in the (, ) channel, as a probe of the nature of the electroweak phase transition in Higgs portal extensions of the Standard Model. In order to maximize the sensitivity in this final state, we develop a new algorithm for the reconstruction of the invariant mass in the presence of neutrinos from the decays, building from a technique developed for the reconstruction of resonances decaying to pairs. We show that resonant di-Higgs production in the channel could be a competitive probe of the electroweak phase transition already with the datasets to be collected by the CMS and ATLAS experiments in Run-2 of the LHC. The increase in sensitivity with larger amounts of data accumulated during the High Luminosity LHC phase can be sufficient to enable a potential discovery of the resonant di-Higgs production in this channel.

a]T. Huang, b,c]J. M. No, a]L. Pernié, d,e]M. Ramsey-Musolf, a]A. Safonov, f]M. Spannowsky, d]P. Winslow \affiliation[a]Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843, USA \affiliation[b]Department of Physics, King’s College London, Strand, WC2R 2LS London, UK \affiliation[c]Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH, UK \affiliation[d]Physics Department, University of Massachusetts Amherst, Amherst, MA 01003, USA \affiliation[e]Kellogg Radiation Laboratory, California Institute of Technology, Pasadena, CA 91125, USA \affiliation[f]Institute of Particle Physics Phenomenology, Physics Department, Durham University, Durham DH1 3LE, UK \ \ \ \ \ \ \

1 Motivation

With the discovery of the Higgs boson at the Large Hadron Collider (LHC) [1, 2], exploring the thermal history associated with electroweak symmetry-breaking (EWSB) has taken on heightened interest. In the Standard Model (SM), EWSB in the early Universe occurs through a crossover transition. In contrast, beyond the Standard Model (BSM) scenarios may lead to a bona fide electroweak phase transition (EWPT). If such a transition occurred and was both first order and sufficiently strong, it could have provided the conditions needed for generating the observed cosmic matter-antimatter asymmetry.

Electroweak baryogenesis (EWBG) (for recent reviews, see [3, 4]) is one of the most widely-studied and experimentally testable scenarios for explaining the origin of the cosmic matter-antimatter asymmetry, characterized by the baryon-to-entropy density ratio as (most precisely) measured by Planck [5]


Successful baryogenesis requires three ingredients in the particle physics of the early Universe, the so-called “Sakharov criteria” [6]: (i) baryon number (B)-violation; (ii) C- and CP-violation; (iii) departure from thermal equilibrium or a breakdown of CPT invariance. The SM contains the requisite B-violation in the guise of electroweak sphalerons, but it fails with regard to the last two criteria. CP-violation in the SM, via the CKM mixing matrix, is too feeble. In the minimal SM, the maximum Higgs mass for a first order EWPT is GeV, as confirmed by a variety of theoretical Monte Carlo simulations [7, 8, 9, 10, 11], while for the observed GeV EWSB occurred through a cross-over phase transition in the early Universe, which would not provide for the necessary out-of-equilibrium conditions.

In contrast, if the observed Higgs boson resides within an extended scalar sector, the nature and properties of the EWPT could differ significantly from those of the SM. In that case, the Universe could have undergone a strong first order EWPT even for a SM-like Higgs boson of mass GeV. The additional scalar degrees of freedom can alter the finite temperature effective potential to make such a transition possible. The simplest realization of this possibility involves the extension of the SM Higgs sector by a single real scalar singlet , the xSM [12, 13, 14, 15, 16]. While the xSM in and of itself is unlikely to be realized in Nature, it embodies the phase transition dynamics associated with more complete models that contain a gauge singlet, such as the next-to-minimal supersymmetric Standard Model (NMSSM) [17, 18], without introducing the complications associated with the other degrees of freedom present in these models. As such, it provides a framework for exploring generic features of singlet-driven phase transitions and the corresponding low-energy phenomenology. To enable a strong first order EWPT in the xSM, the coupling(s) between the new scalar and the SM Higgs doublet need to be sizable (though still perturbative). In general, the xSM gives rise to two neutral scalars with masses that are mixtures of the singlet and neutral component of the doublet. The corresponding phenomenological consequences include reduced SM-like Higgs boson signal strengths [12, 13, 16], deviations of the trilinear Higgs self-coupling from its SM value [19, 16] and resonant di-Higgs production[20, 21].

In this work we focus on resonant di-Higgs production at the LHC: , where () denotes the SM-like (singlet-like) neutral scalars, for GeV. Within the SM, di-Higgs production is non-resonant, and search strategies have been proposed in  [22, 23, 24, 25],  [26, 27],  [26, 28] and  [29, 30, 31] final states, all found to be very challenging due to the smallness of the non-resonant di-Higgs cross section [32, 33, 34, 35, 36]. Hence, ongoing di-Higgs searches by ATLAS and CMS are looking beyond the SM paradigm, also focusing on resonance-enhanced production mechanisms [37, 38, 39, 40].

In this context, two key issues need to be addressed. First, it is important to assess the LHC reach into the viable parameter space for a strong first order EWPT. There have been initial studies in this context in the final state [21], which found that discovery at the LHC may be possible with 100 fb of integrated luminosity for relatively light masses, but a comprehensive analysis has not yet been achieved. Second, in order to achieve maximal LHC sensitivity to the xSM, it is crucial to determine the degree to which different di-Higgs final states provide complementary probes of in different regions of the possible mass range. Here we consider the prospects for LHC discovery/exclusion of resonant di-Higgs production in the channel, which has been initially studied in [41] for low masses ( GeV). We cover the entire mass range GeV TeV, focusing on what is possible to achieve with LHC Run 2. We assess the LHC potential for probing the strong first order EWPT parameter space by defining a set of twelve benchmark xSM parameter choices (corresponding to twelve mass windows in the range GeV), each of them giving the maximum resonant di-Higgs production cross section [] consistent with a strong first order EWPT within its mass window.

For the analysis, we use a Multi Variate Analysis (MVA) discriminator in order to efficiently discriminate the signal from production, the most important SM background (particularly as increases). Conventional experimental techniques do not allow full reconstruction of the resonance mass, which results in diminished discrimination against the leading background. To improve the sensitivity of the analysis, we have deployed a novel technique for the reconstruction of the invariant mass in the process in the presence of neutrinos from the W decays. The proposed method builds on the Missing Mass Calculator (MMC) technique developed for the reconstruction of resonances decaying to pairs [42]. The new technique provides an estimator for using a likelihood constructed over the solutions of the kinematically underconstrained system and, for brevity, is referred to as the Heavy Mass Estimator (HME) in the remainder of the paper.

We find that considering the () final state (and assuming an eventual combination between the CMS and ATLAS experiments) allows to probe into the strong first order EWPT parameter space, defined by [], up to GeV with 300 fb of integrated luminosity, making this channel a promising avenue for analysis during LHC Run 2 and the high luminosity phase of the LHC (HL-LHC).

Our work is organized as follows: Section 2 gives an overview of the xSM, including a summary of current phenomenological constraints and the choice of first order EWPT-viable benchmark points. In Section 3 we discuss the analysis of the channel in detail, including signal and background event generation, object reconstruction, and the algorithm for probabilistic reconstruction of the event kinematics. In Section 4 we apply this analysis to determine the LHC Run 2 and HL-LHC reach. Finally, in Section 5 we offer a summary and outlook.

2 The xSM

2.1 The Model

We consider the most general form for the xSM scalar potential that depends on a Higgs doublet, , and real singlet, (see e.g. [12, 13, 14]):


The and parameters constitute the Higgs portal, providing the only connection between the SM and the singlet scalar . We note that in the absence of and the scalar self-interaction , the potential (2.1) has a symmetry that remains exact if the singlet field does not develop a vacuum expectation value (vev). We however retain both parameters in the current study, as they play a leading role in the EWPT as well as in di-Higgs production at colliders.

Boundedness of the scalar potential from below requires positivity of the quartic coefficients along all directions in field space. Along the () direction, this leads to the bound (), while along an arbitrary direction this implies . After EWSB, with GeV, and we allow for a possible vev for , i.e. , where is taken to be positive without any loss of generality (provided that and can take either sign).

The minimization conditions allow for two of the parameters in (2.1) to be expressed in terms of the vevs and other parameters. For convenience, we choose


For viable EWSB, two conditions must be satisfied: () has to be a stable minimum, which requires


Furthermore, the electroweak minimum must be the absolute minimum, which we impose numerically. After EWSB, the Higgs portal parameters and the singlet vev induce a mixing between the states and . The mass-squared matrix entries are


with the corresponding eigenvalues given by


with by construction. The mass eigenstates are given by


where we identify the more -like state with the Higgs boson observed at the LHC [1, 2] by setting GeV, and where is a singlet-like mass eigenstate. The mixing angle is defined as


By virtue of (7), the couplings of and to SM vector bosons and fermions are universally rescaled w.r.t. the SM Higgs couplings,


with representing a SM final state different from , and , , . In addition to these couplings, the tri-scalar interactions will play an important role in the following discussion of di-Higgs production. Of particular interest are the interactions and , which follow from (2.1) after EWSB, with


2.2 Current Phenomenological Constraints

The singlet-doublet mixing is constrained by measurements of Higgs signal strengths, since all the signal rates associated with Higgs measurements get rescaled by . Currently, the limit from LHC Run 1 data is at 95% C.L. [43]. In addition, ATLAS and CMS searches for a heavy SM-like Higgs boson provide a probe of . For , which we focus on in this work, the decay mode is kinematically allowed, with a partial width given by


Defining as the SM Higgs width evaluated at (as given e.g. in [44]), the total width for the boson is given by


and the signal strength (normalized to the SM value for ) for is


By means of (11)-(13), we can then express the production cross sections (with gauge bosons) and as


with being the SM Higgs LHC production cross section and the SM Higgs branching fraction into for , and given by


We note that heavy Higgs searches in other final states (e.g. ) are much more challenging than the ones discussed above, and are disregarded in what follows. The ATLAS and CMS Collaborations have performed searches for heavy Higgs bosons both in the  [45, 46, 47, 48] and with  [37, 38] and ,  [39, 40]. In Fig. 1 we show the 95% C.L. limits from these searches in the () plane for increasing values of , using (14).

It is important to note that the limits shown in Fig. 1 are derived from present, publicly available (published) analyses, and better analysis techniques can make the search channel more competitive relative to the one, as compared to present bounds. We also stress that a discovery in any one channel would not allow by itself to individually measure and . This highlights the need to explore various channels and final states, like , to correctly interpret a potential discovery of a new state .

Figure 1: Present 95% C.L. excluded regions from Public ATLAS and CMS (blue) [45, 46, 47, 48] and (red) [37, 38, 39, 40] searches in the () plane, for different values of . The colour gradient for corresponds to values of (from lighter to darker) while the colour gradient for corresponds to values of (from lighter to darker). Also shown are the 95% C.L. lower limit on from Higgs signal strength measurements (horizontal solid-black line) [43] and the 95% C.L. lower limit on from EWPO (dashed-black line), both being independent of .

We now turn to the discussion of the constraints on and from electroweak precision observables (EWPO). The effects of the xSM on EWPO may be accurately characterized by its modification of the oblique parameters , , and w.r.t. the SM. From (7), the shift in an oblique parameter can be written entirely in terms of the SM Higgs contribution to that parameter, (which can be found e.g. in [49, 50]), where is either or . These shifts then take the form


The best-fit values for the shifts and standard deviations from the most recent post-Higgs-discovery electroweak fit to the SM by the Gfitter group [51] (for , which is a very accurate approximation in the xSM) are given by


being the covariance matrix in the plane. We then perform a fit to obtain the 95% C.L. allowed region in the () plane


where denote the central values in (17) and , being the and standard deviation from (17). The exclusion limit from is shown in Fig. 1.

2.3 The Electroweak Phase Transition: Benchmarks for Production

The character of the EWPT is understood in terms of the finite temperature effective potential, (see [52] for a review). It is well-known that the standard derivation of suffers from gauge dependence1 [53], and here we employ a high-temperature expansion to restore gauge-independence to our analysis (see [16] for details). Doing so requires considering the Coleman-Weinberg 1-loop effective potential and retaining only the gauge-independent thermal mass corrections to , which are essential for high-temperature electroweak symmetry restoration. This limit is particularly well-suited to the xSM, which generates the barrier between the broken and unbroken electroweak phases required for a first order EWPT at tree-level via the parameters and in (2.1).

In the high-temperature limit, we follow [54, 12] and write the -dependent, gauge-independent (indicated by the presence of a bar) vevs in a cylindrical coordinate representation as


with and . The critical values and are determined by minimizing , while is defined as the temperature at which the broken and unbroken phases are degenerate: . A strong first order EWPT is defined by a sufficient quenching of the sphaleron transitions in the broken electroweak phase (see e.g. [3] for details). The energy of the electroweak sphaleron is proportional to the SU(2)-breaking energy scale, , and as such the approximate criterion for a strong first order EWPT is then .

With these considerations in mind, we implement the xSM in the high-temperature limit in CosmoTransitions [55] to obtain numerically all above quantities characterizing the EWPT and calculate the finite temperature thermal tunneling rate into the electroweak phase; the latter must be sufficiently fast in order to preclude the possibility of the Universe becoming stuck in a false metastable phase. Taking , , , and as our independent parameters (the remaining two are fixed by the values of and ), we perform a MC scan of the xSM parameter space within the following ranges


where the lower bounds on the quartic couplings and ensure vacuum stability. With our choice of independent parameters, , and are fixed by the parameters of the scan. We impose a naïve perturbativity bound on the Higgs portal coupling [15]. We require compatibility with the various experimental constraints discussed in Section 2.2, and demand a strong first order EWPT as described above, together with a sufficient tunneling rate.

(GeV) (GeV) (GeV) (GeV) (GeV) (GeV (GeV)   (pb)
B1 0.961 258 0.68 307 0.52 -266 0.26 -138 0.26 110 -94.6 1.19 0.50
B2 0.976 341 2.42 257 0.92 -377 0.39 -403 0.77 204 -150 0.59 0.74
B3 0.982 353 2.17 265 0.99 -400 0.45 -378 0.69 226 -144 0.44 0.76
B4 0.983 415 1.59 54.6 0.17 -642 3.80 -214 0.16 44.9 82.5 0.36 0.33
B5 0.984 455 2.08 47.4 0.18 -707 4.63 -607 0.85 46.7 93.5 0.26 0.31
B6 0.986 511 2.44 40.7 0.18 -744 5.17 -618 0.82 46.6 91.9 0.15 0.24
B7 0.988 563 2.92 40.5 0.19 -844 5.85 -151 0.08 47.1 104 0.087 0.23
B8 0.992 604 2.82 36.4 0.18 -898 7.36 -424 0.28 45.6 119 0.045 0.30
B9 0.994 662 2.97 32.9 0.17 -976 8.98 -542 0.53 44.9 132 0.023 0.33
B10 0.993 714 3.27 29.2 0.18 -941 8.28 497 0.38 44.7 112 0.017 0.20
B11 0.996 767 2.83 24.5 0.17 -920 9.87 575 0.41 42.2 114 0.0082 0.22
B12 0.994 840 4.03 21.7 0.19 -988 9.22 356 0.83 43.9 83.8 0.0068 0.079
Table 1: Values of the various xSM independent and dependent parameters for each of the benchmark values chosen to maximize the value at the LHC.

From the results of our scan, we define twelve consecutive mass windows, each 50 GeV wide (starting from the production threshold GeV), which together span the range GeV. The upper bound GeV is determined by the fact that the scan does not yield experimentally viable points compatible with a strong first order EWPT above GeV even though it potentially accepts points up to 1 TeV. Among all viable model points within each mass window we select the points which yield the maximum resonant di-Higgs production cross section [] (depending essentially on the values of and , as discussed in Section 2.2) to define twelve xSM strong first order EWPT motivated benchmarks. This benchmark point set, which we refer to in the following as BM, is presented in Table 1. Searches for resonant di-Higgs production in the channel sensitive to this set of benchmarks will be capable of probing into the strong first order EWPT region. In the remainder of the paper we assess the LHC potential to probe such a strong first order EWPT via resonant di-Higgs production in the final state.

3 Resonant Di-Higgs Production: the Channel

As discussed above, in this work we explore the LHC sensitivity to resonant di-Higgs production in the (, ) final state. By exploiting the two largest branching ratios of the 125 GeV Higgs boson we can retain sensitivity to smaller production cross sections, i.e. larger , and develop dedicated reconstruction approaches to suppress SM backgrounds. We require both bosons to decay leptonically (with ) to suppress the otherwise overwhelming background from QCD multi-jet production. The cancellation of momenta of two neutrinos in the decay does not allow to reconstruct the invariant mass of the heavy resonance, which substantially diminishes the LHC sensitivity to resonant di-Higgs production. To improve the sensitivity of the search, we develop a novel technique, called Heavy Mass Estimator, designed to estimate the most likely invariant mass of the heavy state probabilistically. The technique is conceptually similar to the Missing Mass Calculator (MMC) algorithm, which has previously been applied successfully to the mass reconstruction of resonances decaying into pairs [42, 56].

Throughout the remainder of the paper, we assume an LHC center-of-mass-energy of 13 TeV and an integrated luminosity ranging between 300 fb and 3000 fb, expected to be collected between respectively the end of LHC Run 2 data-taking (foreseen in 2022) and the end of the High Luminosity phase of LHC (foreseen in 2035).

3.1 Monte Carlo Generation and Object Reconstruction

For each of the xSM benchmark points in Table 1 we generate our signal (, ) using Herwig++ [57]. The dominant SM background is top-pair () production2, which has been simulated at next-to-leading order (NLO) accuracy with Powheg [58] and then subsequently processed with Herwig++ for parton showering and hadronization to evaluate experimental sensitivity. For simplicity, we restrict our signal and background Monte Carlo generation to , but the subsequent sensitivity analysis takes into account the would-be contributions from final states with two electrons and one electron and one muon, for which we expect very similar efficiencies.

To evaluate the sensitivity achievable with the LHC data in this channel, we use the CMS detector and performance parameters as a benchmark. Assuming similar performance of the CMS and ATLAS detectors, in the combination we double the luminosity delivered per experiment. The simulation of the CMS detector response is performed using Delphes 3.3.0 [59] and the recommended by CMS input card [59, 60], with all reconstructed physical objects, such as tracks, calorimeter deposits, isolated muons, electrons, jets, and missing transverse energy , used in the data analysis. Multiple proton-proton collisions during the same bunch-crossing (pileup) can have a strong impact on hadronic observables, particularly during the high luminosity LHC runs, and we include in our reconstruction the effect of an average of 40 simultaneous proton-proton interactions. A particle-flow algorithm [59] has been successfully deployed in the CMS experiment and is implemented in Delphes parametrically using the information from the tracking system and the calorimeters. The particle-flow method is designed to reconstruct individual particles arising from collision by combining information from relevant subdetectors, to improve the quality of particle identification and the performance of global event reconstruction. Muons are reconstructed within the detector acceptance . Reconstruction and isolation selections follow the CMS definitions developed for particle flow muons [61] and use the medium working points for both. Jets are reconstructed using the anti-  jet algorithm [62] with cone size R=0.4 by clustering the particle-flow tracks and particle-flow towers and we require . The jet area method [59] in Delphes is applied in jet reconstruction to subtract the pileup contribution. The b-tagging efficiency and mis-identification rates are modelled using the Delphes parametrization of the CSV algorithm [63]. Tagging efficiencies and the mistag rates correspond to about 70% and 1.5% for the medium working point and 85% and 10% for the loose working point, respectively. The total transverse energy of a single event is calculated as the 2D-vector sum of the transverse momentum of all particles reconstructed by the CMS particle-flow algorithm. The missing transverse energy is defined as the opposite of the total transverse energy, and it quantifies the transverse energy carried away from neutrinos.

3.2 Invariant Mass Reconstruction for : Heavy Mass Estimator

The cancellation of momenta of the two or more undetected neutrinos in the final state does not allow the reconstruction of the invariant mass of the heavy scalar (and similarly for one of the 125 GeV scalars ) using experimentally measurable quantities. To improve the analysis sensitivity, the HME technique, a MMC-like probabilistic algorithm [42, 56], can be efficiently implemented for the reconstruction of the mass of . To illustrate the implementation of the HME algorithm, we start with an idealized detector, in which properties of all visible particles are perfectly measured and the missing transverse energy is equal to the negative vector sum of all visible particles. The latter assumes that the missing transverse energy measurement is not affected by pileup.

We note that for the production process considered both 125 GeV states are on-shell, whereas one of the two bosons from the decay is typically off-shell (we use the label for the on-shell , ). With these simplified assumptions, the kinematics of the majority of the signal events satisfies the following:


where GeV, , are the x- and y-components of the missing transverse energy vector, represent the various momentum 4-vectors and , are their x- and y-components. The momentum carried by each neutrino is described by three unknown momentum projections, leading to a total of five equations, one bound and six unknowns.

Figure 2: Left: Global likelihood function (solid blue) computed from HME for a single signal event using the xSM benchmark point 3 (B3). The true value of the mass is marked by the red grid bar. Right: HME distribution for B3 (red circle), B6 (blue square), B9 (magenta up triangle) and (green down triangle). All distributions are normalized to unity.

As seen from equations (21)-(25), four constraints reduce the number of unknowns to two, which we choose as the pseudo-rapidity and azimuthal angle of one neutrino. Assigning random values to these two unknowns would then allow one to scan the parameter space of allowed solutions to build a procedure to integrate over the space of solutions consistent with the experimental measured quantities. We refer to a single generation of the two unknowns as iteration if it respects the bound for the invariant mass of the off-shell W (else such single generation is discarded). Each iteration yields an estimator for the mass of :


Furthermore, as not all pairs of values of the unknowns and are equally likely, generating pairs of these values according to a suitably defined probability density function would increase frequency of the estimated mass being close to the true value. Such probability density function (PDF) can be obtained from a Monte Carlo simulation. For each event, we generate thousands of iterations according to the PDF for and , and for each iteration we store the calculated value of building a probability distribution function for , which we refer to as the HME global likelihood function. In the full implementation of the algorithm, the values of and used in  (23)-(25) are generated according to Gaussian functions to account for the width of Higgs and bosons. The addition of these two variables effectively increases the dimensionality of the space in which the scan is performed to four. The introduction of additional probability density functions to account for realistic resolutions of experimental measurements further increases the dimensionality of the space scanned. One of the most essential additions accounts for the b-jet energy mismeasurements which leads to the invariant mass of the two b-jets being on average lower than . We compute and apply an energy correction extracted from the simulation for the leading b-jet, and use equation (25) to correct the energy of the sub-leading b-jet. This procedure simultaneously improves the missing transverse energy estimation used in  (21)-(22), that is finally smeared according to the detector resolution predicted by Delphes. Fig. 2 (left) provides an illustration of a typical HME global likelihood for a single event, which peaks near the true value of the heavy scalar mass . In this analysis, we use the most probable mass from the likelihood as the estimator of the heavy Higgs mass . Note that selecting the peak position of a single event global likelihood as the estimator is the simplest solution, which one could likely improve upon by utilizing more information on the shape of the likelihood or even using the entire distribution in the analysis. Fig. 2 (right) shows the reconstructed mass for various xSM benchmark scenarios described in Table 1.

3.3 Analysis Selection

The experimental signature consists of two energetic leptons, two energetic b-tagged jets and significant missing transverse energy due to neutrinos. As discussed above, the dominant SM background process is , with a very large production cross section (see e.g. [64]). Other potential SM backgrounds are Drell-Yan, single-top, di-boson and production [65], as well as production of the SM Higgs boson (decaying to ) in association with jets (e.g. in vector boson fusion). However, it has been shown in [41] (see also [65]) that basic selection criteria together with mild kinematic cuts on , and the invariant mass of the and di-lepton systems yield all these other backgrounds to be negligible, while maintaining a high signal efficiency. We therefore can safely disregard them in the present work.

Variable Cut
Table 2: Pre-MVA selection.

Initial event pre-selection is performed as follows: we require the presence of two muons3 with opposite sign and GeV, ; if more than two muons are present in the event, the two oppositely charged muons with largest transverse momentum are selected. In addition, at least two b-tagged jets with GeV and are required. At least one of the two b-jet candidates has to be b-tagged using the CSV algorithm at the medium working point, while the other jet is only required to satisfy at least the loose b-jet requirement. If more than two b-jets satisfy all selection criteria, the two b-jets candidate with the invariant mass closer to GeV are selected. Finally we require the missing transverse energy to be GeV. After event pre-selection, we also perform a set of kinematic cuts (pre-MVA selection) summarized in Table 2, which reject approximately 5% of the signal events (for all signal mass points) and about 40% of events.

Figure 3: Kinematic variables (top left), (top middle), (top right), (bottom left), (bottom middle), (bottom right), with distributions normalized to unity.

In order to optimize the background and signal discrimination we have tested several MVA algorithms (Likelihood, LikelihoodMIX, KNN, MLP, BDT, and BDTD) available from the TMVA package (version 4.1.2) [66], choosing the algorithm most performing in terms of background/signal discrimination as a function of : a BDT for low mass (xSM benchmarks B1 to B7), and a Likelihood method for high mass (xSM benchmarks B8 to B12). The training of the MVA has been done independently for each signal mass point4 considered in the analysis using the discriminating variables , , , , , , , , , , and . The variable is computed by measuring the between each lepton and each jet and selecting the smaller among these values. Kinematic distributions after pre-selection for the first six variables above are shown in Fig. 3 for the xSM signal samples B3, B6 and B9 together with the background. The transverse mass variable is defined as


The variable [67, 68] provides a transverse mass estimate in systems where more than one neutrino is present, treating the lepton and the b-jet as a single object. Fig. 4 shows the and distributions (normalized to unity) after the pre-MVA selection. The discriminating power of the variable in di-Higgs final states, already appreciated in [28], is good for all signal samples where the invariant mass of the heavy resonance is greater than twice the top mass.

Figure 4: (left) and (right) distributions after the pre-MVA selection, with distributions normalized to unity.

Table 3 shows the expected event yield with 300 fb of integrated luminosity for the background and each of the xSM signal benchmarks after pre-selection and pre-MVA selection cuts, as well as the final yield of signal () and background () events after the MVA selection in each benchmark scenario. After the MVA-based selection is applied, the invariant mass of the heavy Higgs scalar is reconstructed for the surviving events using the HME probabilistic technique described in Section 3.2. The HME distribution for signal and background is then used for setting upper limits on the signal production cross section.

Process (pre-selection) (pre-MVA cuts) (MVA) (MVA)
B1 395 383 183 6962
B2 395 385 171 27372
B3 318 310 152 26593
B4 137 134 35 1425
B5 104 102 19 193
B6 52 50 12 95
B7 31 30 10 91
B8 22 21 4.5 28
B9 13 12 3.2 23
B10 5 5 1.2 13
B11 3 3 0.8 10
B12 1 1 0.2 4.5
Background 782721 382836
Table 3: Number of signal () and background () events expected collecting 300 fb of integrated luminosity after applying different stages of selection and prior to the final fit using the HME mass estimator distribution.

3.4 Systematic Uncertainties

For the systematic uncertainties in evaluating signal acceptance we assume, based on previous publications presented by the CMS collaboration, a conservative systematic uncertainty of about 10% [69, 70, 71, 72]. This systematic includes the uncertainty on the integrated luminosity, the uncertainty on the trigger efficiency, the lepton identification and isolation, the uncertainty on the parton distribution functions (PDFs) and on the factorization and renormalization scales.

The systematic uncertainties associated with the precision in the knowledge of the background shape and normalization can significantly affect the sensitivity of the analysis. In the absence of real data, we illustrate a possible way to estimate these uncertainties, which we believe to be conservative, as the experimentalists are likely to deploy more sophisticated approaches to reduce the impact of systematic uncertainty on the sensitivity and we also anticipate improvements in the quality of the description of the background arising from theoretical efforts and Monte Carlo tuning. For the purposes of this study, we define a control region dominated by the background, and use it to compare data and simulation to obtain a scale-factor (SF) for correcting the simulation prediction for the contribution. The control region is designed to contain background events with properties and kinematics as in the main signal region.

Since the kinematics has no strong dependence on the di-jet or di-lepton invariant mass, we define control regions by selecting events with the measured di-jet invariant mass greater than 150 GeV or di-lepton invariant mass greater than 100 GeV (see Fig. 3). Using more than one control region allows to cross-check and validate the scale-factors and, if needed, adjust the uncertainty associated with the scale factors. Once the control region is defined we apply the same kinematic selections and perform the MVA training exactly as it is done for the main analysis. We choose the MVA cut that yields the same background rejection as the cut that has been found to yield optimal sensitivity for the same target mass point in the main analysis. For all mass points, the signal contribution remains negligible in the control regions. Finally, the yield of surviving background events is used to derive the uncertainty in the scale factor (in our case, the same events play the role of both the ”data” and the ”prediction” so the mean value of the scale factor is by definition equal to unity). Following this methodology, we estimate the systematic uncertainty on the knowledge of the background normalization to be 1% for the signal samples B1, B2 and B3; 5% for B4; 10% for B5; 12% for B6, B7, B8, B9, B10; and 15% for B11 and B12. In the final sensitivity estimates, for the three lowest mass points we chose to increase the systematic uncertainty for the background normalization from 1% to 3%. This has been driven by the considerations that the lower mass ranges are the most susceptible to the knowledge of the background normalization (this is because the EHM mass for and the lower mass signal samples are the most alike). Furthermore, the kinematic phase space of the control region never fully emulates the phase space of the signal region and so actual data analyses are likely to use several control regions to ensure good control of the background normalization, which is likely to increase the systematic uncertainty.

4 Prospects for LHC Run 2 and HL-LHC

Once the full set of selections is applied to the signal and background samples, and the systematic uncertainties are defined, we compute the expected limits on the resonant di-Higgs production cross section multiplied by the branching fraction ().

Figure 5: The dashed line represent the cross section times branching ratio for each mass sample in the set BM. The continuous black line instead represents the predicted 95% C.L. upper limit on with 300 fb of data collected. On the left the limits are obtained performing a cut-and-count experiment on the final number of signal and background events. The limit in this case is entirely driven by the high yield of the background events and the discontinuities are due to the non-continuous assigned systematic uncertainty of the scale factors for the background normalization. The confidence intervals for the expected limit are given at 68% and 95% coverage probability. On the right limits are computed by fitting the HME distribution.

Figure 6: The dashed line represent the cross section times branching ratio for each mass sample in the set BM. The continuous black line instead represents the excluded at 95% C.L. in case 3000 fb of data collected. The confidence intervals for the expected limit are given at 68% and 95% coverage probability.

For calculations of expected limits, shown in Fig. 5, we adopt the modified frequentist criterion CLs [73]. The chosen test statistic, used to determine how signal- or background-like the data are, is based on the profile likelihood ratio [74]. Systematic uncertainties are incorporated in the analysis via nuisance parameters and are treated according to the frequentist paradigm. Results presented in this paper are obtained using asymptotic formulae [75]. The dashed line represents the cross section times branching fraction expected for each mass point within BM (note that benchmark models are chosen as models yielding highest cross-section for each mass range individually, which affects the smoothness of the theoretical prediction curve). The continuous black line represents the predicted 95% C.L. upper limit on the cross-section times branching fraction . In Fig. 5 (left) the limits are shown using a cut-and-count analysis, i.e. applying the whole selection and counting the number of signal and background events expected at the end of the analysis. In this scenario the Heavy Mass Estimator developed in this study has not been used. This simplified method has substantially higher background contamination, with the limit being entirely driven by the background level. The discontinuity in the limit is an artefact of using a non-continuous systematic uncertainties in the scale factors for background normalization (remember that the limits are calculated only for the discrete set of mass points connected by a line to guide reader’s eye). In Fig. 5 (right) the upper limits are computed by fitting the HME distribution, in which case the background under the signal peak is substantially lower and the use of the fit procedure reduces dependence on the background normalization scale factors. Limits are computed assuming a search with 300 fb of integrated luminosity, using both electrons and muons in the final state, and presuming an eventual combination of the ATLAS and CMS data. Weakening of the limit at 350 GeV is due to the similarity of the HME shapes for the background and signal in that mass range. Above 600 GeV the limit trend changes. In this region, the number of events in the signal region becomes almost constant, as can be seen in the HME distribution in Fig. 2 right, and the limit follows the same trend.

In Fig. 6 the limits are shown assuming 3000 fb of integrated luminosity. We obtained such limits by scaling the number of event in signal and background. We did not re-simulate the pileup scenario. The local expected significance of the analysis is computed generating toy-models following the background hypothesis and the same profile likelihood ratio-based CLs technique that has been used for deriving the limits. The local p-value is then converted into significance presented in Fig. 7. Both the limits and the sensitivity correspond to a combination of , and channels with signal and background selection efficiencies equal to those for the channel studied in this paper, and an eventual combination of the CMS and ATLAS results. The sensitivity is shown assuming 300 fb (blue curve) and 3000 fb (red curve) of integrated luminosity per experiment. We find that 3000 fb of data could potentially be sufficient for a potential discovery for GeV, possibly with the exception of the region around = 350 GeV where fluctuations and the look-elsewhere effect may bring the global significance below the conventional 5 threshold. The confidence interval on the sensitivity central value is given at 68% coverage probability.

Figure 7: The colored central lines represent the sensitivity of the analysis assuming 300 fb of integrated luminosity (blue curve) and 3000 fb (red curve). The confidence interval on the sensitivity central value is given at 68% coverage probability.

Finally, we compare the sensitivity of the channel with that of the and channels. We extrapolate the current, public 13 TeV limits from the CMS search for resonant di-Higgs in  [76] and  [77] to 300 fb of integrated luminosity, assuming a naïve (root-squared luminosity) improvement of the present CMS 95% C.L. limit in both final states, and compare with the limits from Fig. 5 (Right)5. The results are shown in Fig. 8, and indicate that while provides the best limits for low masses while may yield better limits than either or in the high mass region. We nevertheless stress that this comparison of , and sensitivities is to be regarded as only indicative, since it is expected that future sensitivity in the and final states improves better than , but a precise estimate of the comparative sensitivity is out of the scope of current work. Still, the comparison suggests that is indeed a competitive search channel for resonant di-Higgs production at the LHC, particularly for high masses.

Figure 8: 13 TeV LHC projected 95% C.L. limits (solid-black lines) on (in pb) for an integrated luminosity fb and assuming an ATLAS-CMS combination, in the final state (as shown in Fig. 5, Right) and in the and final states (through a naïve extrapolation of the resonant di-Higgs 13 TeV CMS analysis in the  [76] and  [77] final states). In all cases the dark (pale) colored bands correspond to the confidence intervals for the expected limit at 68% (95%) coverage probability.

5 Outlook

Exploring the thermal history of EWSB is an important endeavor for particle physics and one for which high energy collisions at the LHC and future colliders can provide invaluable input. Monte Carlo studies imply that for the mass of the observed Higgs boson, EWSB in the SM occurs through a cross over transition. However, the simplest extension of the SM scalar sector – the xSM – may lead to a decidedly different thermal history. In particular, for suitable choices of model parameters, the xSM can generate a strong first order EWPT, thereby fulfilling one of the key conditions for baryogenesis at the electroweak scale. Among the possible signatures of this possibility is resonant di-Higgs production in LHC collisions, catalyzed by the interaction of the singlet-like scalar with pairs of the SM-like Higgs bosons.

In order to fully probe this possibility, it is important to consider a variety of possible final states associated with the di-Higgs decay products. Here, we have considered the channel, with the -bosons decaying leptonically. The presence of two neutrinos in the final state makes the reconstruction of the decaying Higgs-like boson (and thus, of the parent singlet-like scalar) challenging. To address this challenge, we have developed a new Heavy Mass Estimator technique that allows one to achieve the needed mass reconstruction of the singlet-like scalar. Employing the HME and a MVA analysis of signal and background, we show that one is able to exclude the first order EWPT parameter space associated with maximum resonant di-Higgs production cross section with 300 fb of integrated luminosity for GeV and a statistically significant observation over roughly the same mass range with 3000 fb. The projected sensitivity in the channel exceeds in the high mass region ( GeV) that expected from and channels based on a naïve extrapolation of the present CMS 13 TeV public results for the latter channels, indicating that is a competitive di-Higgs LHC search channel for high invariant masses.

Finally, putting our results in the context of other, prospective future collider probes, we note that part, but not all, of the EWPT-viable xSM parameter space accessible in the channel at the LHC would be accessible with precision Higgs studies at the International Linear Collider with TeV and 1 ab integrated luminosity. Full access would require a circular collider ( GeV and 1 ab[78]. Should the HL-LHC exclude this portion of parameter space, then a comprehensive probe would likely require a future 100 TeV collider.


We thank Andreas Papaefstathiou for support with Herwig++. A.S. and T.H. are supported, in part, by the U.S. Department of Energy grant DE-SC0010813, A.S. and L.P. are supported in part by the Qatar National Research Foundation grant NPRP9-328-1-066. M.J.R.M. and P.W. are supported, in part, under U.S. Department of Energy contract DE-SC0011095. J.M.N. is supported in part by the People Programme (Marie Curie Actions) of the European Union Seventh Framework Programme (FP7/2007-2013), REA grant agreement PIEF-GA-2013-625809 and the European Research Council under the European Unions Horizon 2020 program (ERC Grant Agreement no.648680 DARKHORIZONS). J.M.N., M.J.R.M., M.S., and P.W. thank the Munich Institute for Astro- and Particle-Physics (MIAPP), where a portion of this work has been completed. MS is supported in part by the European Commission through the ”HiggsTools” Initial Training Network PITN-GA-2012-316704. We finally thank the CMS collaboration for making Delphes available.


  1. The value of the EWSB vev at the critical temperature, , is inherently gauge-dependent as it is not an observable. Furthermore the standard method for extracting also introduces a separate and spurious gauge-dependence. The consequence is that the conventional criterion for avoiding baryon washout during a first order EWPT (which defines a “strong” first order EWPT), , inherits both sources.
  2. Other potential (and largely subdominant) backgrounds such as Drell-Yan, diboson and single-top can be disregarded, as briefly discussed in Section 3.3.
  3. We recall that we have restricted our analysis to muons, yet it will apply equally to .
  4. The use of discrete mass values in this work is a simplification; in the actual data analysis, training of the MVA would be performed to optimize sensitivity within ranges of target masses . Effectively, this would split the analysis in several sub-analyses, each optimized for a specific range of target masses .
  5. We note the results from Fig. 5 (Right) assume an eventual combination of CMS and ATLAS. This means that a sensitivity improvement should be added to the CMS and limits for a fair comparison. The comparison also assumes SM branching fractions for , which is indeed the case for the xSM.


  1. ATLAS collaboration, G. Aad et al., Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC, Phys. Lett. B716 (2012) 1–29, [1207.7214].
  2. CMS collaboration, S. Chatrchyan et al., Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC, Phys. Lett. B716 (2012) 30–61, [1207.7235].
  3. D. E. Morrissey and M. J. Ramsey-Musolf, Electroweak baryogenesis, New J. Phys. 14 (2012) 125003, [1206.2942].
  4. T. Konstandin, Quantum Transport and Electroweak Baryogenesis, Phys. Usp. 56 (2013) 747–771, [1302.6713].
  5. Planck collaboration, P. A. R. Ade et al., Planck 2013 results. XVI. Cosmological parameters, Astron. Astrophys. 571 (2014) A16, [1303.5076].
  6. A. Sakharov, Violation of CP Invariance, c Asymmetry, and Baryon Asymmetry of the Universe, Pisma Zh.Eksp.Teor.Fiz. 5 (1967) 32–35.
  7. K. Kajantie, M. Laine, K. Rummukainen and M. E. Shaposhnikov, Is there a hot electroweak phase transition at m(H) larger or equal to m(W)?, Phys. Rev. Lett. 77 (1996) 2887–2890, [hep-ph/9605288].
  8. M. Gurtler, E.-M. Ilgenfritz and A. Schiller, Where the electroweak phase transition ends, Phys. Rev. D56 (1997) 3888–3895, [hep-lat/9704013].
  9. M. Laine and K. Rummukainen, What’s new with the electroweak phase transition?, Nucl. Phys. Proc. Suppl. 73 (1999) 180–185, [hep-lat/9809045].
  10. F. Csikor, Z. Fodor and J. Heitger, Endpoint of the hot electroweak phase transition, Phys. Rev. Lett. 82 (1999) 21–24, [hep-ph/9809291].
  11. Y. Aoki, F. Csikor, Z. Fodor and A. Ukawa, The Endpoint of the first order phase transition of the SU(2) gauge Higgs model on a four-dimensional isotropic lattice, Phys. Rev. D60 (1999) 013001, [hep-lat/9901021].
  12. S. Profumo, M. J. Ramsey-Musolf and G. Shaughnessy, Singlet Higgs phenomenology and the electroweak phase transition, JHEP 08 (2007) 010, [0705.2425].
  13. V. Barger, P. Langacker, M. McCaskey, M. J. Ramsey-Musolf and G. Shaughnessy, LHC Phenomenology of an Extended Standard Model with a Real Scalar Singlet, Phys. Rev. D77 (2008) 035005, [0706.4311].
  14. J. R. Espinosa, T. Konstandin and F. Riva, Strong Electroweak Phase Transitions in the Standard Model with a Singlet, Nucl. Phys. B854 (2012) 592–630, [1107.5441].
  15. D. Curtin, P. Meade and C.-T. Yu, Testing Electroweak Baryogenesis with Future Colliders, JHEP 11 (2014) 127, [1409.0005].
  16. S. Profumo, M. J. Ramsey-Musolf, C. L. Wainwright and P. Winslow, Singlet-catalyzed electroweak phase transitions and precision Higgs boson studies, Phys. Rev. D91 (2015) 035018, [1407.5342].
  17. S. J. Huber, T. Konstandin, T. Prokopec and M. G. Schmidt, Electroweak Phase Transition and Baryogenesis in the nMSSM, Nucl. Phys. B757 (2006) 172–196, [hep-ph/0606298].
  18. J. Kozaczuk, S. Profumo, L. S. Haskins and C. L. Wainwright, Cosmological Phase Transitions and their Properties in the NMSSM, JHEP 01 (2015) 144, [1407.4134].
  19. A. Noble and M. Perelstein, Higgs self-coupling as a probe of electroweak phase transition, Phys. Rev. D78 (2008) 063518, [0711.3018].
  20. M. J. Dolan, C. Englert and M. Spannowsky, New Physics in LHC Higgs boson pair production, Phys. Rev. D87 (2013) 055002, [1210.8166].
  21. J. M. No and M. Ramsey-Musolf, Probing the Higgs Portal at the LHC Through Resonant di-Higgs Production, Phys. Rev. D89 (2014) 095031, [1310.6035].
  22. U. Baur, T. Plehn and D. L. Rainwater, Determining the Higgs boson selfcoupling at hadron colliders, Phys. Rev. D67 (2003) 033003, [hep-ph/0211224].
  23. U. Baur, T. Plehn and D. L. Rainwater, Probing the Higgs selfcoupling at hadron colliders using rare decays, Phys. Rev. D69 (2004) 053004, [hep-ph/0310056].
  24. A. Azatov, R. Contino, G. Panico and M. Son, Effective field theory analysis of double Higgs boson production via gluon fusion, Phys. Rev. D92 (2015) 035001, [1502.00539].
  25. F. Kling, T. Plehn and P. Schichtel, Mad-Maximized Higgs Pair Analyses, Submitted to: Phys. Rev. D (2016) , [1607.07441].
  26. M. J. Dolan, C. Englert and M. Spannowsky, Higgs self-coupling measurements at the LHC, JHEP 10 (2012) 112, [1206.5001].
  27. A. Papaefstathiou, L. L. Yang and J. Zurita, Higgs boson pair production at the LHC in the channel, Phys. Rev. D87 (2013) 011301, [1209.1489].
  28. A. J. Barr, M. J. Dolan, C. Englert and M. Spannowsky, Di-Higgs final states augMT2ed – selecting events at the high luminosity LHC, Phys. Lett. B728 (2014) 308–313, [1309.6318].
  29. D. E. Ferreira de Lima, A. Papaefstathiou and M. Spannowsky, Standard model Higgs boson pair production in the ( )( ) final state, JHEP 08 (2014) 030, [1404.7139].
  30. D. Wardrope, E. Jansen, N. Konstantinidis, B. Cooper, R. Falla and N. Norjoharuddeen, Non-resonant Higgs-pair production in the final state at the LHC, Eur. Phys. J. C75 (2015) 219, [1410.2794].
  31. J. K. Behr, D. Bortoletto, J. A. Frost, N. P. Hartland, C. Issever and J. Rojo, Boosting Higgs pair production in the final state with multivariate techniques, Eur. Phys. J. C76 (2016) 386, [1512.08928].
  32. J. Baglio, A. Djouadi, M. M. Grober, R. an Muhlleitner, J. Quevillon and M. Spira, The measurement of the Higgs self-coupling at the LHC: theoretical status, JHEP 04 (2013) 151, [1212.5581].
  33. D. de Florian and J. Mazzitelli, Higgs Boson Pair Production at Next-to-Next-to-Leading Order in QCD, Phys. Rev. Lett. 111 (2013) 201801, [1309.6594].
  34. R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, P. Torrielli et al., Higgs pair production at the LHC with NLO and parton-shower effects, Phys. Lett. B732 (2014) 142–149, [1401.7340].
  35. S. Borowka, N. Greiner, G. Heinrich, S. Jones, M. Kerner, J. Schlenk et al., Higgs Boson Pair Production in Gluon Fusion at Next-to-Leading Order with Full Top-Quark Mass Dependence, Phys. Rev. Lett. 117 (2016) 012001, [1604.06447].
  36. S. Borowka, N. Greiner, G. Heinrich, S. P. Jones, M. Kerner, J. Schlenk et al., Full top quark mass dependence in Higgs boson pair production at NLO, 1608.04798.
  37. ATLAS collaboration, M. Aaboud et al., Search for pair production of Higgs bosons in the final state using proton–proton collisions at TeV with the ATLAS detector, Phys. Rev. D94 (2016) 052002, [1606.04782].
  38. CMS collaboration, V. Khachatryan et al., Search for heavy resonances decaying to two Higgs bosons in final states containing four b quarks, Eur. Phys. J. C76 (2016) 371, [1602.08762].
  39. CMS collaboration, V. Khachatryan et al., Search for two Higgs bosons in final states containing two photons and two bottom quarks in proton-proton collisions at 8 TeV, Phys. Rev. D94 (2016) 052012, [1603.06896].
  40. ATLAS collaboration, G. Aad et al., Searches for Higgs boson pair production in the channels with the ATLAS detector, Phys. Rev. D92 (2015) 092004, [1509.04670].
  41. V. Martin Lozano, J. M. Moreno and C. B. Park, Resonant Higgs boson pair production in the decay channel, JHEP 08 (2015) 004, [1501.03799].
  42. A. Elagin, P. Murat, A. Pranko and A. Safonov, A New Mass Reconstruction Technique for Resonances Decaying to di-tau, Nucl. Instrum. Meth. A654 (2011) 481–489, [1012.4686].
  43. ATLAS collaboration, G. Aad et al., Constraints on new phenomena via Higgs boson couplings and invisible decays with the ATLAS detector, JHEP 11 (2015) 206, [1509.00672].
  44. LHC Higgs Cross Section Working Group collaboration, J. R. Andersen et al., Handbook of LHC Higgs Cross Sections: 3. Higgs Properties, 1307.1347.
  45. ATLAS collaboration, G. Aad et al., Search for an additional, heavy Higgs boson in the decay channel at in collision data with the ATLAS detector, Eur. Phys. J. C76 (2016) 45, [1507.05930].
  46. ATLAS collaboration, G. Aad et al., Combination of searches for , , and resonances in collisions at TeV with the ATLAS detector, Phys. Lett. B755 (2016) 285–305, [1512.05099].
  47. CMS collaboration, S. Chatrchyan et al., Search for a standard-model-like Higgs boson with a mass in the range 145 to 1000 GeV at the LHC, Eur. Phys. J. C73 (2013) 2469, [1304.0213].
  48. CMS collaboration, V. Khachatryan et al., Search for a Higgs Boson in the Mass Range from 145 to 1000 GeV Decaying to a Pair of W or Z Bosons, JHEP 10 (2015) 144, [1504.00936].
  49. M. E. Peskin and T. Takeuchi, Estimation of oblique electroweak corrections, Phys. Rev. D46 (1992) 381–409.
  50. K. Hagiwara, S. Matsumoto, D. Haidt and C. S. Kim, A Novel approach to confront electroweak data and theory, Z. Phys. C64 (1994) 559–620, [hep-ph/9409380].
  51. Gfitter Group collaboration, M. Baak, J. Cúth, J. Haller, A. Hoecker, R. Kogler, K. Mönig et al., The global electroweak fit at NNLO and prospects for the LHC and ILC, Eur. Phys. J. C74 (2014) 3046, [1407.3792].
  52. M. Quiros, Field theory at finite temperature and phase transitions, Helv. Phys. Acta 67 (1994) 451–583.
  53. H. H. Patel and M. J. Ramsey-Musolf, Baryon Washout, Electroweak Phase Transition, and Perturbation Theory, JHEP 07 (2011) 029, [1101.4665].
  54. M. Pietroni, The Electroweak phase transition in a nonminimal supersymmetric model, Nucl. Phys. B402 (1993) 27–45, [hep-ph/9207227].
  55. C. L. Wainwright, CosmoTransitions: Computing Cosmological Phase Transition Temperatures and Bubble Profiles with Multiple Fields, Comput. Phys. Commun. 183 (2012) 2006–2013, [1109.4189].
  56. M. Spannowsky and C. Wymant, Making the most of missing transverse energy: Mass reconstruction from collimated decays, Phys. Rev. D87 (2013) 074004, [1301.0345].
  57. M. Bahr et al., Herwig++ Physics and Manual, Eur. Phys. J. C58 (2008) 639–707, [0803.0883].
  58. S. Alioli, P. Nason, C. Oleari and E. Re, A general framework for implementing NLO calculations in shower Monte Carlo programs: the POWHEG BOX, JHEP 06 (2010) 043, [1002.2581].
  59. DELPHES 3 collaboration, J. de Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaître, A. Mertens et al., DELPHES 3, A modular framework for fast simulation of a generic collider experiment, JHEP 02 (2014) 057, [1307.6346].
  60. C. Collaboration, CMS, the Compact Muon Solenoid: Technical proposal, .
  61. CMS collaboration, S. Chatrchyan et al., Performance of CMS muon reconstruction in collision events at TeV, JINST 7 (2012) P10002, [1206.4071].
  62. M. Cacciari, G. P. Salam and G. Soyez, The Anti-k(t) jet clustering algorithm, JHEP 04 (2008) 063, [0802.1189].
  63. CMS collaboration, S. Chatrchyan et al., Identification of b-quark jets with the CMS experiment, JINST 8 (2013) P04013, [1211.4462].
  64. CMS Collaboration collaboration, CMS, First measurement of the differential cross section for ttbar production in the dilepton final state at sqrts = 13 TeV, Tech. Rep. CMS-PAS-TOP-15-010, CERN, Geneva, 2015.
  65. CMS collaboration, C. Collaboration, Search for resonant Higgs boson pair production in the final state at , .
  66. A. Hocker et al., TMVA - Toolkit for Multivariate Data Analysis, PoS ACAT (2007) 040, [physics/0703039].
  67. C. G. Lester and D. J. Summers, Measuring masses of semiinvisibly decaying particles pair produced at hadron colliders, Phys. Lett. B463 (1999) 99–103, [hep-ph/9906349].
  68. A. Barr, C. Lester and P. Stephens, m(T2): The Truth behind the glamour, J. Phys. G29 (2003) 2343–2363, [hep-ph/0304226].
  69. CMS collaboration, CMS, A search for pair production of new light bosons decaying into muons, Phys. Lett. B752 (2016) 146, [hep-ph/1506.00424].
  70. CMS collaboration, CMS, Search for a light charged Higgs boson decaying to cs in pp collisions at sqrt(s) = 8 TeV, J. High Energy Phyis. 12 (2015) 1, [hep-ph/1510.04252].
  71. CMS collaboration, CMS, Search for two Higgs bosons in final states containing two photons and two bottom quarks, Phys. Rev. D94 (2016) 052012, [hep-ph/1603.06896].
  72. CMS collaboration, CMS, Search for a low-mass pseudoscalar Higgs boson produced in association with a bb pair in pp collisions asqrt(s) = 8 TeV, Phys. Lett. B758 (2016) 296, [hep-ph/1511.03610].
  73. G. Cowan, K. Cranmer, E. Gross and O. Vitells, Asymptotic formulae for likelihood-based tests of new physics, The European Physical Journal C 71 (2011) 1554.
  74. A. L. Read, Presentation of search results: the cls technique, Journal of Physics G: Nuclear and Particle Physics 28 (2002) 2693.
  75. G. Cowan, K. Cranmer, E. Gross and O. Vitells, Asymptotic formulae for likelihood-based tests of new physics, The European Physical Journal C 71 (2011) 1554.
  76. CMS collaboration, C. Collaboration, Search for resonant Higgs boson pair production in the final state using 2016 data, CMS Public Analysis HIG-16-029 (2016) .