Gravitational waves and electroweak baryogenesis in a global study of the extended scalar singlet model

We perform a global fit of the extended scalar singlet model with a fermionic dark matter (DM) candidate. Using the most up-to-date results from the $\mathit{Planck}$ measured DM relic density, direct detection limits from the XENON1T (2018) experiment, electroweak precision observables and Higgs searches at colliders, we constrain the 7-dimensional model parameter space. We also find regions in the model parameter space where a successful electroweak baryogenesis (EWBG) can be 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 $\mathit{and}$ lower limit on the second scalar mass, the fermion DM mass and the scalar-fermion DM coupling. In agreement with previous studies, we find that our model can simultaneously yield a strong first-order 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.


Introduction
The discovery of the Higgs boson at the LHC [1,2] 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 [3] and subsequent discoveries [4][5][6][7][8][9][10] 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/space-based GW experiments, we can hope to observe such signals [11][12][13].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 DM-nucleon 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) [14][15][16][17], a mechanism that aims to explain the observed matter-antimatter asymmetry via a strong first-order EWPT.In the SM, this phase transition is not first-order [18,19] and thus requires a modification.With an extra scalar, a potential barrier can be generated between the symmetric high-temperature minimum and the EWSB one as the universe cools down [20,21].This leads to a strong first-order EWPT which can be probed using GWs and standard collider searches [22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40].Secondly, the new scalar mixes with the SM Higgs boson and provides a portal for a fermion DM to saturate the observed DM abundance [41][42][43].
Simple DM models with a Higgs portal type interaction are still viable and enjoy a rich interest in the particle physics community [41][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57].In our study, we focus on a singlet fermion DM model which was first introduced in Ref. [58] and subsequently improved in Ref. [59].After the discovery of a SM-like Higgs boson at the LHC, the model was revisited in Ref. [60] in the context of vacuum stability (see also Ref. [61]).Here it was pointed out that the model is stable and perturbative up to the Planck scale for a 125 GeV Higgs boson.In light of EWBG, the model was first studied in Ref. [62] and more recently in Ref. [63].Using a Monte Carlo scan of the model parameter space, the model was shown to realise a strong first-order 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 up-to-date 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 [64], direct detection limits from the XENON1T (2018) experiment [65], EWPO [66] and Higgs searches at colliders [67,68].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 matter-antimatter asymmetry; this was not possible in the Z 2 symmetric case studied in our previous work [35].We also find that our global fit places a strong upper and lower limit on the second scalar mass m H , fermion DM mass m ψ and the scalar-fermion DM coupling g S .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.

Singlet fermion dark matter model
We extend the SM by adding a new real scalar singlet S 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 S. The model Lagrangian is given by [59] where L SM is the SM Lagrangian, (2.4) In general, a linear term in the S field is allowed by symmetry.However, such a term can be removed by performing a constant shift in S which also redefines µ2 S , µ 2 Φ , µ 3 , g S and µ ΦS . 1 In writing the above Lagrangians, we have assumed that these parameters are defined after a constant shift in S. If we set µ 3 = g S = µ ΦS = 0, we can see that the above Lagrangian becomes Z 2 symmetric under S → −S, i.e., it is even in S. In this case, the fermion DM ψ is decoupled and becomes a hidden DM candidate, whereas the scalar S serves as a new DM candidate and reproduces the scalar Higgs portal model [35,[69][70][71][72][73][74][75][76][77][78].
With an extra scalar field, the tree-level scalar potential is given by where V S and V portal can be read directly from Eqs. (2.2) and (2.4) respectively.The SM part of the potential reads where is the SM Higgs doublet and (G ± , G 0 ) are the Goldstone bosons.
In general, both φ and S can develop non-trivial vacuum expectation values (VEVs).At T = 0, these are denoted by v 0 and s 0 respectively, i.e., 0|φ|0 After EWSB, we can expand Φ and S in the unitary gauge as where (ϕ, s) fields represent quantum fluctuations around the T = 0 VEVs.Using the results presented in Appendix A, we arrive at the following EWSB conditions (2.11) The portal interaction Lagrangian in Eq. (2.4) induces a mixing between the ϕ and s fields.Thus, the squared mass matrix is non-diagonal.As shown in Appendix A, its elements are given by 13) The squared mass matrix in Eq. (2.12) can be diagonalised by rotating the interaction eigenstates (ϕ, s) into the physical mass eigenstates (h, H) as where α is the mixing angle.Thus, for small mixing, h is a SM-like Higgs boson, whereas H is dominated by the scalar singlet.
For the tree-level scalar potential in Eq. (2.5) to be bounded from below, the following conditions must be satisfied (see Appendix A) (2.15) After EWSB, the fermion DM Lagrangian in Eq. (2.3) becomes where is the physical fermion DM mass.

