Gravitational waves and electroweak baryogenesis in a global study of the extended scalar singlet model
Abstract
We perform a global fit of the extended scalar singlet model with a fermionic dark matter (DM) candidate. Using the most uptodate results from the Planck measured DM relic density, direct detection limits from the XENON1T (2018) experiment, electroweak precision observables and Higgs searches at colliders, we constrain the 7dimensional model parameter space. We also find regions in the model parameter space where a successful electroweak baryogenesis (EWBG) is viable. This allows us to compute the gravitational wave (GW) signals arising from the phase transition, and discuss the potential discovery prospects of the model at current and future GW experiments. Our global fit places a strong upper and lower limit on the second scalar mass, the fermion DM mass and the scalarfermion DM coupling. In agreement with previous studies, we find that our model can simultaneously yield a strong firstorder phase transition and saturate the observed DM abundance. More importantly, the GW spectra of viable points can often be within reach of future GW experiments such as LISA, DECIGO and BBO.
Keywords:
scalar singlet, fermion dark matter, direct detection, electroweak baryogenesis, gravitational wave signals.ADP1827T1075, KCLPHTH/201854 \arxivnumber
1 Introduction
The discovery of the Higgs boson at the LHC Aad:2012tfa (); Chatrchyan:2012xdj () has finally completed the Standard Model (SM) of particle physics. Not only does it provide a new way to study the properties of the Higgs boson, it also offers a way to investigate the details of electroweak symmetry breaking (EWSB). Meanwhile, a more recent observation of the first gravitational wave (GW) signal Abbott:2016blz () and subsequent discoveries Abbott:2016nmj (); Abbott:2017ylp (); Abbott:2017vtc (); Abbott:2017oio (); TheLIGOScientific:2017qsa (); Monitor:2017mdv (); Abbott:2017gyy () have opened up a new window to probe the early history of our universe. In particular, rather violent events such as the firstorder electroweak phase transition (EWPT) would necessarily leave GW imprints. With the current and future generations of ground or spacebased GW experiments, we can hope to observe such signals Caprini:2015zlo (); Weir:2017wfa (); Caprini:2018mtu (). The existence of dark matter (DM) also offers a way to probe the early history of our universe. With the current generation of direct DM searches, experiments are probing the DMnucleon interaction with increasing sensitivity and placing strong limits on the allowed particle DM models.
Motivated by the above experimental probes that are constantly developing, we revisit an extended scalar singlet extension of the SM in this paper. In particular, we focus on two main features of this model. Firstly, it helps to facilitate electroweak baryogenesis (EWBG) Kuzmin:1985mm (); Cohen:1993nk (); Riotto:1999yt (); Morrissey:2012db (), a mechanism that aims to explain the observed matterantimatter asymmetry via a strong firstorder EWPT. In the SM, this phase transition is not firstorder Arnold:1992rz (); Kajantie:1996qd () and thus requires a modification. With an extra scalar, a potential barrier can be generated between the symmetric hightemperature minimum and the EWSB one as the universe cools down Curtin:2014jma (); Kotwal:2016tex (). This leads to a strong firstorder EWPT which can be probed using GWs and standard collider searches Vaskonen:2016yiu (); Ashoorioon:2009nf (); Enqvist:2014zqa (); Artymowski:2016tme (); Tenkanen:2016idg (); Huang:2016cjm (); Kakizaki:2015wua (); Hashino:2016rvx (); Hashino:2016xoj (); Kobakhidze:2016mch (); Chala:2016ykx (); Choi:1993cv (); Beniwal:2017eik (); Kobakhidze:2017mru (); Cai:2017tmh (). Secondly, the new scalar mixes with the SM Higgs boson and provides a portal for a fermion DM to saturate the observed DM abundance Silveira1985136 (); PhysRevD.50.3637 (); Burgess:2000yq ().
Simple DM models with a Higgs portal type interaction are still viable and enjoy a rich interest in the particle physics community Silveira1985136 (); PhysRevD.50.3637 (); Burgess:2000yq (); Espinosa:2008kw (); Alanne:2014bra (); Falkowski:2015iwa (); Buttazzo:2015bka (); Heikinheimo:2016yds (); Balazs:2016tbi (); Lewis:2017dme (); Ghorbani:2017jls (); Chen:2017qcz (); Kamon:2017yfx (); Ettefaghi:2017vbh (); Baker:2017zwx (); Bernal:2018ins (). In our study, we focus on a singlet fermion DM model which was first introduced in Ref. Kim:2008pp () and subsequently improved in Ref. Baek:2011aa (). After the discovery of a SMlike Higgs boson at the LHC, the model was revisited in Ref. Baek:2012uj () in the context of vacuum stability (see also Ref. Espinosa:2011ax ()). Here it was pointed out that the model is stable and perturbative up to the Planck scale for a GeV Higgs boson. In light of EWBG, the model was first studied in Ref. Fairbairn:2013uta () and more recently in Ref. Li:2014wia (). Using a Monte Carlo scan of the model parameter space, the model was shown to realise a strong firstorder phase transition without conflicting with any bounds from direct DM searches, electroweak precision observables (EWPO) and latest Higgs data from the LHC.
We aim to perform the most comprehensive and uptodate study of the extended scalar singlet model with a fermionic DM candidate. In our global fit, we include the latest results from the Planck measured DM relic density Ade:2015xua (), direct detection limits from the XENON1T (2018) experiment Aprile:2018dbl (), EWPO Haller:2018nnx () and Higgs searches at colliders Bechtle:2013wla (); Bechtle:2013xfa (). We also find regions in the model parameter space where a successful EWBG is viable, compute the resulting GW spectra, and check the discovery prospects of the model at current and future GW experiments. In agreement with previous studies, we confirm that our model with additional couplings to the SM Higgs boson can simultaneously explain the observed DM abundance and matterantimatter asymmetry; this was not possible in the symmetric case studied in our previous work Beniwal:2017eik (). We also find that our global fit places a strong upper and lower limit on the second scalar mass , fermion DM mass and the scalarfermion DM coupling . In addition, the GW spectra of viable points can often be within reach of future GW experiments such LISA, DECIGO and BBO.
The rest of the paper is organised as follows. In section 2, we introduce the extended scalar singlet model with a fermionic DM candidate. After taking note of the free parameters of our model, we describe a set of constraints and likelihoods used in our global fit in section 3. Our model results and conclusions are presented in sections 4 and 5 respectively. Appendices A, B, C and D provide supplementary details for understanding various expressions in the paper.
2 Singlet fermion dark matter model
We extend the SM by adding a new real scalar singlet and a Dirac fermion DM field . The fermion DM is assumed to be living in the hidden sector and communicates with the SM particles only via the new scalar . The model Lagrangian is given by Baek:2011aa ()
(1) 
where is the SM Lagrangian,
(2)  
(3)  
(4) 
In general, a linear term in the field is allowed by symmetry. However, such a term can be removed by performing a constant shift in which also redefines , , , and .^{4}^{4}4The parameter appears in the SM Higgs potential, see Eq. (6). In writing the above Lagrangians, we have assumed that these parameters are defined after a constant shift in . If we set , we can see that the above Lagrangian becomes symmetric under , i.e., it is even in . In this case, the fermion DM is decoupled and becomes a hidden DM candidate, whereas the scalar serves as a new DM candidate and reproduces the scalar Higgs portal model Cline:2013gha (); Beniwal:2015sdl (); He:2016mls (); Escudero:2016gzx (); Wu:2016mbe (); Banerjee:2016vrp (); Casas:2017jjg (); Beniwal:2017eik (); Athron:2017kgt (); Hoferichter:2017olk (); Athron:2018ipf ().
With an extra scalar field, the treelevel scalar potential is given by
(5) 
where and can be read directly from Eqs. (2) and (4) respectively. The SM part of the potential reads
(6) 
where
(7) 
is the SM Higgs doublet and are the Goldstone bosons.
In general, both and can develop nontrivial vacuum expectation values (VEVs). At , these are denoted by and respectively, i.e.,
(8) 
After EWSB, we can expand and in the unitary gauge as
(9) 
where fields represent quantum fluctuations around the VEVs. Using the results presented in Appendix A, we arrive at the following EWSB conditions
(10)  
(11) 
The portal interaction Lagrangian in Eq. (4) induces a mixing between the and fields. Thus, the squared mass matrix
(12) 
is nondiagonal. As shown in Appendix A, its elements are given by
(13) 
The squared mass matrix in Eq. (12) can be diagonalised by rotating the interaction eigenstates into the physical mass eigenstates as
(14) 
where is the mixing angle. Thus, for small mixing, is a SMlike Higgs boson, whereas is dominated by the scalar singlet.
3 Constraints and likelihoods
In light of the recent discovery of a SMlike Higgs boson at the LHC Aad:2012tfa (); Chatrchyan:2012xdj (), we set
(18) 
Thus, the model is completely described by the following 7 free parameters
(19) 
The remaining parameters in Eqs. (2), (4) and (6) can be expressed as (see Appendix B)
(20)  
(21)  
(22)  
(23)  
(24) 
To study the phenomenology of our model, we implement the extended scalar singlet and fermion DM model in the LanHEP_v3.2.0 Semenov:2014rea () package. For the calculation of the fermion DM relic density and Higgs decay rates, we use micrOMEGAs_v4.3.5 Belanger:2014vza () which relies on the CalcHEP Belyaev:2012qa () package.
We make parameter inferences by adopting a frequentist approach and performing 7dimensional scans of the model parameter space using the Diver_v1.0.4 Workgroup:2017htr () package.^{5}^{5}5http://diver.hepforge.org The combined loglikelihood used in our global fit is
(25) 
where

