Dark matter and electroweak phase transition in the mixed scalar dark matter model

We study the electroweak phase transition in the framework of scalar singlet-doublet mixed dark matter model, where the particle dark matter candidate is the lightest neutral Higgs that comprised of the CP-even component of inert doublet and a singlet scalar. The dark matter can be dominated by the inert doublet or singlet scalar depending on the mixing. We present several benchmark models to investigate the two situations after imposing several theoretical and experimental constraints. An additional singlet scalar and the inert doublet drive the electroweak phase transition to be strongly first order. A strong first order electroweak phase transition and a viable dark matter candidate can be accomplished in two benchmark models simultaneously, where a proper mass splitting among the neutral and charged Higgs masses is needed.


I. INTRODUCTION
The dark matter (DM) relic abundance and baryon asymmetry of the universe (BAU) are two fundamental problems puzzling the cosmologists and physicists. With accumulating of direct and indirect detection experiments, the weakly interacting massive particle (WIMP) is still one of the most popular DM candidates, which requires new physics beyond the standard model of particle physics. On the other hand, the electroweak baryogenesis (EWBG) provides a very attractive mechanism to explain the BAU problem. One key ingredient of EWBG mechanism is the strong first order electroweak phase transition (EWPT), which prevents the baryon asymmetry from being washed out by the sphaleron process [1][2][3] and helps to understand the electroweak symmetry breaking. The model that can provide a strong first order EWPT might drive deviation of the triple Higgs coupling from the SM prediction, which can be detected at future high energy proton-proton collider [3,4]. It could be very intriguing if new physics can explain DM relic abundance and accomplish a strong first order EWPT.
In the simplest scalar extended dark matter models, the DM-Higgs quartic coupling leads to the interaction between the DM and nuclear. Except the so-called Higgs funnel regime (with m DM ∼ m h /2), the scalar singlet dark matter suffers severely constraints from the upgrade of direct detection experiments, such as PandaX-II [5], LUX [6], and XENON [7], etc. Another popular DM model is the Inert Doublet Model (IDM), where only one doublet receives a VEV in the two-Higgs doublet model. The IDM has been studied extensively as a framework to provide a scalar doublet DM [8][9][10][11][12][13][14][15][16][17][18]. In this model, the low mass and intermediate mass regions are almost excluded by the collider constraints and the DM direct detection limits, the remnant space is the high mass region ( > 500 GeV) which is quite difficult to probe in the direct detection experiments [19]. Then it is very natural to consider if extending of such models can offer a possibility to avoid current and even future direct detection bounds. DM phenomenology aside, The strength of the EWPT can generally be enhanced by adding a scalar to the SM Higgs potential through cubic or quartic couplings with the SM-like Higgs field. The role played by the scalar dark matter during the universe cooling down has been studied extensively, including singlet as well as inert doublet types DM, e.g. see Ref. [20][21][22][23][24][25][26][27][28][29][30][31].
A scenario with the DM being the singlet-doublet mixed scalar is very appealing, which has been studied previously [32][33][34][35][36]83] 1 . In such a model, the IDM is extended with a SM singlet scalar. The inert doublet and the singlet are odd in a new Z 2 discrete symmetry whereas the SM sector is even. After electroweak symmetry breaking, the electroweak charged doublet and the extra singlet mix into the scalar DM. The mixing comes from the renormalizable couplings to the SM-like Higgs field. The DM-Higgs boson triple and quartic couplings in such model are controlled by the singlet-doublet scalars mixing. As a consequence, the DM annihilation and spin-independent DM-nucleon scattering processes can be characterized by the mixing. The thermal DM relic abundance can be achieved with an appropriate magnitude of the mixing. The parameter space of this kind of model has the advantage of less constrained by the spin-independent direct detection experiments.
In this work, we study dark matter phenomenology and EWPT in the scalar singlet-doublet mixed DM model. Different from Refs. [34,35], although the DM self-interactions have no impact on the DM annihilation and DM-nuclear scattering, we keep them in our model since these terms might affect the vacuum barrier during the Universe cooling down. In the model, besides the contribution of the inert doublet to the vacuum barrier, there is an additional contribution comes from the singlet scalar which is also determined by the mixing. Therefore, a proper mixture also characterizes the realization of the strong first order EWPT.
In this model, we find that the additional singlet along with the mixing effects make the situation of EWPT different from the IDM, two benchmarks with moderate dark matter mass are allowed to gain the correct relic abundance and to achieve a strong first order EWPT. Some features of the model are found and listed bellow: 1. The vacuum stability and T parameter constrain the singlet-doublet mixing non-trivially.
2. DM relic density and the DM-nuclear scattering cross section highly depend on the mixing of the singlet and doublet scalars.
3. Co-annihilation effects play an important role near the degenerate masses region of the dark sector, which lead to the cancellations between the annihilation processes to yield correct DM relic density.

