A simultaneous study of dark matter and phase transition: two-scalar scenario

The simplest extension of the Standard Model by only one real singlet scalar can explain the observed dark matter relic density while giving simultaneously a strongly first-order electroweak phase transition in the early universe. However, after imposing the invisible Higgs decay constraint from the LHC, the parameter space of the single scalar model shrinks to regions with only a few percentage of the DM relic abundance and when adding the direct detection bound, e.g. from XENON100, it gets excluded completely. In this paper, we extend the Standard Model with two real gauge singlet scalars, here s and s′, and show that the electroweak symmetry breaking may occur via different channels. Despite very restrictive first-order phase transition conditions for the two-scalar model in comparison to the single scalar model, there is a viable space of parameters in different phase transition channels that simultaneously explains a fraction or the whole dark matter relic density, a strongly first-order electroweak phase transition and still evading the direct detection bounds from the latest LUX/XENON experiments while respecting the invisible Higgs decay width constraint from the LHC.


Introduction
There have been numerous attempts in different mainstreams from modifying the theory of gravity to extending the Standard Model (SM) of elementary particles to accommodate the problem of the missing mass or the dark matter (DM). The most successful example of the latter is the ΛCDM model which incorporates the existence of a cosmological constant (responsible for the accelerating expansion of the universe) and a cold dark matter (CDM) as new species of particle(s) living in a dark sector. Within the ΛCDM model, the weakly interacting massive particle (WIMP) has been specially a successful DM paradigm. A WIMP candidate of dark matter can be embedded in various extensions of the standard model with Z 2 or larger symmetry groups in the hidden (dark) sector. In this paper, we investigate the scalar extension of the SM with the Z 2 discrete symmetry group needed for stabilizing the dark matter candidate in the so-called freeze-out mechanism.

JHEP12(2019)077
The first and the simplest of such models is the extension of the SM with only one real single scalar (only one degree of freedom) which has been studied vastly in the literature, see e.g.  in which various phenomenological aspects such as the dark matter relic density extracted from WMAP and Planck [24][25][26][27], the Higgs invisible decay width from the LHC experiments [28,29], the upper bound on the dark matter elastic scattering cross section off nuclei by XENON100, XENON1T, LUX [30][31][32], gamma rays from annihilation of dark matter interpreted by Fermi-LAT data [33,34], and the theoretical aspects such as perturbativity, vacuum stability, electroweak phase transition, gravitational waves have been investigated. Despite having a small space of parameter, the model is remarkably successful in addressing a subset of the aforementioned constraints. The real scalar field in the single scalar dark matter model plays the role of both the DM candidate and the DM-SM mediator through the Higgs portal.
The status of this model has been reported by the GAMBIT Collaboration in [35]. According to the GAMBIT, taking into account all the direct and indirect constraints (without imposing the conditions for the first-order phase transition) the model remains alive whether the singlet scalar stands for only a fraction or the whole dark matter relic abundance. The viable parameter space with couplings of order unity lies in the DM mass between the Higgs mass and 300 GeV or above 1 TeV where for the latter the scalar field can constitute all the DM content. On the other hand, the real singlet scalar model is also capable of giving a first-order electroweak phase transition (EWPT) from the symmetric phase to the broken phase of the SU(2) electroweak gauge symmetry group. Some papers that have studied also the phase transition in the single scalar model are [36][37][38][39][40][41][42][43][44][45][46]. It has been pointed out that the dark matter constraints are strongly in conflict with the firstorder phase transition conditions (see e.g. [47]), while we have shown in [48] that in fact the observed dark matter relic density and the first-order electroweak phase transition are consistent, but the parameter space significantly gets reduced only after imposing the invisible Higgs decay constraint and gets completely excluded after considering the bounds from the direct detection experiments.
It was shown in [49] that by adding another real singlet scalar to the single scalar model, the problem of the restrictive direct detection constraints gets resolved non-trivially. The presence of a term ss H † H in the Lagrangian used in [49], with s and s being the two real singlet scalars, is a key term in making the two-scalar model significantly different from simply summing up two single scalar models. The idea of two-scalar extension of the SM has been explored also in [23] to address the DM and in [50,51] for the Higgs inflation and the electroweak phase transition. Models with a complex scalar field or a composite Higgs have been studied in [52][53][54][55][56][57][58][59][60] and models with multi-scalar extension of the SM are found in [51,61]. Among the scalar extensions of the SM, the question of the electroweak phase transition has been also studied extensively in two Higgs doublet models (2HDM) .
In this paper we investigate in detail the question of the strongly first-order electroweak phase transition and the problem of dark matter simultaneously in the two-scalar model. Having two real scalars, in addition to the Higgs doublet scalar field H, the fields configuration space becomes three dimensional which in turn makes the structure of the phase transitions richer. Let us assume the vacuum expectation values (VEV) of the Higgs and JHEP12(2019)077 the two extra scalars i.e. the VEVs of (H, s, s ), by (v sym , w 1 , w 1 ) in the symmetric phase and (v brk , w 2 , w 2 ) in the broken phase. In high temperatures that the electroweak symmetry group SU(2) × U(1) is not broken, the VEV of the Higgs is vanishing, therefore throughout the paper we set v sym = 0. To stabilize the dark matter candidate, here chosen to be the scalar s, by a discrete symmetry group Z 2 , we need to set the VEV of the dark matter to zero after the phase transition, i.e. w 2 = 0. Therefore the general form of the transition from symmetric to broken phase is (0, w 1 , w 1 ) → (v brk , 0, w 2 ). We analyze different scenarios depending on possible values of w 1 , w 1 , v brk and w 2 to give a strongly first-order phase transition. We then combine the analytic conditions of the EWPT with the constraints from the direct and indirect dark matter searches. Despite the strong bounds from the first-order EWPT and the direct detection constraints which excludes completely the single scalar model, the two-scalar model evades remarkably all these constraints at the same time and predicts viable dark matter models.
The paper is organized as follows. In section 2 we show analytically that there are different channels of the EWPT and obtain the necessary conditions for the EWPT to be of the first-order type. The section is divided into two subsection with two two-scalar models; one without the s-s cross-coupling terms and the other including these terms. Then in section 3 we elaborate the DM relic density and direct detection constraints. In section 4 we numerically search for the viable space of parameters combining the strongly first-order EWPT conditions, the observed dark matter relic density, the direct detection constraints and the limit of the invisible Higgs decay width. We also compare the results with the single scalar model exposed to the first-order EWPT and the DM direct and indirect bounds. We conclude and summarize in section 5. In appendices A and B we bring the details of finding the minima of the scalars configuration (H, s, s ) and the deepest minimum condition respectively.