: loglikelihood for the Planck measured DM relic density, see subsection 3.1;

: loglikelihood for the direct detection limits from the XENON1T (2018) experiment, see subsection 3.2;

: loglikelihood for the EWBG constraint, see subsection 3.3;

: loglikelihood for the electroweak precision observables (EWPO) constraint, see subsection 3.4;

: loglikelihood for the direct Higgs searches performed at the LEP, Tevatron and the LHC, see subsection 3.5;

: loglikelihood for the Higgs signal strength and mass measurements performed at the LHC, see subsection 3.5.
Here denotes the free parameters of our model. These are uniformly sampled over their ranges shown in Table 1 in either flat or logarithmic space.
In the following subsections, we outline the details of all constraints and likelihoods used in our global fit.
Parameter  Minimum  Maximum  Prior type 
log  
flat  
flat  
log  
flat  
log  
log 
3.1 Thermal relic density
From the Planck satellite’s observation of the temperature and polarization anisotropies in the cosmic microwave background (CMB), a strong bound on the presentday abundance of the DM particles can be extracted. The latest results indicate Ade:2015xua ()
(26) 
where is the density parameter, is the critical mass density and is the reduced Hubble constant.
In our model, the Dirac fermion is the DM candidate. Its relic density is mainly determined by an channel annihilation into SM particles via an exchange. Annihilation into , and final states are also possible via the  and channels. Due to a mixing between the interaction eigenstates , the decay rates go as
(27)  
(28) 
where is a general SM final state, e.g., quarks, leptons or gauge bosons. Depending on the mixing angle , three scenarios are possible.