4.
A proper mass splitting between the Z 2 odd neutral Higgs and charged Higgs is needed to obtain a strong first order EWPT.
The paper is organized as follows. In Section II, we present the mixed scalar dark matter model. Theoretical constraints and LEP bounds on scalar masses are addressed in Section III. The Section IV is devoted to the comments on DM relic density and direct detection in the model. In the Section V, we present the approach of estimation of the strength of EWPT with thermal effective potential. In Section VI, we give a comprehensive analysis of DM phenomenology and EWPT in several benchmark models. After that, we address the collider prospects of the benchmark models in Section VII. We conclude in the Section VIII.

II. THE MIXED SINGLET-DOUBLET SCALAR DARK MATTER MODEL
We extend the Inert Doublet Model with an additional real singlet scalar (S ) and assume the Z 2 symmetry to inert doublet Φ I and the S , the tree level potential of the model is where Φ H is the SM Higgs doublet and Φ I is the inert doublet. The doublets are given as where Φ H develops a VEV v = 246 GeV whereas the inert doublet and the real singlet scalar S do not produce any VEV. The Z 2 symmetry of Φ I , S remains unbroken. Goldstones G + and G 0 are absorbed in W ± , Z bosons after spontaneous symmetry breaking. All the parameters in Eq. 1 are assumed to be real. By minimizing the potential in Eq. 1 we obtain the relation and the mass terms for the SM Higgs and the charged scalars can be obtained as Working in the basis of real neutral scalars, (S , H 0 , A), where H 0 and A are the real and imaginary components of the neutral component of Φ I . The mass squared matrix is Considering one of the mixed states of the CP-even scalars S and H 0 is the lightest particle among the three eigenstates, then it could serve as the DM candidate, and the DM mixing is induced by the upper left 2 × 2 block of the mass matrix Eq. 5.

III. THEORETICAL AND EXPERIMENTAL CONSTRAINTS
A. Perturbative unitarity In the high energy limit, the quartic contact interaction terms contribute to the tree-level scalar-scalar scattering matrix dominantly. The s-wave scattering amplitudes are constrained by the perturbative unitarity limits, which requires that the eigenvalues of the S -matrix M must be smaller than the unitarity bound: The perturbative unitarity of the general two-Higgs-Doublet model were first studied in Ref. [41,42]. In this work we extend the formalism of Ref. [43,44] for the IDM to states containing an extra singlet S . The initial states are classified according to their total hypercharge Y (0, 1 or 2), weak isospin σ (0, 1 2 or 1) and discrete Z 2 charge X.
For simplicity, we only list here the initial states with hypercharge Y = 0 and σ = 0 which differ from the 2HDM initial states [43,44]: where the first three states are even under Z 2 and the last two states are odd. In the following, we present all the scattering matrices of the model: Then the eigenvalues Λ X Yσi of the above scattering matrices (where i = ± or 1, 2, 3) can be calculated as and the Λ even 00 1,2,3 correspond to the three roots of the polynomial equation

B. Vacuum stability conditions
To make sure the potential is bounded from below in the limit of large field values, we impose copositivity criteria on the quartic couplings of the tree level potential as explored in Ref. [45]. The vacuum stability conditions for the model are given as and for λ 4 > 0 case, if λ 4 < 0, one obtains λ 1 λ 2 λ 6 + λ 1 λ 8 + λ 2 λ 7 + λ 6 (λ 3 + λ 4 − |λ 5 |) It should be noted that, with loop corrections being included, one may expect more viable parameter spaces allowed by the vacuum stability conditions.

C. Electroweak precision test
We should also make sure that electroweak precision test data is respected when we study the phenomenology of the model. As we know, electroweak precision test can be expressed in terms of three measurable quantities, called S , T and U, which parameterize the contributions from beyond standard model physics to electroweak radiative corrections [46]. The most relevant one is the T parameter, which characterize the spoil extent of the global S U(2) symmetry. Since the scalar singlet S mixes with the H 0 in this model, the self-energies of the Goldstone bosons get corrected by the diagrams with virtual inert particle pairs and S . Following the method in Ref. [10], we first obtain In the limit of custodial SU(2) symmetry, i.e., m H ± = m A , Eq. 31 can be recast as which sets bounds on the magnitude of the mixing between S and H 0 , as well as the mass splitting among neutral and charged Higgses in the model. The T parameter can be obtained as T = ∆ρ/α EM following the notation of Ref. [10]. The Gfitter fit to the electroweak data [47]: T = 0.09 ± 0.13 has been used to constrain the parameters.