Constraints and likelihoods
In light of the recent discovery of a SM-like Higgs boson at the LHC [1,2], we set Thus, the model is completely described by the following 7 free parameters The remaining parameters in Eqs.(2.2), (2.4) and (2.6) can be expressed as (see Appendix B) ) To study the phenomenology of our model, we implement the extended scalar singlet and fermion DM model in the LanHEP v3.2.0 [79] package.For the calculation of the fermion DM relic density and Higgs decay rates, we use micrOMEGAs v4.3.5 [80] which relies on the CalcHEP [81] 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 [82] package. 2 The combined log-likelihood used in our global fit is where • ln L Ωh 2 (θ): log-likelihood for the Planck measured DM relic density, see subsection 3.1; • ln L XENON1T (θ): log-likelihood for the direct detection limits from the XENON1T (2018) experiment, see subsection 3.2; • ln L vc/Tc (θ): log-likelihood for the EWBG constraint, see subsection 3.3; • ln L EWPO (θ): log-likelihood for the electroweak precision observables (EWPO) constraint, see subsection 3.4; • ln L HB (θ): log-likelihood for the direct Higgs searches performed at the LEP, Tevatron and the LHC, see subsection 3.5; • ln L HS (θ): log-likelihood for the Higgs signal strength and mass measurements performed at the LHC, see subsection 3.5.
Here θ ≡ (m H , s 0 , µ 3 , λ S , α, m ψ , g S ) 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.

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 present-day abundance of the DM particles can be extracted.The latest results indicate [64] Ω DM h 2 = 0.1188 ± 0.0010, ( where Ω DM = ρ DM /ρ c is the density parameter, ρ c = 3H 2 0 M 2 p is the critical mass density and h = H 0 /(100 km s −1 Mpc −1 ) is the reduced Hubble constant.
In our model, the Dirac fermion ψ is the DM candidate.Its relic density is mainly determined by an s-channel annihilation into SM particles via an h/H exchange.Annihilation into hh, HH and hH final states are also possible via the t-and u-channels.Due to a mixing between the interaction eigenstates (ϕ, s), the decay rates go as where X is a general SM final state, e.g., quarks, leptons or gauge bosons.Depending on the mixing angle α, three scenarios are possible.
1. α = 0: In this case, h is a SM-like Higgs boson, whereas H is a scalar singlet.Thus, the only allowed final states from the fermion DM annihilation are hh, HH and hH.
2. α = π/4: 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.
3. α = π/2: In this case, h is a scalar singlet, whereas H is a SM-like Higgs boson.Similar to the α = 0 case, the only allowed final states from the fermion DM annihilation are hh, HH and hH.With two scalar mediators h and H, the annihilation rate of the fermion DM into SM particles is enhanced when m ψ ∼ m h,H /2.At these two resonances, the fermion DM relic density Ω ψ h 2 drops rapidly with increasing scalar-fermion DM coupling g S .For the fermion DM to account for the observed DM abundance, i.e., Ω ψ h 2 = Ω DM h 2 , smaller values of g S 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 [69,70,83] as where Ω DM h 2 = 0.1188 is the Planck measured central value in Eq. (3.9).Consequently, the indirect and direct detection rates must be scaled by f 2 rel and f rel respectively. 3In regions of the model parameter space where f rel > 1, 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. 4In the latter case, we instead use a Gaussian likelihood function as an upper limit and require the parameter points to satisfy f rel ≤ 1.The results for both of these scenarios will be discussed in more detail in section 4.

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 E R on the order of a few keV.As most radioactive elements and high-energy 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 DM-quark interaction proceeds via a t-channel exchange of h/H particles.With two neutral scalar mediators (h, H), the resulting DM-nucleus interaction is nuclear spin-independent (SI).The SI DM-nucleus cross-section is given by where µ ψN = m ψ m N /(m ψ + m N ) is the DM-nucleus reduced mass and Z (A − Z) are the number of protons (neutrons) in the target nucleus N .The dimensionful parameters (G p , G n ) are the effective DM-nucleon couplings.These are given by (see Appendix C) where N ∈ (p, n), is the Higgs-nucleon coupling and are the hadronic matrix elements.For isospin conserving couplings (G p G n ), the DM-nucleus cross-section in Eq. (3.13) is enhanced by a factor of A 2 .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 N .For this reason, direct detection experiments rely on heavy target materials with large Z to better constrain the DM-nucleon cross-section σ ψN SI .In our model, it is given by where µ ψN = m ψ m N /(m ψ + m N ) is the DM-nucleon reduced mass, m N = 939 MeV and f N = 0.3 [69].Currently, the best upper limits on the SI DM-nucleon cross-section comes from the XENON1T (2018) experiment [65].To constrain the model parameter space from the XENON1T experiment, we use a one-sided Gaussian likelihood function, i.e., we require the parameter points to satisfy5 where σ XENON1T is the 90% C.L. upper limit from the XENON1T experiment and is the effective SI DM-nucleon cross-section.The scaling of σ ψN SI by f rel is done to suppress signals when f rel < 1.In regions of the model parameter space where f rel > 1, 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. [93].

Electroweak baryogenesis (EWBG)
In our model, the VEV of the new scalar S does not initially have to be zero.Thus, the transition pattern can be ( φ , s ) = (0, s i ) → (v, s).At low temperatures, the latter minimum evolves slowly to become the electroweak minimum at T = 0, i.e., ( φ , s ) = (v 0 , s 0 ).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 [14][15][16][17].In addition, it can generate a significant gravitational wave (GW) signal [94].
Using the standard notation, we define a strong first-order phase transition by where v is the Higgs VEV at temperature T .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 v/T [20,[95][96][97][98].
To find regions in the model parameter space where a successful EWBG is viable, we first find the minima of the effective potential V eff (φ, S, T ) (see Appendix D) numerically, and compute the critical temperature T c at which the initial and symmetry breaking minima are degenerate.This allows us to compute the dimensionless parameter v c /T c (the Higgs VEV v c at the critical temperature T c ) and constrain parts of the 7-dimensional model parameter space, i.e., parameter points are excluded if they lead to a too weak phase transition.Specifically, we use a one-sided Gaussian likelihood function and require the parameter points to satisfy Loop functions f T (x) (solid blue) and f S (x) (dashed red).
as a conservative limit.A theoretical uncertainity of 5% on the resulting v c /T c values is assumed to obtain a smooth likelihood function.The actual uncertainty can be much larger as the value of v c /T c required to facilitate EWBG is not yet settled.
In addition, parameter points are also excluded if they exhibit any of the following three features.
1. Incorrect minimum at T = 0: This situation arises when the electroweak vacuum ( φ , S ) = (v 0 , s 0 ) is not the true minimum of the potential at T = 0.
2. Runaway directions in the potential : This occur when the φ and S field values in the symmetric or broken phase are too large, or if the potential in Eq. (2.5) is unbounded from below in the general φ and S directions, i.e., when λ ΦS ≤ −2 λ Φ λ S .
We also perform a complete analysis of the phase transition in this model by following our previous work on the Z 2 symmetric case, i.e., scalar Higgs portal [35] and a very recent update on the calculation of the phase transition dynamics [99].In particular, we find the percolation temperature T p 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.

Electroweak precision observables (EWPO)
With an extra scalar, our model can induce corrections to the gauge boson self-energy diagrams.Its effect on the electroweak precision observables (EWPO) can be parametrised by the oblique parameters S, T and U [100].The γγ and γZ self-energies (Π γγ and Π γZ respectively) are not modified as the new scalar is electrically neutral.Thus, only the W and Z boson self-energies are subject to corrections.
In our model, the oblique parameters are shifted from their SM values by [59] (3.26) These are also plotted in Fig. 2. From Eqs. (3.22)-(3.24), it is evident that and the following correlation matrix To constrain the model parameter space from the EWPO, we use the following likelihood function [102] ln where ∆O i denotes the central values for the shifts in Eq. (3.28), Σ 2 ij ≡ σ i ρ ij σ j is the covariance matrix, ρ ij is the correlation matrix in Eq. (3.29) and σ i are the associated errors in Eq. (3.28).

Higgs searches at colliders
Due to a mixing between the interaction eigenstates (ϕ, s), the coupling strengths between the mass eigenstates (h, H) and SM particles are modified with respect to the SM expectation.The effective squared couplings of (h, H) to SM particles are [63] where X refers to a SM quark, lepton or gauge boson.For the loop-induced processes, these are given by [103] g hYY where YY ∈ (γγ, γZ, gg, ggZ).With modified branching ratios of h/H 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 [67] package.From the model predictions for the two scalar masses, total decay widths, branching ratios, and effective squared couplings defined in Eqs.(3.31) and (3.32), 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 cross-section 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 (h, H), the signal strengths are given by [63] )  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.
To constrain the model parameter space from the Higgs signal strength and mass measurements, we use the HiggsSignals v1.4.0 [68] package.Assuming a Gaussian probability density function (p.d.f.) for the two scalar masses, we compute a chi-square χ 2 HS using the peak-centered method. 6In this method, χ 2 HS 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 [109].Following the assignment, a χ 2 µ 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 mass-resolution such as the h → γγ decay mode), a corresponding χ 2 m is also evaluated by comparing the model predicted and observed Higgs boson mass.Thus, the total χ 2 HS is given by7 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 χ 2 HS .The predicted signal strengths of the two scalars are added incoherently, assuming negligible interference effects.Finally, the computed χ 2 HS is used to define a Higgs signal strength log-likelihood as Thus, a large χ 2 HS indicates a large deviation between the model predicted signal strength and the best-fit value for a fixed Higgs boson mass, and vice versa.