: In this case, is a SMlike Higgs boson, whereas is a scalar singlet. Thus, the only allowed final states from the fermion DM annihilation are , and .

: This corresponds to a maximal mixing between the interaction eigenstates. All of the allowed final states from the fermion DM annihilation are shown in Fig. 1.

: In this case, is a scalar singlet, whereas is a SMlike Higgs boson. Similar to the case, the only allowed final states from the fermion DM annihilation are , and .
With two scalar mediators and , the annihilation rate of the fermion DM into SM particles is enhanced when . At these two resonances, the fermion DM relic density drops rapidly with increasing scalarfermion DM coupling . For the fermion DM to account for the observed DM abundance, i.e., , smaller values of are required to compensate for the enhanced DM annihilation rate into SM particles.
In order to address the strong possibility of a multicomponent dark sector, we define a relic abundance parameter Cline:2012hg (); Cline:2013gha (); Beniwal:2015sdl () as
(29) 
where is the Planck measured central value in Eq. (26). Consequently, the indirect and direct detection rates must be scaled by and respectively.^{6}^{6}6In our study, we do not include any indirect detection limits as the fermion DM annihilation rate into SM particles is wave suppressed Agrawal:2010fh (). However, when a pure pseudoscalar, parityviolating interaction term is introduced, the resulting indirect detection limits can be sizeable Esch:2013rta (); Bagherian:2014iia (); Franarin:2014yua (); Kim:2016csm (); Ghorbani:2014qpa (); Balazs:2015boa (); Athron:2018hpc (); Kim:2018uov (). In regions of the model parameter space where , parameter points are robustly excluded by the relic density constraint.
We investigate both possibilities of our model to either account for all or part of the observed DM abundance. In the former case, we use a Gaussian likelihood function with a central value equal to the Planck measured one and a combined uncertainty equal to the Planck measured uncertainty with a 5% theoretical error.^{7}^{7}7A possible source of theoretical uncertainty is in our relic density calculations as performed in micrOMEGAs. In the latter case, we instead use a Gaussian likelihood function as an upper limit and require the parameter points to satisfy . The results for both of these scenarios will be discussed in more detail in section 4.
3.2 Direct detection
Direct detection experiments aim to measure the recoil of a nucleus from an elastic scattering off a DM particle. Such an event generates a typical recoil energy on the order of a few keV. As most radioactive elements and highenergy cosmic rays induce nuclear recoils with energies well above this value, direct DM searches must be conducted in deep underground laboratories to shield them from potential background sources.
In our model, the DMquark interaction proceeds via a channel exchange of particles. With two neutral scalar mediators , the resulting DMnucleus interaction is nuclear spinindependent (SI). The SI DMnucleus crosssection is given by
(30) 
where is the DMnucleus reduced mass and are the number of protons (neutrons) in the target nucleus . The dimensionful parameters are the effective DMnucleon couplings. These are given by (see Appendix C)
(31) 
where ,
(32) 
is the Higgsnucleon coupling and
(33) 
are the hadronic matrix elements.
For isospin conserving couplings , the DMnucleus crosssection in Eq. (30) is enhanced by a factor of . This is expected as the matrix element for a SI interaction involves a coherent sum over the individual protons and neutrons in the target nucleus . For this reason, direct detection experiments rely on heavy target materials with large to better constrain the DMnucleon crosssection . In our model, it is given by
(34) 
where is the DMnucleon reduced mass, MeV and Cline:2013gha ().
Currently, the best upper limits on the SI DMnucleon crosssection comes from the XENON1T (2018) experiment Aprile:2018dbl (). To constrain the model parameter space from the XENON1T experiment, we use a onesided Gaussian likelihood function, i.e., we require the parameter points to satisfy^{8}^{8}8The official XENON1T (2018) limits are only available for DM masses up to TeV. Beyond TeV, we perform a linear extrapolation of the limit due to the reduced DM number density.
(35) 
where is the 90 C.L. upper limit from the XENON1T experiment and
(36) 
is the effective SI DMnucleon crosssection. The scaling of by is done to suppress signals when . In regions of the model parameter space where , parameter points are already ruled out by the relic density constraint.
We also include a theoretical uncertainty of 5 in our analysis. This can easily arise from the uncertainties associated with the nuclear physics, DM halo and velocity distribution parameters. For a recent review, see Ref. Workgroup:2017lvb ().
3.3 Electroweak baryogenesis (EWBG)
In our model, the VEV of the new scalar does not initially have to be zero. Thus, the transition pattern can be . At low temperatures, the latter minimum evolves slowly to become the electroweak minimum at , i.e., . The initial transition can break the electroweak symmetry by tunnelling through a potential barrier to the broken phase minimum. This transition can proceed via nucleation of bubbles of the broken phase which results in a departure from thermal equilibrium Kuzmin:1985mm (); Cohen:1993nk (); Riotto:1999yt (); Morrissey:2012db (). In addition, it can generate a significant gravitational wave (GW) signal Grojean:2006bp ().
Using the standard notation, we define a strong firstorder phase transition by
(37) 
where is the Higgs VEV at temperature . However, one has to keep in mind that the calculation of the baryon asymmetry remaining after the transition is quite complicated. This leads to a slightly different exact lower bound on Quiros:1999jp (); Funakubo:2009eg (); Curtin:2014jma (); Katz:2014bha (); Fuyuto:2014yia ().
To find regions in the model parameter space where a successful EWBG is viable, we first find the minima of the effective potential (see Appendix D) numerically, and compute the critical temperature at which the initial and symmetry breaking minima are degenerate. This allows us to compute the dimensionless parameter (the Higgs VEV at the critical temperature ) and constrain parts of the 7dimensional model parameter space, i.e., parameter points are excluded if they lead to a too weak phase transition. Specifically, we use a onesided Gaussian likelihood function and require the parameter points to satisfy
(38) 
as a conservative limit. A theoretical uncertainity of 5 on the resulting values is assumed to obtain a smooth likelihood function. The actual uncertainty can be much larger as the value of required to facilitate EWBG is not yet settled.
In addition, parameter points are also excluded if they exhibit any of the following three features.