First-order phase transition
The strongly first-order phase transition in the early universe is one of the three Sakharov conditions [87] for the electroweak baryogenesis. In high temperature of the early universe, the electroweak symmetry group is unbroken and rest in its symmetric phase, SU(2) ×U(1), with the Higgs VEV vanishing, but as the universe expanded, i.e. at lower temperatures, the vacuum acquires a non-vanishing VEV and the symmetry is broken into U(1) electromagnetic gauge group. In the SM framework, a strong first-order phase transition gives the Higgs mass an upper limit, m H 70 GeV [88,89], which is in conflict with the measured Higgs mass at the LHC being 125 GeV. This motivates the extension of the SM which among numerous possible extensions the addition of a real singlet scalar is the simplest. However as it has been shown in [48], the viable space of parameters survived from the DM and the EWPT constraints, gets excluded mostly by taking into account the invisible Higgs decay constraint. Here we investigate the idea of extending the SM by two real scalars and examine the model against the simultaneous consideration of the DM and the EWPT along side the imposition of the direct and indirect probes. In the two-scalar model we restrict ourselves to only terms with dimensionless couplings, therefore terms such as JHEP12(2019)077 s 3 or ss 2 are not present. Moreover, we analyze the model in two parts, once with the s-s cross-coupling terms for the scalars s and s , i.e. s 2 s 2 and ss 3 and s 3 s , and once without these s-s cross-coupling terms.
It can be seen from eq. (A.1) that at very high temperature, T → ∞, the only extremum of the thermal effective potential is the point (v = 0, w = 0, w = 0) in VEV space. However, with the expansion of the universe as the temperature decreases, the non-zero local minima for the scalars come into existence. So in principle as the universe cools down from very high to very low temperature, any of the scalar fields h, s and s may undergo more than once a transition from a symmetric phase to a broken phase. In this paper by electroweak phase transition we mean the transition from a vanishing VEV into a non-zero VEV for the Higgs scalar field. We are not interested here in considering the scenarios of the symmetry breaking in the dark sector. Therefore, in a symbolic transition from (v sym , w 1 , w 1 ) → (v brk , w 2 , w 2 ) in the VEV space, the parameters w 1 , w 1 are the VEVs of the scalars s, s at temperature T c + δT 1 and w 2 , w 2 are the VEV's of the scalars at T c − δT 2 for some arbitrary δT 1 and δT 2 and for T c being the critical temperature at which the phase transition in the Higgs sector triggers. The phase transition may continue until T c − δT 2 approaches zero or it may end before the zero temperature. On the other hand, δT 1 can be arbitrarily small so that it is enough for the VEVs, w 1 and w 1 , to exist before or very close to the critical temperature. We will follow this strategy throughout the paper.