D. LEP bounds
The LEP bounds on the scalar masses of this model can be considered as the same as in IDM, and the only difference is that the CP-even neutral Higgs H is an admixture of doublet and singlet. The precise measurements of the W and Z widths lead to the following lower limits on the scalar masses to ensure that decay channels of W ± → H/χH ± , AH ± and Z → H/χA, H + H − are kinematically forbidden. The production of the charge Higgs pairs e + e − → H + H − at LEP-II [48] sets which does not apply when the scalar mass is larger than m Z /2 providing the mass splitting m H ± − m H is smaller than 5 GeV [49].
For the neutral Higgses of the model, the constraints come from the pair production process of e + e − → H/χA followed by the cascade decay A → H/χZ → H/χ ff , which can be obtained through a reinterpretation of the neutralino production search at LEP-II. The analysis of Ref. [50] shows that the limit on the neutral Higgs is max(m A ,m H )≥ 100 GeV for the IDM, which can be roughly applied to our model when the mixing of H 0 and S is small, or can be interpreted as max(m A ,m χ )≥ 100 GeV when the mixing of H 0 and S is large, which can be figured out from Eq. 53 in the Section VII D .

IV. DM RELIC DENSITY AND DIRECT DETECTION
The thermal DM relic density should match the observed DM density in the Universe of Ωh 2 = 0.1199 ± 0.0022 [51], or at least of not over-producing such density via thermal production. To an approximation, the relic density is given by where x f = m D /T f 20 with T f being the typical freeze-out temperature of a WIMP [52], m D = m χ or m H is the DM mass, M Pl is the Planck mass, g * is number of relativistic degrees of freedom, σv is the thermally averaged cross section for DM pair annihilation into the SM particles (i.e. ff , W + W − , ZZ, hh). For the mixed DM scenario, the DM annihilation is characterized by the masses of dark matter and the mixing of DM constituents. To illustrate the DM annihilation properties with different channels in this model, we collect several benchmarks in distinct DM mass regions in the following subsection. The relic density is calculated by implementing package micrOMEGAs [53].
In the direct detection experiment of dark matter, the nuclear recoil spectrum is directly related to the DM-nuclei scattering cross-section [54], which is given by: where q is the momentum transfer, µ r = (m nucl m D )/(m nucl + m D ) and v is the relative velocity. The couplings of DM to the proton and neutron, f p and f n , can be expressed as where m N is the nucleon mass, f N T q is the form factor of the nucleon (see Table I) and λ DDqq is the effective coupling of the DM particle to a q-flavor quark component in the nucleon. In this model, DM-quark interaction derives from t-channel exchange of the SM Higgs h. Thus, in the limit of zero momentum transfer, the Higgs can be integrated out and the effective coupling becomes with the coupling a hDD given in Eq. 15 and Eq. 16. The direct detection rates in our calculation also have been evaluated using package micrOMEGAs.  Direct searches for dark matter by the LUX [6] and PandaX-II [5] Collaborations 2 have recently come up with the most stringent limits on the spin-independent elastic scattering of DM off nucleons. In comparing the calculated scattering cross section to the limits from the experiments, we rescale the cross section by the ratio of the predicted and observed relic density to account for the reduced flux of dark matter particles in the detectors when the relic density is undersaturated, i.e., σ scaled = σ D−nucl · Ω DM h 2 /Ωh 2 , with Ωh 2 being the observed relic density, we take the central value 0.1199 in our calculations.