Incorrect minimum at : This situation arises when the electroweak vacuum is not the true minimum of the potential at .

Runaway directions in the potential: This occur when the and field values in the symmetric or broken phase are too large, or if the potential in Eq. (5) is unbounded from below in the general and directions, i.e., when .

Nonperturbative couplings: This situation arises when . In this case, our 1loop treatment of the effective potential is not reliable.
We also perform a complete analysis of the phase transition in this model by following our previous work on the symmetric case, i.e., scalar Higgs portal Beniwal:2017eik () and a very recent update on the calculation of the phase transition dynamics Ellis:2018mja (). In particular, we find the percolation temperature at which the phase transition truly completes. This is used to compute the GW signals arising from the phase transition, and discuss the potential discovery prospects of the model at current and future GW experiments. For more details, see section 4.
3.4 Electroweak precision observables (EWPO)
With an extra scalar, our model can induce corrections to the gauge boson selfenergy diagrams. Its effect on the electroweak precision observables (EWPO) can be parametrised by the oblique parameters , and Peskin:1991sw (). The and selfenergies ( and respectively) are not modified as the new scalar is electrically neutral. Thus, only the and boson selfenergies are subject to corrections.
In our model, the oblique parameters are shifted from their SM values by Baek:2011aa ()
(39)  
(40)  
(41) 
where for , is the boson mass, and . The loop functions and are given by Grimus:2008nb ()
(42)  
(43) 
These are also plotted in Fig. 2. From Eqs. (39)–(41), it is evident that
(44) 
Thus, for large , is required, whereas large mixing angles are compatible with the EWPO constraint provided .
Using the SM reference as GeV and GeV, the most recent global electroweak fit gives Haller:2018nnx ()
(45) 
and the following correlation matrix
(46) 
To constrain the model parameter space from the EWPO, we use the following likelihood function Profumo:2014opa ()
(47) 
where denotes the central values for the shifts in Eq. (45), is the covariance matrix, is the correlation matrix in Eq. (46) and are the associated errors in Eq. (45).
3.5 Higgs searches at colliders
Due to a mixing between the interaction eigenstates , the coupling strengths between the mass eigenstates and SM particles are modified with respect to the SM expectation. The effective squared couplings of to SM particles are Li:2014wia ()
(48) 
where refers to a SM quark, lepton or gauge boson. For the loopinduced processes, these are given by Djouadi:2005gi ()
(49) 
where . With modified branching ratios of into SM particles, the scalar sector of our model can be constrained using the direct Higgs searches performed at the lepton (e.g., LEP) and hadron (e.g., Tevatron, LHC) colliders.
To constrain the model parameter space from the direct Higgs searches performed at the Tevatron and the LHC, we use the HiggsBounds_v4.3.1 Bechtle:2013wla () package. From the model predictions for the two scalar masses, total decay widths, branching ratios, and effective squared couplings defined in Eqs. (48) and (49), HiggsBounds computes and compares the predicted signal rates for the search channels considered in multiple experimental analyses. By comparing the predicted signal rates against the expected and observed crosssection limits from the direct Higgs searches, it determines whether or not a given parameter point is excluded at 95 C.L..
For the two physical scalars , the signal strengths are given by Li:2014wia ()
(50)  
(51) 
In the absence of invisible and cross Higgs decay modes, scales as . However, when these decay modes are kinematically allowed, they suppress the signal strength with respect to the SM expectation. Thus, the scalar sector of our model can also be constrained using the Higgs signal strength and mass measurements performed at the LHC.
Experiment  Channel  Obs. signal strength  Ref. 
ATLAS  Aad:2015gba ()  
ATLAS  Aad:2015gba ()  
ATLAS  Aad:2015gba ()  
ATLAS  Aad:2015gba ()  
ATLAS  Aad:2015gba ()  
CMS  Chatrchyan:2013iaa ()  
CMS  Chatrchyan:2013mxa ()  
CMS  Khachatryan:2014ira ()  
CMS  Chatrchyan:2014vua ()  
CMS  Chatrchyan:2014vua () 
To constrain the model parameter space from the Higgs signal strength and mass measurements, we use the HiggsSignals_v1.4.0 Bechtle:2013xfa () package. Assuming a Gaussian probability density function (p.d.f.) for the two scalar masses, we compute a chisquare using the peakcentered method.^{9}^{9}9A theoretical mass uncertainty of zero is assumed for both scalars as is fixed, whereas is a free model parameter. In this method, is evaluated by assigning, for each signal (or peak) observed in multiple experimental analyses (see Table 2), a combination of the two Higgs bosons from our model provided their masses lie within the experimental resolution of an analysis Stal:2013hwa (). Following the assignment, a is evaluated by comparing the signal strength measurement for the peak to the model predicted signal strength. When a mass measurement is also available (e.g., from channels with a good massresolution such as the decay mode), a corresponding is also evaluated by comparing the model predicted and observed Higgs boson mass. Thus, the total is given by^{10}^{10}10For more details on the functional form of individual chisquares, see Ref. Bechtle:2013xfa ().
(52) 
In situations where more than one Higgs boson can contribute to a signal (as in our case), an optimal assignment of the Higgs bosons to the signals is achieved by minimising the overall . The predicted signal strengths of the two scalars are added incoherently, assuming negligible interference effects. Finally, the computed is used to define a Higgs signal strength loglikelihood as
(53) 
Thus, a large indicates a large deviation between the model predicted signal strength and the bestfit value for a fixed Higgs boson mass, and vice versa.
4 Results
We perform scans of our 7D model parameter space using Diver_v1.0.4 Workgroup:2017htr () with lambdajDE = true, NP = 50,000 and convthresh = 10. To efficiently sample all parts of the parameter space (even the degenerate ones), we also run several targeted scans and combine the output chains to obtain highquality profile likelihood (PL) plots.
We present our model results in the form of 1 and 2dimensional PL plots. For a model parameter where , a 1D PL is defined as
(54) 
Thus, is a function of only, i.e., all other parameters are profiled out. Similarly, a 2D PL is defined as
(55) 
Thus, is a function of and only. Using Eqs. (54) and (55), we can define a PL ratio Cowan:2010js () as
(56) 
where is the bestfit point, i.e., a parameter point that maximises the total likelihood function . Using Wilks’ theorem Wilks:1938dza (), Eq. (56) can be used to construct contours corresponding to C.L. regions.
In the following subsections, we present our model results in the form of 1D and 2D PL plots. These are generated using the pippi_v2.0 Scott:2012qh () package.
4.1 EWBG only
We start by finding regions in the model parameter space where a successful EWBG is viable. This is achieved by performing a 7D scan of the model using only the loglikelihood, i.e.,
(57) 
where is defined in subsection 3.3. The resulting 2D PL plots are shown in Fig. 3. In the dark blue regions where the PL ratio , the dimensionless parameter and a successful EWBG is viable. To understand the results in more detail, we go over each panel in Fig. 3 onebyone.