Results
We perform scans of our 7D model parameter space using Diver v1.0.4 [82] with lambdajDE = true, NP = 50,000 and convthresh = 10 −5 .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 high-quality profile likelihood (PL) plots.
We present our model results in the form of 1-and 2-dimensional PL plots.For a model parameter θ i where i = 1, . . ., 7, a 1D PL L PL (θ i ) is defined as Thus, L PL (θ i ) is a function of θ i only, i.e., all other parameters are profiled out.Similarly, a 2D PL L PL (θ i , θ j ) is defined as Thus, L PL (θ i , θ j ) is a function of θ i and θ j only.Using Eqs.(4.1) and (4.2), we can define a PL ratio [110] as where θ ≡ ( θ1 , . . ., θn ) is the best-fit point, i.e., a parameter point that maximises the total likelihood function L(θ).Using Wilks' theorem [111], Eq. ( 4.3) can be used to construct 1σ (2σ) contours corresponding to ∼ 68.3% (95.4%)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 [112] package.

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 v c /T c loglikelihood, i.e., ln L(θ) = ln L vc/Tc (θ), (4.4) where ln L vc/Tc (θ) 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 Λ ≡ L/L max = 1, the dimensionless parameter v c /T c ≥ 0.6 and a successful EWBG is viable.To understand the results in more detail, we go over each panel in Fig. 3 one-by-one.
For a fixed m H and s 0 , Eq. (4.5) has 3 degrees of freedom.As µ 3 , λ S and α are profiled over, it is non-trivial to predict the exact shape of the upper limit on m H as a function of s 0 .The upper limit also weakens as |s 0 | increases.The net result is that for m H 5 TeV, |s 0 | 50 GeV is required to facilitate EWBG.
When α = 0, π, the above condition is satisfied for all values of m H . Thus, a successful EWBG is viable at large values of m H . On the other hand, when α = π/2, Eq. (4.6) imposes the strongest upper limit on m H , namely m H 1.23 TeV.As α → 0, π, the upper limit on m H becomes weaker, as is evident from the plot.
3. (m H , µ 3 ) and (m H , λ S ) planes: In these two planes, all possible combinations of (m H , µ 3 ), (m H , λ S ) and profiled out parameters give v c /T c ≥ 0.6.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 v c /T c likelihood is weakly dependent on µ 3 and λ S as expected from Eq. (3.5).For instance, at large values of µ 3 or λ S which would give |λ ΦS | ≥ 4π or λ ΦS ≤ −2 λ Φ λ S , small values of s 0 can be chosen to avoid such situations.m ψ .This is expected as the contribution from m ψ to the effective potential appears only at 1-loop order.
For m H 5 TeV and m ψ 3.2 TeV, no combination of the profiled out parameters can keep |λ ΦS | < 4π.On the other hand, when m ψ 3.2 TeV, the positive contribution from H to the effective potential can be cancelled out by the negative contribution from ψ, thereby giving v c /T c ≥ 0.6 and Λ = 1.
5. (m H , g S ) plane: For g S 5.62, all values of m H and profiled out parameters give v c /T c ≥ 0.6, and maximise the v c /T c likelihood.However, values of g S > 5.62 lead to runaway directions in the potential as the contribution from g S in the 1-loop 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 m H generally push up the EWSB minimum and cause it to not become the global minimum at T = 0.However, this effect can be counteracted by choosing a large value for m ψ which gives a large negative contribution to the effective potential.One exception is g S > 5.62 which always generates runaway directions in the effective potential.For the remaining model parameters, namely (m H , s 0 , µ 3 , λ S , α, m ψ ), 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.

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 log-likelihood function in Eq. (3.8).In practice, we consider two scenarios in which the fermion DM accounts for either a small fraction (f rel ≤ 1) or all (f rel = 1) of the observed DM abundance.In the former case, we use a relic density likelihood that is one-sided Gaussian, whereas in the latter, we use a Gaussian likelihood.For more details, see subsection 3.1.