V. RESEARCH STRATEGY OF ELECTROWEAK PHASE TRANSITION
The strength of electroweak phase transition can be parametrized by v c /T c , and generally the strong first order condition v c /T c > 1 is needed to prevent the washout of baryon asymmetry [55] 3 . For more details on the electroweak baryogengesis, we refer to Ref. [3], and the interplay between the v c /T c and the generation of baryon asymmetry in the framework of EWBG can be found in Ref. [56].
The effective potential at finite temperature can be used to explore the Universe cooling history. To perform the numerical analysis of vacuum structure at finite temperature T , one needs the effective potential 7 with V 0 , V CW and V T being the tree-level, one-loop temperature-independent and -dependent Coleman-Weinberg potentials, respectively. The tree-level potential V 0 could be obtained from Eq. 1 after the field expansion. There could be multi-step phase transition depends on the vacuum structure [56,[58][59][60]. In our case, we focus on the one-step phase transition for simplicity. More complicate scenarios are left to future studies. In the Landau gauge (ξ = 0) and MS scheme, the temperature-independent one-loop Coleman-Weinberg potential is given by [61][62][63][64]: Here, all fields coupling to the Higgs are summed, n i is the number of degrees of freedom for bosonic and fermionic fields, and renormalization-scheme-dependent constants C i = 1/2 for transverse gauge bosons and 3/2 for for longitudinal polarizations of gauge bosons and other particles; m 2 i (h) are the field-dependent squared masses for all particles. Here, the counterterms have been absorbed into V CW implicitly.
The temperature-dependent effective potential is obtained as [64,65] with the J functions defined as To accomplish the analysis of the phase structure as a function of T analytically, high-temperature expansions are always used: where a b = 16a f = 16π 2 exp(3/2 − 2γ E ). It should be noted that the T 2 terms in the above expressions can drive symmetry restoration at high temperatures, and the non-analytic (m 2 ) 3/2 term in Eq. 45 is responsible for the barrier between the symmetric (at high temperature) and broken phases(at the critical temperature). To avoid the breakdown of perturbation theory induced by symmetry restoration at high temperatures, one needs to perform a resummation of daisy diagrams [66] by adding finitetemperature corrections to the boson masses in Eq. 43: with m(T ) 2 computed from the infrared limit of the corresponding two-point function [67] and are presented in the Appendix A, this approach is usually called Parwani's method in literatures.
Here, the thermal masses contribution of m S S (T ) gives an additional contribution to the (m 2 ) 3/2 term 4 of Eq. 45 in comparison with the IDM case. In Fig. 1 we depict the fraction m S S (T ) 3/2 / i (m 2 i ) 3/2 to illustrate the contribution in various parameter spaces. We find that the ratio is more than 10%, from which one could expect the vacuum barrier can be altered. Therefore the mixed dark matter and the mixing of CP-even neutral scalars might affect the universe cooling process, which will be explored in the following section.

VI. DM AND EWPT PHENOMENOLOGY
In the low mass region of dark matter, the Ref. [25][26][27] found the dark matter relic abundance and the strong first order EWPT can be accomplished in the Higgs funnel regime in the IDM. The bounds of the current direct detection on the high mass region is much weaker, which is caused by the lower number density for high mass DM. Therefore, we alternatively focus on the intermediate dark matter mass region. We investigate if the contributions of dark matter and the other CP-even scalar to the vacuum barrier help to realize the strong first order EWPT.

A. Benchmarks for DM from different mass regions
In this subsection, we present several benchmarks for distinct DM mass regions. After imposing the Higgs collider bounds by using HiggsBounds [68][69][70], and checking rate measurements with HiggsSignals [71], all the parameter setting in each benchmark model agree with the constraints aforementioned. The quartic coupling λ 2 , λ 6 , λ 8 are not involved in the DM annihilation or DM-nucleus scattering processes, so we fix them only by following the constraints from the unitarity and vacuum stability conditions.

The original IDM case
At the beginning, we briefly study the DM phenomenology in the original IDM for comparison. By setting λ 6 = λ 7 = λ 8 = θ = 0, we restore to the IDM scenario. Two benchmark points are shown in Table II.  Firstly, we take m H to be smaller than the gauge bosons, the main DM annihilation channel is HH → h → bb. To avoid an overlarge cross section which leads to the relic under-abundance, we expect the DM-Higgs coupling a hHH = 2vλ L to be small enough. As the DM mass increasing, a much smaller λ L is required, as could be found from IDM-BP1 and IDM-BP2 models where the DM masses are in the Higgs funnel region. When the mass of the DM becomes larger than the gauge boson's, the annihilations to the bosons are much efficient, the DM is under-abundant even the coupling λ L → 0. Then the degenerate case is the only viable scenario, for example, a parameter set can be m H = 550GeV, m A = m H ± = 555GeV, and also λ L → 0 5 .

χ as the DM candidate
In the singlet-doublet scalar mixed model, we first investigate the possibility of χ being the DM candidate. We identify two benchmarks given in Table III, where the correct DM relic density can be matched and the magnitude of the cross sections of scattering off nucleus are below the LUX (or PandaX-II) strictly upper limits. When m χ < m V (V = W ± , Z), the dominant annihilation channel is χχ → h → bb, then the DM-Higgs coupling a hχχ characterizes the relic density as well as the DM-nucleus scattering. For m χ > m V (BP2), the annihilations to gauge bosons could be dominant which lead to the DM under-abundant (as discussed in the original IDM case). Due to the DM-gauge boson quartic couplings λ χχVV are proportional to sin 2 θ, we need a smaller θ in order to avoid the larger cross section of χχ → VV. Above the hh threshold with a small θ, the channel χχ → hh is the dominant one.     Table IV, where χ is the DM candidate. Here, considering the stability constraints, the λ 2 that does not take part in the DM phenomenology is taken to be a larger value in comparision with other benchmark scenarios in Table III and Table V. For the BP3, the Higgs invisible decay channel is open, the masses of both H and χ are smaller than m h /2 in comparison with the BP1, three decay channels h → χχ, h → Hχ, h → HH should be considered. For the BP4 and BP5, the decay channel H → χh is involved in the decay chains of e + e − → HA. We will discuss the constraints from collider experiments on parameter space and possible signals searching at colliders in the Sec. VII.

H as the DM candidate
For a comparison between the scenarios of χ and H being the DM candidates, here we investigate the case of m H < m χ by employing the benchmarks in Table V. The anatomy of dark matter phenomenology in different benchmark models are discussed in the following.  We exchange the two scalars' masses in comparison with BP1. The DM-Higgs triple coupling (here it is a hHH ) plays an important role in the DM relic density estimation since the annihilations to the gauge bosons are kinematically forbidden. A simply exchange of λ L and λ 7 compared to BP1 would satisfy the relic abundance and the direct detection cross section limits of DM, as indicated by Eq. 14, Eq. 15 and Eq. 16.

BP7
When m H > m V , the pair production of WW(ZZ) in the final state of DM pair annihilation processes becomes increasingly important, because the four-point interaction through gauge couplings and s-channel process start contributing significantly. As a result the cross section for DM pair annihilating to gauge bosons becomes very large, such that the estimated thermal DM relic density may be systematically below the observed DM relic density with combinations of model parameters.
The DM-gauge boson quartic couplings λ HHVV are proportional to cos 2 θ, so a smaller θ drives much more efficient annihilations to the gauge bosons and therefore DM relic density is not abundant. We take a θ = 1.1 in the benchmark model BP7. The main component of the DM is still the singlet scalar, which is the same as the BP2 by taking the χ to be the DM candidate.

Degenerate case
In the IDM scenario, when scalars H, A and H ± are nearly mass-degenerate, the co-annihilation processes become important. There is a cancellation taking place between the t/u channel contributions and the four-vertex diagram [10,19]. In this manner, the WIMP depletion rate can be balanced by varying the mass splitting between DM and other scalars and the parameter λ L to obtain the correct mixture of transverse and longitudinal gauge bosons in the final state. These solutions are always found for small mass splittings of the odd particles, and require some tuning of the value of λ L . In practice the maximal allowed mass splitting is of the order O(10) GeV [19].
In benchmark models of BP-c1 and BP-c2, the cancellation also happens augmented by the mixing between the singlet and doublet components of the DM particle. For benchmark BP-c2, the difference from the BP-c1 case is the odd particle χ do not degenerate with H. The mass splitting is much bigger which varies µ soft dramatically, so the λ L should be tuned to adjust the DM-Higgs couping a hhH to be smaller enough, which is required to evade the direct detection limits.

B. Electroweak phase transition assisted by the mixed dark matter
We exam all the benchmarks given above to find which can yield a strong first order EWPT. There are two benchmark models left which are summarized in Table VI along with the corresponding critical temperatures and the strengths of the phase transition. In the following we present the combined results of both DM and EWPT by varying sensitive parameters near the benchmark points BP5 and BP-c2. The χ − H mixing effects are shown in Fig. 3 and Fig. 5, and other parameters effects are presented in Fig. 4 and Fig. 6. Before devoting to the numerical analysis of DM and EWPT, we firstly consider the electroweak precision test constraints. The T parameter exclusion of two benchmark scenarios are shown in Fig. 2. The two plots in top panel imply that the mixing angle should be small for the DM mass around 200 GeV. For BP-c2, parameter space is less constrained by the T parameter. We also perform the constraints from the stability conditions, the blue region in Fig. 3, Fig. 4, Fig. 5, and Fig. 6 are excluded. One can see that the parameter space is constrained severely. Moreover, a large variation of the quartic coupling λ 2 , λ 6 , λ 8 only leads to tinny affect on the one-step phase transition case as numerically checked by us. Therefore we focus on the survey of the phase transition in the parameter spaces of λ L , λ 7 , θ, m H,χ with others being fixed.

EWPT
As aforementioned in Sec. V, the vacuum barrier can be affected by the DM mass slightly. We find that the DM mass dependence of the strength of EWPT is tinny in this benchmark as explored in Ref. [56], and the whole mass range can provide a strongly first-order EWPT (red-dotted region), as depicted in the parameter spaces of m χ − θ (bottom-right of Fig. 3) and m χ −λ 7 (bottom-right of Fig. 4). We observe that the strength of the EWPT subtly depends on the mass splitting ∆M = |m H − m A,H ± | and the mixing angle, and the larger θ leads to a strong first order EWPT when the mass splitting ∆M is large. As depicted in the top-left panel of Fig. 3 and top panels of Fig. 4, the realization of the strongly first-order EWPT prefers a narrow range of m H with the fixed values of θ, λ 7 . Furthermore, a limited range of the coupling λ L is required by the condition v c /T > 1 as shown in top-right panel of Fig. 3 and bottom-left panel of Fig. 4.

DM relic density
In this benchmark model, the scalar χ plays the role of DM. As demonstrated in the λ L − θ, λ 7 − θ, m χ − θ planes of Fig. 3, the mixing angle θ induces efficient change of the magnitude of the DM relic density, this is because the main annihilation channel is χχ → hh with a θ-dependent four-vertex coupling 2λ 7 cos 2 θ + 2λ L sin 2 θ. We find that the DM relic requires mixing angle to be smaller than 0.5, and then we notice that the main contribution comes from 2λ 7 cos 2 θ part. Thus the DM relic density sensitively depends on the parameters λ 7 .
In the top-left panel of Fig. 3 and top panels of Fig. 4, when m H becomes larger than about 400 GeV, the main annihilation channel changes to be χχ → tt which is proportional to the DM-Higgs coupling a hχχ . Due to the coupling a hχχ can be affected a lot by the m H since which is involved in the parameter µ soft , therefore, one can expect the predicted DM relic abundance is affected by the variation of m H . In the top-left panel of Fig. 3 (top panels of Fig. 4), with the increasing of m H one gets an increasing (decreasing) a hχχ , and therefore a narrower (broader) viable region for the correct DM relic density.

DM direct detection
Considering the current restrictive limit from LUX (and PandaX-II), we study the constraints on the parameter space of the model. The gray regions in the figures are excluded by the DM direct detection limits.
In the intermediate mass region (around 200 GeV), we find that LUX limits constrain the parameter θ and λ 7 rigorously, due to the DM-Higgs coupling a hχχ ∼ 2vλ 7 cos 2 θ. Only a small region evades the constraint in each plane of Fig. 3 and Fig. 4. The narrow allowed parameter spaces are due to the cancellation effects in the DM-Higgs coupling a hχχ , as can be seen from Eq. 15 when m χ < m H . In the bottom-right plots of Fig. 3 and Fig. 4, we preform the exclusion of θ, λ 7 by varying the DM mass. Although the limit is much strict, there is still room for obtaining the correct DM relic density and the strong first order EWPT.

H-DM benchmark model: BP-c2
The combined results of Benchmark BP-c2 are presented in the Fig. 5 and Fig. 6. In this benchmark model, the scalars H, A and H ± are nearly mass-degenerate.

EWPT
To obtain the correct DM relic density the scalars in the inert doublet should be degenerate when the DM mass m H are larger than the gauge bosons', arising from the cancellation between different annihilation channels. However, the strength of the EWPT v c /T c > 1 cannot be obtained due to the small mass spitting ∆M, as found in Ref. [25]. In this work, the existence of the extra singlet scalar χ provides an additional contribution to the vacuum barrier at finite temperature, as explored in Sec. V. Therefore we get the chance to yield a strong enough EWPT as well as the correct DM relic abundance simultaneously 6 .
We find that with the increase of the mass of the singlet scalar, the mixing angle needs to be smaller to meet the strong first order EWPT condition, as shown in the top-left plot of Fig. 5. And also the criteria severely constrain the coupling λ L as shown in the top-right plots of Fig. 5 and bottom plots of Fig. 6.

DM relic density
In the top-left panel of Fig. 5 and top panels of Fig. 6, the mass difference of m 2 χ −m 2 H changes the DM-Higgs coupling a hHH efficiently (see Eq. 14 and Eq. 16), leading to the variation of relic abundance. The correct DM relic density constrains the mass of the singlet scalar χ to be smaller than 500 GeV.
From the top-right of Fig. 5 and bottom panels of Fig. 6, we find that the relic density highly constrains on the parameters λ L to be small (around zero). The reason is that the coupling a hHH ∼ 2vλ L cos 2 θ in this case, and annihilation cross section of DM pair can be suppressed by a smaller λ L in order to avoid the relic underabundance (as also discussed in Subsec. VI A 3).
In the large mass spitting region of m χ and m H , the term −µ soft sin 2 2θ becomes the dominant part of the coupling a hHH . So, as the a 2 hHH increasing, the annihilation cross section of DM becomes larger which causes the DM relic density deficiency, see the m χ − θ, m χ − λ 7 and m H − λ 7 plots in Fig. 5 and Fig. 6. Also, the DM relic density behaviors in m χ − λ 7 and m H − λ 7 planes are symmetrical with respect to λ 7 , this is because the annihilation cross section HH → hh is proportional to λ 2 7 .

DM direct detection
At the relative large DM mass region, the direct detection is less constraining. However, one still needs small λ L to suppress the Higgs portal coupling a hHH , LUX results rule out the most space of parameter λ L in the λ L −θ, m χ −λ L , λ L −λ 7 spaces. The exclusion in the bottom-right panels of Fig. 5 and Fig. 6 with respect to the varying DM mass tells that plenty parameter space is needed to probe in the underground experiments of the next generation.

A. Invisible Higgs decay
For the case of m h > 2m H,χ , the invisible decay of the SM like Higgs needs to be taken into account: ATLAS search [74] sets bounds on the invisible Higgs decay branching ratio Br(h → inv) < 28% at the 95% CL, which has been considered in our choice of DM scenarios with m χ,H ≤ m h /2, i.e., in the benchmark models BP1, BP3, and BP6.

B. Diphoton rate deviation
As in IDM, the charged state H ± leads to an additional contribution to the loop-induced h → γγ and γZ rates. [75,76]. The h → γγ rate can be obtained as [75][76][77][78][79][80] here, A SM ≈ −6.56+0.08i is the leading contribution coming from W bosons and top quarks for m h = 125 GeV. The contribution coming from H ± can enhance or suppress the h → γγ rate relative to the SM depending on the sign of λ 3 , which is characterized by the mass splitting of m H − m H ± as well as the mixing of the singlet and inert doublet scalars. For our benchmarks where strong first order EWPT and correct DM relic abundance can be accomplished, an tinny enhancement of h → γγ (around O(10 −3 )) is obtained due to λ 3 < 0. Other benchmark scenarios favor λ 3 > 0, which may lead to a large suppression of the rate, around O(10 −2 ). And the rate can be diluted by the invisible decay of SM-like Higgs when m H,χ < m h , as also noticed in Refs. [25,26,30].
C. Triple Higgs couplings v.s. EWPT As mentioned in the Section I, the EWBG stands out from many baryogenesis mechanisms. Based on the mechanism, the triple Higgs couplings might be modified which can be probed at the future hadron colliders, such as SPPC [4]. However, the strong first order phase transition not necessarily requires the enhancement or suppression of the triple Higgs coupling comparing with SM prediction, as studied in the singlet assisted case in Ref. [81]. Since there is no mixing between the SM-like Higgs and other neutral Higgses, we expect the model does not cause large modification of the triple Higgs couplings in the allowed benchmark scenarios. Generally, we can expect a tinny loop correction to the triple Higgs couplings (C hhh ) from new scalar sectors. The ratio of C hhh of our model with respect to the SM can be obtained as with C hhh = d 3 V e f f (T = 0)/dh 3 . The plots of ratio r 3h in Fig. 7 indicates that the deviation from the SM prediction is around O(10 −2 − 10 −1 ) in the benchmark models BP5 and BP-c2. When the triple Higgs vertices being modified by the magnitude of O(10 −1 ), which falls in the projection of the ILC [82] 7 , one can expect the benchmark modes being detected through the process of e + e − → Zhh. When the deviation is ∼ O(10 −2 ), one may expect the benchmark models to be probed by the SPPC [4,84] through the process of pp → hh.

D. Lepton collider and hadron collider prospects
At first, due to the mixing of χ and H in the model, one can expect the pair production process of H(χ)A at leptonic colliders [10,85], the corresponding cross section is where g H(χ)AZ = cos θ(sin θ) and is the space function of the two-body phase space. Unfortunately, for our two benchmark models (BP5 and BP-c2) that can address both the strongly first order EWPT and DM relic abundance, the future leptonic colliders, such as CEPC [86], and FCC-ee [87], lack the ability to search the H(χ)A pair productions, since such processes are kinematically unaccessible for the projected design of the centre-of-mass energy. We plot the cross sections predictions at √ s = 1000 GeV at ILC for the two benchmark models in Fig. 8. With the increasing of the mixing angle θ, one obtain the increasing (decreasing) of the cross section of e + e − → χ(H)A. As in the IDM, the added scalars may alleviate the Higgs natrualness problem since more bosonic fields contributions are added to the Veltman Condition [10,88]. While for the naturalness problem to be probed at leptonic colliders, one needs to expect a sizable modification of wave-function renormalization of the SM Higgs field after these extra Higgses and dark matter fields (H, χ, A) are integrated out [89][90][91]. We leave the detailed study on interplay between naturalness probe and EWPT to the future study.
At the Hadronic Collider, the dominant production channels of the extra Higgs pairs in the model can be 8 and followed by the decay channels being We explore the possible prospects of the future search of the two benchmark models: BP5 and BP-c2. We use MadGraph5 [92] to estimate cross sections of different channels, the cross section of the H ± A, H ± H and H ± χ channels of Eq. 54 are estimated to be 2.9 fb, 2.1 fb and 0.6 fb at 14 TeV, the cross section of the HA, χA and H + H − channels of Eq. 55 are estimated to be 1.6 fb, 0.5 fb and 0.7 fb. As in the usual IDM models, the multi-lepton plus missing energy signatures [18,93,94] can be expected in the final states at the detector level. While, due to the smallness of the cross sections, these channels might be unreachable by LHC run-2 at 14 TeV or even SPPC. More detailed studies are left to the future projects. As for the BP-c2, the scalar mass relation of m H ± = m A > m H indicate that the followed decay of Eq. 57 is replaced by H ± → HW * , the macroscopic decay lengths of which can be cτ ∼ 1mm [60] for the mass splitting between H ± (A) and H being 2 GeV. And the bottom-right panels of the Fig. 5 and Fig. 6 illustrate that m H can be very close to and even degenerate with the m H ± ,A for θ > 0.5 and λ 7 > 0.07, which means one can expect cτ > 1mm and therefore the displaced decay of H ± → HW * and A → HZ * can be searched through a mono-jet with a soft displaced vertex [95][96][97]. This signature can make the benchmark model different from the IDM case [60]. When the λ L ≤ −0.2, one have m H ∼ m A,H ± for dark matter relic abundance undersaturated case of BP5. Then, one can also expect the displaced vertex signature.

VIII. CONCLUSION AND DISCUSSION
In this work, we study the dark matter and electroweak phase transition in the framework of mixed scalar dark matter model. In the model, the dark matter candidate is the lightest Z 2 -odd mixed particle coming from the mixing of the singlet scalar and the CP-even neutral component of the inert doublet. We investigate the scenarios where the dark matter particle is mostly the CP-even neutral component of the inert doublet (the benchmark scenario with H being the DM candidate) or the singlet (with χ being the DM candidate), which is determined by the mixing of the two. After imposing theoretical constraints of unitarity, stability and electroweak precision test(in particular the T parameter), we explore several benchmark models classified by the mass of the dark matter particle. The dark matter annihilation process can be affected a lot by the mixing and co-annihilation effects. Two benchmark models in which the strongly first order electroweak phase transition can be accomplished have been investigated.
The effects of quartic scalar couplings (λ L , λ 7 ), DM mass, mass splittings among different extra scalars, and mixing angle on the phase transition as well as dark matter phenomenology have been studied. Certain mass differences between the mixed scalars χ, H and other scalars in the inert doublet (H ± and A) are required to successfully realize the strongly first order electroweak phase transition. The DM direct detection limits can be evaded through tuning the mixing parameter. Lepton and hadron collider prospects have also been addressed in the model.
The unbroken Z 2 symmetry in the model precludes spontaneous and explicit CP violation arising from the Higgs potential, because of which the model fails to accommodate successful electroweak baryogenesis to explain the baryon asymmetry of the universe. One remedy method could be the introduction of effective high dimensional operators [98,99], which is beyond the scope of this study.