plane: For TeV, all values of and some combination of 5 profiled out parameters (namely , , , and ) give and maximise the loglikelihood, thus everywhere. Due to the dependence of in Eq. (22), large values of should lead to runaway directions, , and/or nonperturbative coupling, . With , a large contribution from to can be suppressed. However, this choice of makes in Eq. (20) nonperturbative as its contribution appears as . Ultimately, the solution is to choose a small value for as its contribution in Eq. (22) appears as . In addition, small values of can also help in keeping . Thus, for TeV, large values of can facilitate EWBG.
For TeV and GeV, the white region () is disfavoured as it leads to . This is expected as the contribution from in Eq. (22) is dominant at large values. With large , no choice of , and can keep . In fact, the requirement translates into an upper limit on as a function of , , and . Using Eq. (22), we get
(58) For a fixed and , Eq. (58) has 3 degrees of freedom. As , and are profiled over, it is nontrivial to predict the exact shape of the upper limit on as a function of . The upper limit also weakens as increases. The net result is that for TeV, GeV is required to facilitate EWBG.

plane: Similar to the plane for TeV, some combination of the profiled out parameters gives for all values of . However, when TeV and , the Higgs quartic coupling in Eq. (20) becomes nonperturbative. In fact, the requirement translates into the following upper limit on as a function of
(59) When , the above condition is satisfied for all values of . Thus, a successful EWBG is viable at large values of . On the other hand, when , Eq. (59) imposes the strongest upper limit on , namely TeV. As , the upper limit on becomes weaker, as is evident from the plot.