Scenario I:
The resulting 2D PL plots from our 7D scans are shown in Fig. 4. For m H m h /2 = 62.6 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 h → HH is kinematically allowed and dominant in this region for all values of the mixing angle α, it reduces the SM-like Higgs signal strength µ h with respect to SM expectation, see Eq. (3.33).This translates into a large χ 2 µ in Eq. (3.36) and is thus disfavoured.
To understand the remaining set of results in more detail, we go over each panel in Fig. 4 one-by-one.
1. (m H , s 0 ) plane: For m H 4 TeV, the parameter planes are ruled out by the EWBG constraint as they either lead to runaway directions, λ ΦS ≤ −2 λ Φ λ S , or nonperturbative couplings, |λ Φ |, |λ ΦS | ≥ 4π.Although, some combinations of the profiled out parameters can give a successful EWBG at large values of m H (see Fig. 3), they are often not compatible with the remaining constraints.This is especially true for the EWPO constraint which only depends on m H and α.For large m H , α 0, π is required in order to satisfy the EWPO constraint.

(m H , α) plane:
We see that the model is allowed by all constraints for a range of low and high m H values provided α 0, π.This is expected as the second scalar H is decoupled in this regime and gives no new contribution to the observed Higgs signal strengths.However, when m H m h , the two scalars are indistinguishable from the point of direct Higgs searches and Higgs signal strength measurements.As is evident in Eqs.(3.14)For m ψ ∈ [m h /2, m H /2], 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 g S which are favoured by the XENON1T limit (as it gives a small DM-nucleon cross-section σ ψN SI ) but disfavoured by the relic density constraint (as it gives f rel > 1) and vice versa.
The diagonal band corresponds to the second resonance m ψ m H /2. Similar to the first resonance m ψ m h /2, all points in this band are allowed by the relic density and XENON1T limits.As g S is profiled over, small values of g S can easily give f rel ≤ 1 and σ eff SI ≤ σ XENON1T .On the other hand, when m H 4 TeV, parameter points are disfavoured by the EWPO and EWBG constraints.For m ψ 3.2 TeV, the region is robustly excluded by the combined constraints.
5. (m H , g S ) plane: In this plane, a lower limit on g S comes from the DM relic density constraint as smaller values of g S lead to an overabundance of the fermion DM in the universe today.This lower limit becomes weaker as m H increases.For m H 4 TeV, the coupling λ ΦS becomes non-perturbative, thus this region is disfavoured.Similarly, values of g S 3.2 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 m H , m ψ and g S . 8From these plots, it is evident that the combined constraints impose an upper and lower limit on m H , m ψ and g S , namely m h /2 m H 4 TeV, 32 GeV m ψ 3.2 TeV, 5.6 × 10 −3 g S 3.5. (4.7) In Fig. 6, we show the key observables such as the fermion DM relic density (left panel) and effective SI DM-nucleon cross-section (right panel).These can be compared against  the Planck measured value and XENON1T limit.It is evident that all of the sampled points satisfy f rel ≤ 1 and σ eff SI ≤ σ XENON1T .We also show the projected sensitivity of the LUX-ZEPLIN (LZ) [113] experiment.Intriguingly, the LZ experiment will probe 2 orders of magnitude smaller DM-nucleon cross sections than the XENON1T experiment.Due to  the two resonances m ψ m h,H /2 and the ability to profile over α, the direct detection cross-section in our model can be significantly suppressed to avoid bounds from current and future direct search experiments.