Model without s-s cross-coupling terms
The potential of the model possessing two extra real scalars beyond the SM along side the Higgs doublet, and without the s-s cross-coupling terms reads, where H † = 1 √ 2 (0 v +h) denotes the Higgs doublet scalar field after the symmetry breaking, and s, s are the two real singlet scalar fields. The dominant one-loop thermal effective potential is given by, 1

JHEP12(2019)077
The thermal effective potential is obtained by summing up eq. (2.1) and eq. (2.2), 2 that explicitly is given by eq. (A.1) ignoring the cross-coupling terms. The potential is invariant under Z 2 transformation if applied for both scalars at the same time, s → −s, s → −s . It means only through both scalars s and s the Z 2 symmetry is reserved, therefore the lighter scalar could be assumed as the dark matter candidate. In this paper, we take the scalar s to be the dark matter particle. The effect of the thermal correction is only in the mass term of the tree-level potential. So in the total effective potential instead of the coefficients µ 2 h , µ 2 s and µ 2 s we deal with T -dependent masses µ 2 h (T ), µ 2 s (T ) and µ 2 s (T ) which are defined in (A.2).
The VEVs of the scalar fields (h, s, s ) can take different values before and after the EWPT. What is important to have in mind, is that after the EWPT, the VEV of the Higgs particle should be non-zero and the VEV of the lighter scalar field which is the DM candidate must be vanishing. Therefore, the most general structure of the VEVs would be (v = 0, 0, w = 0). It is shown in appendix A that after the EWPT if we choose one of the scalar's VEV to be zero, the other scalar must take a vanishing VEV as well. So the only possibility for the VEVs after the EWPT is (v = 0, 0, 0).

Phase transition scenarios. As seen in appendix
λ h , 0, 0) in which µ 2 h (T ), µ 2 s (T ) and µ 2 s (T ) are defined in eq. (A.2). Note these are only some selected solutions and in general there are more complicated expressions for w 2 and w 2 . We analyze all these four possible transitions one by one to figure out which can be of first order type.

Phase transition
In this scenario only the Higgs particle undergoes a non-zero VEV while the other two scalars keep the Z 2 discrete symmetry in all low and high temperatures. In order for (0, 0, 0) to be a local minimum, eq. (A.6) must be satisfied. This would leave us with a set of constraints on µ 2 (T )'s,