and planes: In these two planes, all possible combinations of , and profiled out parameters give . Thus, the PL ratio is roughly flat and equal to 1 everywhere; hence, we do not show these planes in Fig. 3. In fact, the likelihood is weakly dependent on and as expected from Eq. (22). For instance, at large values of or which would give or , small values of can be chosen to avoid such situations.

plane: For TeV, all values of give . As does not appear directly in Eqs. (20) and (22), the likelihood is weakly dependent on . This is expected as the contribution from to the effective potential appears only at 1loop order.
Figure 3: 2D profile likelihood (PL) plots from a 7D scan of our model using only the electroweak baryogenesis (EWBG) constraint. The contour lines mark out the (68.3) and (95.4) C.L. regions. In regions where , a successful EWBG is viable as (see text for more details). The parameter planes and are not shown as they are unconstrained by the EWBG constraint. For TeV and TeV, no combination of the profiled out parameters can keep . On the other hand, when TeV, the positive contribution from to the effective potential can be cancelled out by the negative contribution from , thereby giving and .

plane: For , all values of and profiled out parameters give , and maximise the likelihood. However, values of lead to runaway directions in the potential as the contribution from in the 1loop corrections become large. This pushes the broken phase minimum too far away from the origin and suppresses the vacuum decay probability.
In summary, it is not difficult to facilitate a successful EWBG in our model. For any specific model parameter, usually some combination of the remaining parameters give viable points even if the parameter in question causes problems. For instance, large values of generally push up the EWSB minimum and cause it to not become the global minimum at . However, this effect can be counteracted by choosing a large value for which gives a large negative contribution to the effective potential. One exception is which always generates runaway directions in the effective potential. For the remaining model parameters, namely , the 1D PL ratio is roughly flat and equal to 1 for all parameter values. Thus, we do not show the 1D PL plots for our model parameters.
4.2 Global fit
With some intuition on the choice of free model parameters that can facilitate a successful EWBG, we present results from a global fit of our model using the total loglikelihood function in Eq. (3). In practice, we consider two scenarios in which the fermion DM accounts for either a small fraction or all of the observed DM abundance. In the former case, we use a relic density likelihood that is onesided Gaussian, whereas in the latter, we use a Gaussian likelihood. For more details, see subsection 3.1.
4.2.1 Scenario I:
The resulting 2D PL plots from our 7D scans are shown in Fig. 4. For GeV, the parameter planes are ruled out by the observed Higgs signal strength measurements, EWPO and direct Higgs searches performed at the LEP experiment. As the decay channel is kinematically allowed and dominant in this region for all values of the mixing angle , it reduces the SMlike Higgs signal strength with respect to SM expectation, see Eq. (50). This translates into a large in Eq. (53) and is thus disfavoured.
To understand the remaining set of results in more detail, we go over each panel in Fig. 4 onebyone.