Scenario II:
In this subsection, we present results from our global fit assuming f rel = 1.The only difference with respect to the f rel ≤ 1 case is our use of a Gaussian likelihood function for the Planck measured DM relic density.In this case, not only small values of g S are disfavoured by the relic density constraint (as they give f rel > 1), large values of g S are also disfavoured (as they give f rel 1).The resulting 2D and 1D PL plots from our 7D scans are shown in Figs.7 and 8 respectively.In comparison with Figs. 4 and 5 respectively, the shape of the 1 and 2σ C.L.  In general, we find that our model can easily satisfy all constraints provided α 0, π.This is expected as the new scalar H is decoupled in this regime and provides no new contribution to the observed Higgs signal strength measurements.Thus, the allowed final states from the fermion DM annihilation are hh, HH and hH, which gives f rel = 1.
An important point to note is that in the Z 2 symmetric case [35], the scalar Higgs portal coupling cannot simultaneously explain the observed DM abundance and matterantimatter asymmetry.This is expected as large values of the portal coupling are required to yield a successful EWBG, whereas small values are required to satisfy the direct detection limits.In contrast, our model contains additional couplings (e.g., µ 3 and µ ΦS ) between the new scalar S and SM Higgs boson; these couplings can aid in generating a strong first-order EWPT.As µ 3 and µ ΦS does not significantly affect the phenomenology of the fermion DM (which is mostly determined by g S and α), the fermion DM can easily saturate the observed DM abundance.These two features together allow the model to simultaneously explain the observed DM abundance and matter-antimatter asymmetry.