JHEP12(2019)077
and a similar set of conditions must hold for (v 2 = The conditions on µ 2 h (T ) in eqs. (2.5) and (2.6a) are clearly inconsistent, which means that the two minima cannot coexist. Therefore, the first-order phase transition from (0, 0, 0) to (v, 0, 0) is not possible.

Phase transition
In appendix A, we see that (v = 0, w = 0, w = 0) with w 2 = µ 2 s (T ) λs is an extremum of the potential. Here we examine the transition from (v = 0, w 2 = µ 2 s (T ) λ h , w = 0, w = 0). In other words, the mediator scalar, s , gets always vanishing VEV before and after the phase transition. Then at high temperature, the Higgs has a zero VEV and the DM particle, s, has a non-zero VEV. This situation is closely related to the real single scalar dark matter model with the difference that here there is an additional real singlet scalar with a vanishing VEV. The minimum conditions for the point (v = 0, w = 0, w = 0) is given in eq. (2.6) and those for the VEV point (v = 0, w = 0, w = 0) can be extracted from eq. (A.6) in appendix A, The critical temperature below which the universe starts a transition from vanishing Higgs VEV to non-zero VEV, is given by the following expressions, For T T c it is necessary that both (0, w, 0) and (v, 0, 0) be local minima of the potential. Furthermore, the point (v, 0, 0) in the VEV space must be as well a global minimum for temperature below T c . It can be shown that eq. (2.7) holds for all values of the temperature in 0 < T < T c , if it holds at T = 0 and T = T c , which leads to,

JHEP12(2019)077
Then eq. (2.6) holds for all T T c if it holds only at T = 0 and T = T c , The last condition that must be considered in this scenario is that the minimum (v, 0, 0) should be the global one for the temperatures below the critical temperature. That is, from eq. (B.2) for all T < T c , Equivalently, one can translate this constraint in where use has been made of eq. (2.7a) and the following equality at T = T c from the definition of the critical temperature, This scenario is very similar to the last one, with the difference that here the DM candidate scalar, s, always takes zero VEV but the heavier scalar, s , goes from non-zero VEV before EWPT at high temperature to zero VEV at temperatures lower than the critical temperature. The local minimum conditions for the VEVs at low temperature after the EWPT, i.e. for (v, 0, 0) are those given in eq. (2.6). The conditions for above the critical temperature are given by eq. (2.7), but with an interchange in the scalar fields, i.e. s ↔ s . The critical temperature similarly is obtained, (2.14) For the VEV point (v, 0, 0) to be the deepest minimum after the phase transition we have, In this scenario both scalars s and s have non-zero VEVs above the critical temperature and both get zero VEV after the phase transition takes place. As is discussed in appendix B, the critical temperature can be obtained from eq. (B.1), The local minima conditions for the VEVs (v = 0, w = 0, w = 0) before the EWPT are now more involved, where must be satisfied at least for all T T c . These conditions at T = 0 yields, In order to have a first order transition from symmetric phase to broken symmetry phase of the Higgs vacuum, the VEV set (v, 0, 0) must be a global minimum. This condition is obtained via eqs. (B.2) and is given by, JHEP12(2019)077

Model including s-s cross-coupling terms
In the previous subsection, we ignored the s-s cross-coupling terms, i.e. the interaction terms consisting only the singlet scalars, s and s . If we include also these terms, the total tree-potential would be the sum of the potential in eq. (2.1) and the s-s cross-coupling terms, Note that we have considered only the cross-coupling terms with dimensionless couplings. The one-loop thermal potential in this scenario has the same form as in eq. (2.2), however the coefficients c h , c s and c s are now different from eqs. (2.3) as now there are more one-loop Feynman diagrams for thermal mass corrections, As seen in eqs. (2.23b) and (2.23c), the coupling λ ss appears in the thermal corrections. The reason is that the one-loop thermal mass correction for the scalar s (s ), in addition to the Higgs field, includes as well the scalar s (s) in the loop. However, the couplings λ ss and λ ss , although playing a role in first-order phase transition conditions, but they do not enter directly in the mass thermal corrections.
Phase transition scenarios. Finding a complete set of the extrema (v, w, w ) (with v, w and w being the VEV's of h, s and s respectively) from eqs. (A.4), for the general case of totally non-vanishing λ ss , λ ss and λ ss , is possible but the solutions are very lengthy. Therefore, in the following subsections we consider only the solutions which with no loss of generality are simpler and can also be compared with the phase transitions in model without s-s cross-coupling terms.

Phase transition
The phase transition here is from . This extremum solution of the potential is possible for at least two different choices of the s-s cross-couplings, i.e. for λ ss = 0 , λ ss = 0 , λ ss = 0 and for λ ss = 0 , λ ss = 0 , λ ss = 0 . Here we derive the first-order phase transition conditions for the first set of the s-s cross-couplings above which turns out to be the same as the other set of coupling. The local minimum conditions using eqs. (A.5) and (A.6) are, where similar to the lines in section 2.1, it is enough that eqs. (2.24) satisfy for T = 0 and Similarly the local minimum conditions for the VEV set (v, 0, 0) with v 2 = µ 2 h λ h , must be driven from eqs. (A.5) and (A.6). It turns out that these conditions for the model with the s-s cross-coupling terms is the same as those for the model without the s-s cross-coupling terms in subsection 2.1, i.e, in eqs. (2.6).
The critical temperature is obtained from the degeneracy condition in eq. (B.1) and is given by, After the phase transition, the minimum in the broken phase needs to be a global minimum which is translated into, In this scenario the phase transition is from (v = 0, w = µ 2 s (T ) λ h , w = 0, w = 0). It means that the DM candidate takes non-zero VEV before the phase transition and its VEV flips to zero after the phase transition to retain the Z 2 symmetry. Again there are two sets of the s-s cross-couplings for which the VEV set before the phase transition is an extremum solution to the potential in eq. (2.22): λ ss = 0 , λ ss = 0 , λ ss = 0 and for λ ss = 0 , λ ss = 0 , λ ss = 0 . For both sets of the couplings, the local minimum conditions where T c , the critical temperature, for this scenario is given by, In order for the minimum (v, 0, 0) to be deeper than (0, w, 0), the following condition must be held, where the degeneracy condition on the potential at the critical temperature has been used. Again, the local minimum conditions for (v, 0, 0) is given by eqs. (2.6).

Phase transition
This type of the phase transition from (v = 0, w = 0, w = 0) to (v = 0, w = 0, w = 0), as seen in the appendix A.2, is possible for a choice of the s-s cross-couplings being λ ss = 0 , λ ss = 0 , λ ss = 0 with w and w given by, (2.32) The λ h is given by eq. (2.6). The critical temperature reads, (2.35) Now the minimum (v, 0, 0) after the electroweak symmetry breaking must be deepest minimum and therefore eq. (B.2) must be satisfied, and w, w are given by eq. (2.32). The above inequality leads to the following condition, where a is the parameter defined in eq. (2.35). One of the Sakharov conditions for the baryogenesis is the washout criterion which guarantees an appropriate sphaleron rate to have a strongly first-order phase transition. This condition is translated into an inequality as v c /T c > 1. In all the numerical computation and for each phase transition scenario, we impose also the washout criterion.

Dark matter constraints
In the last section we elaborately studied the simplest possible phase transitions that may occur for going from the symmetric phase to the broken phase of the Higgs vacuum. In this section we discuss the dark matter constraints for the two-scalar model and in the next section we combine these constraints with those of the strongly first-order phase transition in the last section and represent the final results.
From the last section, having observed that the only possibility for the VEVs of the scalars after the EWPT is the VEV structure (v, 0, 0), let us study the mixing among the scalars after the electroweak symmetry breaking. The mass matrix after the Higgs particle gets its non-zero VEV is not diagonal, although both scalars s and s are at their zero VEVs. Diagonalizing the mass matrix can be done by a rotation in the (s, s ) field configuration, where φ and φ are new scalars with diagonalized mass eigenvalues, and θ is the mixing angle and is given by, in which, are the masses of the scalars s and s which are related to the VEV of the Higgs scalar, v. The extremum condition of the potential in eq. (2.1) at (v, 0, 0) gives µ 2 h = v 2 λ h . The entries of the diagonalized mass matrix after the rotation from the field configuration (s, s ) into (φ, φ ) read, Therefore, in two-scalar model we have two singlet scalar WIMPs, φ and φ where we assume the DM candidate is the φ field being the stable WIMP, and the heavier WIMP, φ , is unstable and can decay to the SM particles plus the light WIMP through an intermediate Higgs, i.e. φ → φ + SM. Let us define the mass difference of the two scalars as δ = m φ − m φ . For the mass splitting of O(GeV) and beyond, the life time of the heavy WIMP will be much smaller than the age of the universe and therefore cannot have effective contribution to the present DM relic density [49]. One of the important constraints on the DM models comes from the observed DM relic density given by WMAP and Planck experiments. In this work we assume the present DM relic density is produced thermally via the so-called freeze-out mechanism taken place at some specific freeze-out temperature, T f , in the early universe [92]. Let us discuss briefly how this mechanism works. At high temperatures we believe that the WIMPs and the SM particles are in thermal equilibrium in an expanding universe. It means that the WIMPs annihilation into the SM particles and, the WIMPs productions take place at the same rate. As the universe expands, there is an epoch after which the universe expansion rate surpasses the WIMPs annihilation rate such that it becomes infrequent for the WIMPs to meet each other for the annihilation to happen. At this point in time, the temperature is low enough so that the SM particles possess insufficient kinetic energy to produce WIMPs. This is the epoch in which DM particles decouple from the SM particles, thenceforth the DM density remains constant asymptotically in the comoving volume.
The freeze-out temperature and hence the DM relic density depend strongly on various types of the WIMP interactions with the SM particles. The dominant contributions are due to DM annihilation cross section and the less important contributions come from coannihilation processes. In the latter case, the DM candidate along with the heavier WIMP annihilate into the SM particles. As an example, the relevant (co)annihilation The change in the number densities of the two WIMPs in terms of the temperature are controlled by two coupled Boltzmann equations. Instead of solving the two coupled equations which is not a simple task, one can solve a single Boltzmann equation with an effective DM cross section incorporating both annihilation and coannihilation cross sections [93,94]. If we take the total number density as n = n φ + n φ , the effective Boltzmann equation reads, where the effective cross section is defined as Here, σ φφ , σ φ φ and σ φφ denote respectively, the DM annihilation to the SM particles, the heavier WIMP annihilation to the SM particles and the coannihilation to the SM particles. The effective number of degrees of freedom is g eff = 1+(1+ δ/m φ ) 3/2 e −δ/T , and the Hubble constant in the Boltzmann equation is denoted by H. The thermal averaging of the effective cross section multiplied by the relative DM velocity at temperature T is defined as σ eff v .

JHEP12(2019)077
φ φ h qq Figure 2. The DM-quark direct detection scattering cross section is shown at the leading order in perturbation theory.
The second important constraint that one should impose on the parameter space is the stringent exclusion limits from the DM-nucleon elastic scattering cross sections. These limits are provided by dark matter direction detection experiments among them we exploit here the latest updates of LUX [95] and XENON1T [31]. The underlying interaction which leads to DM-nucleon elastic scattering is given by an effective Lagrangian describing the DM-quark interaction, where the effective coupling c q is obtained in terms of the relevant couplings in the Lagrangian, the mixing angle, the quark mass, and the Higgs mass as follows, There is a standard method by which one can promote the quark-level effective Lagrangian to hadron-level interaction at zero-momentum transfer [96,97]. This can be achieved if we replace the quark current by a nucleon current up to a low energy effective factor c N as, For the DM-proton scattering cross section we use these scalar couplings, f p u = 0.0153, [98]. The final formula for the DM-proton SI elastic scattering cross section is, where µ p is the reduced mass of the DM and the proton.

Numerical results
In this section we impose simultaneously all the dark matter constraints from section 3 and the strongly first-order phase transition from section 2, in our numerical computations. In order to compute numerically the DM relic density we apply the package MicrOMEGAs [99] which requires the implementation of our model into the program LanHEP [100]. The WMAP [24] and Planck [25] measurements of the cosmic microwave background (CMB) strongly constrain the mean density of cold dark matter (CDM). The recent Planck result yields Ω CDM h 2 = 0.12 ± 0.001 [26]. In our analysis we assume that the scalar DM candidate fully or partially saturates the observed relic density such that The DM phenomenology of the present model is fully studied in [49]. However, we recap some main results therein. We recall that in the simplest extension to the SM, with a singlet scalar DM candidate, except the resonance region the rest of the parameter space   is excluded by the recent direct detection (DD) bounds. One of the characteristics that the two-scalar DM model inherits and is absent in the single scalar model is manifested by the regions in the parameter space which evade the current DD upper limits.
In the single scalar model, the DM-nucleon scattering cross section and the annihilation cross section are both proportional to a single coupling constant. Regions in the parameter space with large enough coupling constant giving rise to the correct relic abundance, have large DD scattering cross section which are excluded by the present DD experiments.
In our extended scalar model, when we have two particles in the final state, there is a DM annihilation process with a heavy WIMP mediated in t-or u-channel, see the top-left diagram in figure 1. The presence of this process is critical in the analysis, because this process enters a contribution with a coupling other than that in the DD cross section. Therefore it becomes plausible to find viable regions in the parameter space with small coupling for dark matter elastic scattering cross section and hence small DD cross section, and at the same time large enough dark matter annihilation coupling to induce the correct DM relic abundance.
We consider two models in the following analysis. In the first case the s-s crosscoupling terms are absent in the Lagrangian, i.e. λ ss = λ ss = λ ss = 0 as discussed in subsection 2.1 to obtain the first-order EWPT in the model without s-s cross-coupling terms. The independent free parameters are λ hs , λ hs , λ s , λ s , m s , m s and the mixing angle θ. The coupling constant λ hss is given in terms of the mixing angle and WIMP masses in eq. (3.2). In the second scenario studied in subsection 2.2, the s-s cross-coupling terms are included and the dimension of the parameter space is increased. The set of the free parameters in this case is λ hs , λ hs , λ s , λ s , λ ss ,λ ss , λ ss , m s , m s , and the mixing angle θ. In all phase transition scenarios discussed in subsections 2.1 and 2.2 we perform a full scan with 1.6 × 10 8 samplings over the parameter space in the following parameter intervals: 10 GeV < m φ < 5 TeV, m φ = m φ + δ, 1 GeV < δ < 1 TeV, 0 < λ hs , λ hs < 1, 0 < λ s , λ s < 2, 0 < sin θ < 1, and when relevant, 0 < λ ss , λ ss , λ ss < 2.
In the model without the s-s cross-coupling terms there are four scenarios for the electroweak phase transition. The first scenario 2.1.1, does not give rise to a first-order phase transition because of an internal inconsistency in the first-order conditions. For zero and one. Finally the lower-right plot shows that the critical temperature is of order 10 GeV. Note that it has been assumed that the phase transition takes place above the DM freeze-out temperature. In table 1 a list of benchmarks have been represented. The maximum percentage of DM relic density that can be accounted by the scalar s, is ∼ 25% for the DM with mass of ∼ 210 GeV and T c ∼ 9 GeV. It should be noted also that from the ratio v c /T c in table 1, it is obvious that the phase transition is very strong.
The model with s-s cross-coupling terms consists of three scenarios that for each one we have found a viable space of parameters. In scenario 2.2.1 i.e. for a phase transition from (0, 0, w ) to (v, 0, 0) as it is seen in figure 4  The benchmarks represented in all the tables give at least a fraction of the DM relic density and at the same time are consistent with a strong first-order phase transition while evading the restrictive direct detection bounds e.g. from LUX2017/XENON1T and survive also from the invisible Higgs decay constraint. Despite the very restrictive constraints from the first-order phase transition and the direct detection bounds, we observe that the twoscalar model predicts models of dark matter that remarkably evades all the constraints JHEP12(2019)077 simultaneously. A quantity which is measured by the CMS and the ATLAS is called signal rate. Since in our model the SM Higgs is not mixed with the singlet scalars, this measured quantity puts no constraint on the model parameter space. For the same reason, there is no constraint due to exclusion limits from non SM-like Higgs boson search. Concerning the electroweak precision measurements, in our model the singlet scalars do not mix with the SM gauge bosons and therefore electroweak precision observables are not affected by the singlet scalars.

Conclusion
In this paper an extension to the SM with two real singlet scalar (dubbed as two-scalar scenario) denoted here by s and s has been investigated to examine whether the model is capable to accommodate simultaneously several constraints from thermal processes such as the relic density of dark matter and the strongly first-order electroweak phase transition in the early universe to constraints from the direct detection experiments and the invisible Higgs decays bound at the LHC. It is known from the literature that the single scalar extension of the SM fails to explain simultaneously the following constraints: the observed relic density, the first-order EWPT and the invisible Higgs decay width limit. However in [48] it was shown that two sets of conditions from the DM relic density and the firstorder EWPT are not in fact in conflict in the single scalar model but the space parameter shrinks to regions with a few percentage of the DM relic density when the invisible Higgs decay constraint is imposed. We have shown in two-scalar model that there can be different phase transition channels from the symmetric phase to broken phase of the Higgs vacuum. Despite very restrictive constraints from the first-order EWPT conditions which is more restrictive than the EWPT condition in the single scalar model, some of the channels in the two-scalar model can explain a fraction or the whole observed DM relic density, and at the same time the strongly first-order EWPT, the direct detection bounds from LUX/XENON1T and the invisible Higgs decay constraint. We have also represented the benchmarks for each phase transition scenario showing the viable range of the DM mass and all the corresponding parameters.

Acknowledgments
Arak University is acknowledged for financial support under contract no.98/664.

A Minima in 3-dimensional VEV space
The most general three-level potential we have considered in this paper consists of two extra singlet scalars in addition to the Higgs field. Taking into account the thermal contributions (the one-loop contribution is negligible) we have, Let us assume that the extremum of this potential is located at (v, w, w ), then the first derivatives at this point is vanishing, Eq. (A.3) leads to the following set of equations, As seen from eq. (A.4a), both v = 0 and v = 0 are extrema of the potential. Let us also define the second derivatives of the potential at the extremum point (v, w, w ) as the following, V ss ≡ ∂ 2 V ∂s 2 (v,w,w ) = −µ 2 s (T ) + 3λ s w 2 + λ hs v 2 + λ ss w 2 + λ ss w 2 , (A.5b) V s s ≡ ∂ 2 V ∂s 2 (v,w,w ) = −µ 2 s (T ) + 3λ s w 2 + λ hs v 2 + λ ss w 2 + 2λ ss ww , (A.5c) V hs ≡ ∂ 2 V ∂h∂s (v,w,w ) = 2λ hs vw + λ hss vw , (A.5d) V hs ≡ ∂ 2 V ∂h∂s (v,w,w ) = 2λ hs vw + λ hss vw , (A.5e) V ss ≡ ∂ 2 V ∂s∂s (v,w,w ) = 1 2 λ hss v 2 + 2λ ss ww + λ ss w 2 + λ ss w 2 . (A.5f) The conditions for the point (v, w, w ) to be a local minimum are, A.1 Model without s-s cross-coupling terms The first case we have considered in this paper is when there is no s-s cross-coupling terms in the potential in eq. (A.1), i.e., the case λ ss = λ ss = λ ss = 0. Eqs. (A.4) then is simplified as, v −µ 2 h (T ) + λ h v 2 + λ hs w 2 + λ hs w 2 + λ hss ww = 0 , (A.7a) −µ 2 s (T )w + λ s w 3 + λ hs wv 2 + 1 2 λ hss w v 2 + λ ss ww 2 = 0 , (A.7b) −µ 2 s (T )w + λ s w 3 + λ hs w v 2 + 1 2 λ hss wv 2 + λ ss w w 2 = 0 . In v = 0 class if w = 0 then we must have also w = 0, and vice versa. Therefore the solutions in this class are only in the following forms, After the electroweak symmetry breaking the Higgs vacuum expectation value is non-zero, but if the scalar s wants to be the DM candidate it must take zero VEV after the EWPT (or to be more accurate after the DM freeze-out). Therefore, the only vacuum structure of the two-scalar model after the EWPT is (v 2 = µ 2 h (T ) λ h , 0, 0). The extremum (v, w, w ) must be also local minimum, at least in some temperature intervals, as has been discussed throughout the paper. The second derivatives at the extremum point (v, w, w ) for the model without the s-s cross-coupling terms read, V hh = −µ 2 h (T ) + 3λ h v 2 + λ hs w 2 + λ hs w 2 + λ hss ww ,

B Critical temperature and deepest minimum condition
Let us represent the VEVs of the scalars before the phase transition i.e. in the symmetric phase, as (v sym , w 1 , w 1 ) and after the phase transition i.e. in the broken phase as (v brk , w 2 , w 2 ). Note by the symmetric and broken phase we mean only in the electroweak symmetry group SU(2) and we do not in general consider the symmetry status of other scalar field in the theory. The critical temperature is defined as the temperature at which the symmetric and broken minima of the thermal effective potential become degenerate, therefore, V eff (sym, w 1 , w 1 ; T c ) = V eff (v brk , w 2 , w 2 ; T c ) . (B.1) In order for the local minimum (v sym , w 2 , w 2 ) to be also the global one, it must be deeper than the local minimum (v sym , w 1 , w 1 ), i.e., ∆V eff (T ) ≡ V eff (0, w 1 , w 1 ; T ) − V eff (v brk , w 2 , w 2 ; T ) > 0 , (B.2) which must hold for all T < T c . To require the condition (B.2) to satisfy for T < T c it is enough that the T -derivative of ∆V eff (T ) be negative at T c .

JHEP12(2019)077
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.