plane: For TeV, the parameter planes are ruled out by the EWBG constraint as they either lead to runaway directions, , or nonperturbative couplings, . Although, some combinations of the profiled out parameters can give a successful EWBG at large values of (see Fig. 3), they are often not compatible with the remaining constraints. This is especially true for the EWPO constraint which only depends on and . For large , is required in order to satisfy the EWPO constraint.

plane: We see that the model is allowed by all constraints for a range of low and high values provided . This is expected as the second scalar is decoupled in this regime and gives no new contribution to the observed Higgs signal strengths. However, when , the two scalars are indistinguishable from the point of direct Higgs searches and Higgs signal strength measurements. As is evident in Eqs. (31) and (44), direct detection and EWPO constraints respectively are also relaxed in this regime. The net result is that all values of are allowed when .

and planes: These parameter planes are mostly unconstrained by our global fit except for (excluded by the Higgs signal strength measurements) and TeV (ruled out by the EWBG constraint). For TeV, large values of are required to facilitate a successful EWBG.

plane: For , the fermion DM can only annihilate into light SM quarks, thereby giving . On the other hand, is constrained by the DM relic density and XENON1T limits. When , all values of up to TeV are allowed by the Planck measured relic density and XENON1T limits; this region appears in the plot as a horizontal band. In this band, small values of can yield a fermion DM relic density and DMnucleon crosssection that is compatible with the Planck measured value and XENON1T limit respectively.
For , the region is disfavoured by either the Planck measured relic density or XENON1T limit. This is generally expected from an incompatibility between small values of which are favoured by the XENON1T limit (as it gives a small DMnucleon crosssection but disfavoured by the relic density constraint (as it gives ) and vice versa.
The diagonal band corresponds to the second resonance . Similar to the first resonance , all points in this band are allowed by the relic density and XENON1T limits. As is profiled over, small values of can easily give and . On the other hand, when TeV, parameter points are disfavoured by the EWPO and EWBG constraints. For TeV, the region is robustly excluded by the combined constraints.
Figure 4: 2D PL plots from a global fit of our model assuming . The contour lines mark out the (68.3) and (95.4) C.L. regions. 
plane: In this plane, a lower limit on comes from the DM relic density constraint as smaller values of lead to an overabundance of the fermion DM in the universe today. This lower limit becomes weaker as increases. For TeV, the coupling becomes nonperturbative, thus this region is disfavoured. Similarly, values of are disfavoured by the EWBG constraint as they lead to runaway directions in the potential, see Fig. 3.
In Fig. 5, we show the 1D PL plots for the parameters , and