Gravitational wave signals
The computation of gravitational wave (GW) signals requires a detailed study of the dynamics of the phase transition (PT).Luckily, the analysis of bubble nucleation is to some extent generic, and the steps required are always similar, albeit using a different scalar potential.In our model, the main difficulty is that the transition always involves both scalar fields, and finding a correct tunnelling path in the 2D field space is always necessary.We tackle this problem in the same way as in Ref. [35].In particular, we use the method described there to find the appropriate tunnelling path and the bubble solutions which drive the transition in each case.The main drawback of this calculation is that it is computationally expensive when compared to all other constraints discussed in section 3. Thus, we first identify interesting points in the model parameter space from our global fit, and check the detailed PT dynamics and GW signals afterwards.
For each viable point from our global fit assuming f rel ≤ 1, we compute the tunnelling path and corresponding action as in Ref. [35].Next, we compute the fraction of volume converted to the true vacuum [99] to accurately calculate the bubble percolation temperature T p .This allows us to identify cases in which too strong supercooling renders percolation impossible; as temperature drops, the false vacuum energy dominates the expansion of the universe and an inflationary phase begins [99,114,115].We find that a significant number of interesting points are excluded as the decay is too suppressed and the transition never successfully finishes.This happens because the extended parameter space with respect to the simple scalar potential in Ref. [35] allows for a formation of a large tree-level barrier which can persist even at T = 0 and suppress the vacuum decay probability.In our 7D scans, we only use the approximation involving the critical temperature T c at which the symmetric and EWSB minima are degenerate.Such points are perfectly valid and predict v c /T c > 1 as required for a successful EWBG.However, after a more detailed analysis, we find that roughly 50% of all points remain viable and their GW spectra have large enough amplitudes to be shown in our plots.
For the viable parameter points, we calculate the ratio of the released latent heat to the energy density of the plasma background, α GW 9 , and the size of bubbles carrying the most energy at percolation R MAX , which we then convert to the more familiar inverse time of the phase transition β/H = (8π) 1/3 v w /(HR MAX ) [99,116].These two parameters are essential for computing the GW spectra [11,94].We also assume that the speed of bubble walls is close to the speed of light (v w /c ≈ 1) which is valid for the very strong first-order PTs that we are mostly interested in.Our calculation of the GW spectrum is based on Ref. [99].In particular, we do not include the signal contribution from collisions of bubbles [117][118][119][120] as the bubbles reach equilibrium with the surrounding plasma and most of the energy is pumped into fluid shells around them [121].The two remaining sources are sound waves in the plasma [122][123][124][125] and turbulence [126][127][128][129] ensuing after the sound waves period.We also check the condition for the sound waves to last more than one Hubble time which was assumed to hold while obtaining the GW spectra in references above.We show this criterion in the (α GW , β/H) plane (see Refs. [99,125] for more details) along with our results in Fig. 9, assuming v w /c ≈ 1 which results in the largest allowed parameter space.We find that no parameter points are consistent with this criterion.This implies that the standard formula for sound wave spectra [11] is probably overestimating the true signal.On the other hand, the turbulence signal will be stronger than the standard estimate as the turbulent motion begins more promptly after the PT.In Fig. 10, we show the resulting GW spectra of viable points as sourced by sound waves (top panel) and turbulence (bottom panel), and their dependence on the percolation temperature T p .As the condition to reliably generate a GW signal from sound waves is not fulfilled, a dedicated numerical simulation would be necessary to ascertain the spectral shape.We expect that the final results will lie somewhere between these two figures.In these figures, we discarded points with almost identical GW parameters to avoid plotting many overlapping results.Thus, we only show 100 representative lines out of 10,000 GW spectra as computed from our results.We also show the detection prospects of Laser Interferometer Space Antenna (LISA) [130], Deci-hertz Interferometer Gravitational-wave Observatory (DECIGO) and Big Bang Observer (BBO) [131].The current and future sensitivity bands of LIGO [132][133][134], the European Pulsar Timing Array (EPTA) [135], the Square Kilometre Array (SKA) [136], Cosmic Explorer (CE) [137] and the Einstein Telescope (ET) [138,139] fall in a different frequency range than that of the viable points.Thus, these experiments give no hope for detection of any of our results.
In summary, we find that the GW spectra of viable points that are interesting from the point of view of baryogenesis can lie within reach of future GW experiments such as LISA, DECIGO and BBO.However, the uncertainty of the sound wave spectrum can have a dramatic impact on the results.In the overly optimistic case of the standard sound wave signal, roughly 15% of all of our viable points would be detectable by LISA while for the most pessimistic turbulence-only spectrum, this number falls below half a percent.We also confirm that pulsar timing arrays and terrestrial experiments (e.g., LIGO) are not sensitive to frequencies that result in a GW signal from a strong EWPT.Notably, our results are qualitatively very similar to the Z 2 symmetric ones in Ref. [99], despite our general non-Z 2 symmetric potential.This leads us to believe that our current knowledge of the Higgs boson properties (most notably, a constraint on the mixing angle α) is enough to significantly constrain viable potentials, and bring them closer to the Z 2 symmetric case.

Conclusions
In this paper, we have performed the most comprehensive and up-to-date study of the extended scalar singlet model with a fermionic DM candidate.After performing a 7D scan Sound Waves  of the model using only the EWBG constraint, we found regions in the model parameter space that can facilitate a successful EWBG.From our 1D PL plots, we showed that a successful EWBG is viable in all parts of the model parameter space provided g S 5.62.
After building intuition from the EWBG only results, we performed a global fit of our model using the constraints from the Planck measured DM relic density, direct detection limits from the XENON1T (2018) experiment, electroweak precision observables (EWPO) and Higgs searches at colliders.This allowed us to constrain parts of the 7D model parameter space.In particular, our global fit placed an upper and lower limit on m H , m ψ and g S , namely m h /2 m H 5 TeV, 32 GeV m ψ 3.2 TeV and 5.6 × 10 −3 g S 3.2 GeV.Moreover, we confirmed that our model can simultaneously yield a strong first-order phase transition and saturate the observed DM abundance.This is an important feature which is missing in the Z 2 symmetric case.
From the viable points that satisfied all of the available constraints, we computed the GW spectra, and checked the discovery prospects of the model at current and future GW experiments.In doing so, we found that the GW spectra of viable points can be within reach of future GW experiments such as LISA, DECIGO and BBO.We checked that the condition for sound waves to be a long-lasting source of GWs is not satisfied for any of our results.This implies that the standard sound wave spectrum used in the literature likely overestimates the true signal, whereas the turbulence signal can be stronger than the standard prediction as the turbulence sets in quicker after the end of the phase transition.Unfortunately, the overall result will still likely be a reduction of the overall spectrum, thereby reducing the discovery prospects.Specifically, in our results we find that 15% of our viable points would be within the reach of LISA if the final spectrum was close to the standard sound wave prediction.However, this number falls down to less than half a percent in the most pessimistic case of only a turbulence-sourced GW signal.
After EWSB, the φ and S fields acquire their VEVs in Eq. (2.8).Thus, the following partial derivatives ∂V tree ∂G 0 , ∂V tree ∂G − , ∂V tree ∂G + , ∂V tree ∂φ , ∂V tree ∂S , must vanish at the EWSB minimum ( φ | T =0 , S | T =0 ) = (v 0 , s 0 ).This gives A simple rearrangement gives us the following EWSB conditions Now, we compute the second-order partial derivatives at the EWSB minimum.The only non-zero ones are given by Using Eqs.(A.2) and (A.3), these expressions can be simplified to After EWSB, the φ and S fields can be expanded as As ∂V tree /∂φ = ∂V tree /∂ϕ and ∂V tree /∂S = ∂V tree /∂s, a mass-term for the real scalar fields where is the squared mass matrix.Using Eqs.(A.5)-(A.7),the matrix elements are given by (A.11)For the EWSB minimum to be a stable (i.e., not a saddle point) solution of Eq. (A.1), the symmetric 5 × 5 Hessian matrix H must be positive-definite.At the EWSB minimum, it is given by .
The above matrix is guaranteed to be positive-definite if the determinant (eigenvalue) of the 2 × 2 sub-matrix is non-zero (positive).These two requirements give To study the bounds of the tree-level scalar potential, Eq. (A.1) can be expressed in terms of the φ and S fields as Depending on the chosen direction in the (φ, S) plane, three scenarios are possible.
1. Pure φ direction: In this case, the potential depends only on the φ field.It is bounded from below provided λ Φ > 0.
2. Pure S direction: In this case, the potential depends only on the S field.It is bounded from below provided λ S > 0.
3. General φ and S directions: At large φ and S field values, the quartic terms in Eq. (A.13) dominate.In this case, the potential can be approximated by Thus, the potential is bounded from below provided λ S > 0 and λ ΦS > −2 λ Φ λ S .

B Mass eigenstate basis
The squared mass matrix M 2 is real and symmetric.It can be diagonalised by an orthogonal matrix O. Thus, we define the mass eigenstates (h, H) as The interaction eigenstates (ϕ, s) are given by Now, we consider the following matrix product where D = diag(m 2 h , m 2 H ) is a diagonal squared mass matrix for the physical mass eigenstates.Thus, the last equality in Eq. (B.2) requires The left-hand side of the above expression expands to By equating these expressions to the elements of the diagonal matrix D, we get By computing the inverse of the above 3 × 3 matrix (i.e., by taking α → −α), we get From the above matrix product and Eq.(A.11), we get

C Dark matter-nucleon coupling
The interaction eigenstates (ϕ, s) can be written in terms of the mass eigenstates (h, H) as ϕ s = cos α sin α − sin α cos α h H .
In a typical direct detection experiment, the momentum transfer q is roughly on the order of a few MeV.Assuming that the mediator masses m h/H are well above this value,10 i.e., m 2 h/H q 2 , we can safely approach direct detection in the context of an effective field theory (EFT) and integrate out the scalar mediators [141].Thus, we can write down an effective DM-quark interaction Lagrangian as where is the effective DM-quark coupling.
In order to promote a DM-quark interaction to a DM-nucleon one, the quark contents of a nucleon must be taken into account.For a scalar mediator (as in our model), the quark Yukawa couplings generally scales with the mass of an interacting fermion.Thus, the dominant contribution comes from the strange quark content of a nucleon and from gluons via heavy quark loops.These contributions are parametrised by the hadronic matrix elements as f (N ) where N ∈ (p, n).For a pure scalar interaction, these matrix elements parametrise the contribution of a quark mass m q to the total mass of a nucleon m N .For more details on these parameters and recent estimates, see Ref. [93] and references therein.Using the heavy quark expansion [142], the contribution from gluons via heavy quark loops can be expressed in terms of the lighter quarks as G q m q f is the Higgs-nucleon coupling [93].Thus, the effective DM-nucleon interaction Lagrangian can be written as where is the effective DM-nucleon coupling [84].
For a SI DM-nucleon interaction, the DM-nucleus interaction is a coherent sum over the total number of protons Z and neutrons (A − Z) in the target nucleus N .Thus, the SI DM-nucleus cross-section is given by

D Effective potential
We include the following 1-loop corrections to the zero temperature potential in the cutoff regularisation and on-shell scheme [20,143] V 1-loop (φ, S) = W, Z, t, ψ i=φ, S, χ where n {φ, S, χ, W, Z, t, ψ} = {1, 1, 3, 6, 3, −12, −4}.The subscript "0" implies that the particle masses are calculated at the T = 0 minimum, i.e., (φ, S) = (v 0 , s 0 ).The φ and S field dependent masses are given in Appendix B, whereas the rest are given by The final important correction comes from resumming the multi-loop contributions to the boson longitudinal polarizations which are infrared divergent [18,144].These are incorporate by supplementing the scalars and longitudinal polarizations of the gauge bosons with thermal mass corrections, in particular, by expanding Eq. (D.3) to the leading order in m 2 /T 2 [18].For our model, they are given by Π φ (T ) = Π χ (T ) = T For the φ and S fields, the corrected masses are the eigenvalues of the following squared mass matrix where M 2 is defined in Eq. (A.10).For the Z and γ fields, namely m 2 Z/γ + Π Z/γ (T ), the mass corrections are the eigenvalues of the following squared mass matrix    In other cases, we simply use the following substitution Finally, the effective potential V eff (φ, S, T ) is given by V eff (φ, S, T ) = V tree (φ, S) + V 1-loop (φ, S) + V T (φ, S, T ), (D.9) where V tree (φ, S) is the tree-level scalar potential in Appendix A.

1 .
(m H , s 0 ) plane: For m H 1.3 TeV, all values of s 0 and some combination of 5 profiled out parameters (namely µ 3 , λ S , α, m ψ and g S ) give v c /T c ≥ 0.6 and maximise the v c /T c log-likelihood, thus Λ = 1 everywhere.Due to the dependence of s 0 in Eq.(3.5), large values of |s 0 | should lead to runaway directions, λ ΦS ≤ −2 λ Φ λ S , and/or nonperturbative coupling, |λ ΦS | ≥ 4π.With α = π/2, a large contribution from m H to λ ΦS can be suppressed.However, this choice of α makes λ Φ in Eq. (3.3) nonperturbative as its contribution appears as m 2 H sin 2 α.Ultimately, the solution is to choose a small value for λ S as its contribution in Eq. (3.5) appears as −λ S s 2 0 .In addition, small values of µ 3 can also help in keeping |λ ΦS | < 4π.Thus, for m H 1.3 TeV, large values of |s 0 | can facilitate EWBG.For m H 1.3 TeV and |s 0 | 50 GeV, the white region (Λ = 0) is disfavoured as it leads to |λ ΦS | ≥ 4π.This is expected as the contribution from m H in Eq. (3.5) is dominant at large values.With large |s 0 |, no choice of µ 3 , λ S and α can keep |λ ΦS | < 4π.In fact, the requirement |λ ΦS | < 4π translates into an upper limit on m H as a function of s 0 , µ 3 , λ S and α.Using Eq. (3.5), we get

2 . 1 . 3
(m H , α) plane: Similar to the (m H , s 0 ) plane for m H 1.3 TeV, some combination of the profiled out parameters gives v c /T c ≥ 0.6 for all values of α.However, when m H TeV and α = 0, π, the Higgs quartic coupling λ Φ in Eq. (3.3) becomes non-perturbative.In fact, the requirement |λ Φ | < 4π translates into the following upper limit on m H as a function of α

Figure 3 .
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 1σ (68.3%) and 2σ (95.4%)C.L. regions.In regions where Λ ≡ L/L max = 1, a successful EWBG is viable as v c /T c ≥ 0.6 (see text for more details).The parameter planes (m H , µ 3 ) and (m H , λ S ) are not shown as they are unconstrained by the EWBG constraint.

Figure 5 .Figure 6 .
Figure 5. 1D PL plots for m H (left), m (center) and g S (right) assuming f rel ≤ 1.The respective plots for s 0 , µ 3 and λ S are not shown as they are unconstrained by our global fit.

Figure 8 .
Figure 8. Same as Fig. 5 except for the f rel = 1 case.

Figure 9 .
Figure 9. Viable points from our global fit assuming f rel ≤ 1 in terms of the parameters α GW and β/H.The grey region shows where the sound waves last more than a Hubble time (assuming v w /c ≈ 1 which results in the largest allowed area) and reliably produce a GW signal.

Figure 10 .
Figure 10.Gravitational wave (GW) spectra of viable points as sourced by sound waves (top panel) and turbulence (bottom panel) along with their dependence on the percolation temperature T p .Current sensitivity bands of LIGO and EPTA, as well as detection prospects of LISA, DECIGO, BBO and SKA are also shown for comparison (see text for more details).

1 2 λ
ΦS S 2 , m ψ = µ ψ + g S S. (D.2)The finite temperature corrections to the effective potential are given by V T (φ, S, T ) = n i

Table 1 .
Ranges and priors for the free parameters of our model.All parameters are uniformly sampled over their ranges in either flat or logarithmic space.For the mixing angle α, all constraints are symmetric under α → −α, thus we only scan over α ∈ [0, π].
However, when these decay modes are kinematically allowed, they suppress the h/H signal
and (3.27), direct detection and EWPO constraints respectively are also relaxed in this regime.The net result is that all values of α are allowed when m H m h .3.(mH , µ 3 ) and (m H , λ S ) planes: These parameter planes are mostly unconstrained by our global fit except for m H m h /2 (excluded by the Higgs signal strength measurements) and m H 4 TeV (ruled out by the EWBG constraint).For m H 1.3 TeV, large values of λ S are required to facilitate a successful EWBG.4. (m H , m ψ ) plane: For m ψ m h /2, the fermion DM can only annihilate into light SM quarks, thereby giving f rel > 1.On the other hand, m ψ m h /2 is constrained by the DM relic density and XENON1T limits.When m ψ m h /2, all values of m H up to ∼ 4 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 g S can yield a fermion DM relic density and DM-nucleon cross-section that is compatible with the Planck measured value and XENON1T limit respectively.