Two-step strongly first-order electroweak phase transition modified FIMP dark matter, gravitational wave signals, and the neutrino mass

We study a dynamical freeze-in production of the dark matter considering the electroweak phase transition history of the Universe. The kinematical thresholds of the decay and scattering processes for the dark matter production can be altered by the temperature dependent thermal masses of particles, which might lead to enhancement or reduction of the dark matter relic abundance. The second-stage strongly first-order electroweak phase transition (SFOEWPT) triggered by the hidden scalars can be probed at colliders and the gravitational wave detectors. The two-step SFOEWPT modified late decay FIMP dark matter is accomplished with a Dirac neutrino mass explanation in the scotogenic model.


Introduction
The dark matter (DM), as an essential ingredient of the standard cosmological model, is supported by cosmological observations [1] and can be explained as particles beyond the standard model (SM). The non-gravitational property of the DM is not known to us [2]. The explanations of the DM motivates the Weakly Interacting Massive Particles (WIMPs), where the DM weakly interacts with the SM sector but strongly enough to be probed at various direct detection experiments [3][4][5], indirect detection experiments [6][7][8][9][10] and colliders, such as LHC, and future CEPC, ILC, and FCC-ee. However, the lack of the signals of WIMPs at these detectors motivates people to go beyond the WIMPs paradigm [11]. In the WIMPs paradigm, the DM particles live in the thermal bath at high temperature through the interaction with the SM particles, and the DM relic abundance is produced by the freeze-out mechanism after the DM particles annihilation rate falls below the Universe expansion rate at temperature around T f o ∼ m DM /x f o with x f o ∼ 26 [12]. One alternative approach can be the Feebly Interacting Massive Particle (FIMP), where the DM particles never enter into thermal equilibrium with the SM plasma due to the super-weak or feebly interaction rate between the DM particles and the SM particles [13,14,16]. In this paradigm, a basic assumption is that the initial DM abundance (at reheating epoch) is negligible 1 , and the present DM abundance is produced by the so-called freeze-in mechanism, through the SM particles (in the thermal bath) decaying or annihilating to the DM particles at high temperature around T f i ∼ m DM /x f i with x f i ∼ 1 [14][15][16]. There also exist some other non-thermal DM production mechanisms such as incomplete reheating [21], asymmetric reheating [22][23][24], and Dodelson-Widrow mechanism [25], however, we do not consider in this work.
In the FIMP scenario, the T f i ∼ O(10 2 ) GeV for the weak scale DM, this temperature is in the ballpark of the electroweak phase transition temperature of the Universe. One can expect that the thermal effects from the electroweak phase transition change the production of the FIMP DM by modifying the kinematical threshold of the DM decay and annihilation processes, recent studies can be found in Refs. [26][27][28] 2 . The electroweak phase transition provides an explanation of the electroweak symmetry breaking of the SM at the early Universe, which can be tested by probing the triple Higgs couplings at high energy hadronic colliders [30] and the future linear colliders, such as CEPC, ILC, and FCC-ee. On the other hand, supposing the electroweak phase transition is a strongly first order phase transition (SFOEWPT), the space-based gravitational wave detectors can probe the gravitational wave signals from the phase transition.
Recently, the hidden sector triggered multi-step phase transition arouse peoples interests, see Ref. [34][35][36][37] for the combined studies of the WIMP DM and the two-step SFOEWPT. Wherein, the first-step is a second-order phase transition, and the secondstep is a SFOEWPT. There also exist a broad class of well-motivated DM models with more than one hidden sectors, wherein the DM particle is one component of the hidden sectors [38,39]. In this work, we study the two-step SFOEWPT scenario with the FIMP DM particle being one component of the hidden sectors. The non-zero neutrino mass is 1 The population of DM for the initial negligible DM abundance is related to the deep understanding of the early universe. Due to the FIMP DM never enter into the thermal equilibrium, primordial perturbations in the DM density spectrum would not be washed out and may leave an imprint in the CMB, especially when the FIMP DM is the scalar inflaton [16,17]. This feature is supposed to limit the FIMP models [16]. The interaction strength between the inflaton and the DM sector might be bounded by the DM relic density [18][19][20][21]24]. We left the connection of this work with the inflation and the reheating to a separate publication. 2 For the phase transition modified WIMPs scenario, we refer to Ref. [29].
well established by the oscillation experiments, which requires for new physics beyond the SM [31]. 3 Therefore, it would be inspiring if the neutrino mass can be addressed together with the two-step SFOEWPT modified FIMP DM with the DM being a part of hidden sectors. This constitute the essential ingredient of this work. The scotogenic model proposed by E. Ma [43] falls into the context of the DM models, where the Dirac neutrino mass is radiatively generated with the agumented of two hidden sectors (including an inert scalar doublet and a scalar singlet). 4 The DM particles in the model can be the lightest of the heavy neutrinos or the neutral scalars from the dark sector. 5 This work is organized as follows: We firstly review the scotogenic model in Section 2. The phase transition patterns and the implications for the FIMP DM production are explored in Section 3. The DM phenomenology is studied in Section 4. The gravitational wave signals predictions and related collider interplay are discussed in Section 5 and Section 6. We conclude with Section 7.

The scotogenic model
We revisit the scotogenic model that generates the Dirac neutrino mass with dark matter at one-loop level [43], the Yukawa interactions for the radiative Dirac neutrino mass generation are given by where the introduced inert doublet (η) together with an additional real singlet scalar (S) are all odd under the Z 2 symmetry, the tree level potential of the model is In the global minimum of the Electroweak (EW) vacuum, the doublets are given as where Φ develops a VEV v = 246 GeV. The Z 2 symmetry of η, S remains unbroken in the EW vacuum. Goldstones G + and G 0 are eaten by W ± , Z bosons after the spontaneous symmetry breaking (SSB). We consider CP conserving situation in this work and all soft masses and scalar quartic couplings in Eq. 2.2 are assumed to be real, the parameter relations and relevant theoretical and experimental constraints are given in Appendix A. In Ref. [46], a one-step SFOEWPT is reachable due to a sizable dimensional-six operator of Higgs after integrating out the heavy scalars in the scalar potential. While in this study we focus on the two-step phase transition since the dimensional-six operator can only be obtained at loop level, see Ref. [47,48], which is different from our previous study on mixed scalar WIMPs DM [44] where a sizable interaction rate between the SM Higgs sector and the dark scalars are allowed.
In the EW vacuum, the mixing between the singlet and the inert doublet neutral scalars S, H 0 connects the gauge eigenstates to the mass eigenstates χ and H by Here, the mixing angle is For parameter relations between mass eigenstates and gauge eigenstates, see Appendix A.1.
with the Feynman diagram given in the Fig. 1. For the scotogenic model proposed in Ref. [42], the dangerous of flavor-changing charged-lepton radiative decays can be avoided by assuming a superweak Yukawa interaction of f ακ , which prefers a FIMP DM [14,40,49]. Indeed, a superweak interaction of N 1 with scalars leads to the decouple of N 1 from the neutrino mass matrix, and there are only two light neutrinos acquire non-zero masses, see Ref. [40]. For the Dirac neutrino mass generation in the scotogenic model proposed by Ref. [43], the extra scalar are introduced to accomodate the second yukawa interaction of the Eq. 2.1. There is no problem of the flavor-changing charged-lepton radiative decays. In this work, we would explore the N 1 as the FIMP DM, in which case the N 1 decouples from the neutrino mass matrix. Furthermore, the smallness of the neutrinos masses can be ensured by the smallness mixing between S and H 0 in the model, which naturally leads to the production of χ through freeze-in mechanism, which then decay to the DM particle N 1 after the freeze-in. This study applies to the Global U (1) B−L scenario of Ref. [43] and the millicharge situation of the gauged U (1) B−L scenario.

The electroweak phase transition dynamics
Supposing H, A, H ± and S are all much heavier than the SM-like Higgs, one can safely integrate out the dark scalars, and therefore the tree-level potential with dimensional-six operators is given by with the c 6 mainly coming from the one-loop level, given by Ref. [47,48] respectively for the dark scalar and inert doublet cases. Therefore, to obtain a strongly first order phase transition with one-step pattern, a relatively larger quartic coupling is required [47]. This may result in perturbativity problem [35,36,47,48] even at the model phenomenological valid energy scales. Another reason that motivates us to focus on the two-step pattern phase transition is to account for the production of scalar χ and the DM N 1 through freeze-in mechanism.

Electroweak phase transition types and vacuum structures
Firstly, we analyse the vacuum structure at zero temperature, which is highly related with the possible phase transition types, i.e., one-step and two-step. The tree-level scalar potential of the classical fields is given by In the space of (h, H 0 , S), the minimization conditions, give rise to the possible minima of the potential localized at We note that the possible minima where h and S(or/and H 0 ) have VEVs is precluded because both the inert doublet and the real singlet scalar are all odd under the Z 2 symmetry.
The EW vacuum localized at (± −µ 2 Φ /λ 1 , 0, 0) should be the global minimum, which ensures the possibility of the phase transition type of O → C in Fig. 2. There also could be three other local minima in the direction of (0, H 0 , 0), (0, 0, S ) for the case where the mixing of H 0 , S is negligible, and in the subspace of (0, H 0 , S ) when the mixing of H 0 , S is non-negligible which however is not favored by the neutrino mass generation.
To ensure the EW vacuum to be the global vacuum, the following conditions need to be fulfilled which result in the following two bounds on the parameter spaces, In the spirit of the gauge invariant [37,50], the finite-temperature potential can be described by Eq. 3.2 with the substitutions of Here, we neglected all other Yukawa couplings apart from the top quark Yukawa coupling y t . All three directions of h, H 0 , S can induce minima due to the thermal corrections at finite temperature depending on these scalar quartic couplings. The two-step phase transition can occur through O → B → C (see Fig. 2) when the vacuum along the direction of H 0 appears earlier than that of the h direction, which corresponds to the parameter constraints of 11) or occurs through O → A → C (see Fig. 2) when the Z 2 symmetry in the direction of S breaks earlier than the SSB of EW vacuum in the direction of h, which results in, In this work, we study the above three patterns of phase transition. For the one-step pattern, the phase transition can occur directly from the EW symmetry and Z 2 symmetry phase to the EW symmetry broken phase with Z 2 symmetry. For the two-step pattern of O → A → C, the parameter spaces are restricted by Eqs. (3.7, 3.12), see the left plot of Fig. 3. For the two-step pattern O → B → C, Eqs. (3.7,3.11) restrict the parameters as shown in the right panel Fig. 3.
The bounce configuration of the nucleation bubble ( the bounce configuration of the multi-fields that connects the EW broken vacuum (h− vacuum, true vacuum) and the false vacuum( here it can be Z 2 broken vacuum H 0 (S)− vacuum for the two-step scenarios)) can be obtained by extremizing through solving the equation of motion for φ b (it is h and H 0 /S for two-step scenarios), with the boundary conditions of The phase transition completes at the nucleation temperature when the thermal tunnelling probability for bubble nucleation per horizon volume and per horizon time is of order unity [51][52][53],  Before the study of DM dynamics, we first explore the scalar masses evolving with temperature. The thermal masses (see Appendix B for details) are crucial for the thermal decay width and scattering cross section which determine the amplitude and possible evolution history of the DM relic abundance. From Fig. 4 (the parameters are chosen as in Fig. 11), one can find that in the phase transition scenario of O → C, the thermal masses m H,χ (T ) almost equal to the physical masses at T = 0, therefore, we do not expect a larger deviation of the thermal modified FIMP from the case without taking into account the thermal effects. When the universe goes through the first-stage second-order phase transition and the second-stage first-order phase transition with the dropping of the temperature, m χ (T ) and m H,H ± (T ) show a highly dependence on the temperature for the O → A → C and O → B → C phase transition scenarios, see Fig. 5 (with parameters are chosen as the top-right panel of Fig. 12 to ensure all dark scalars live in thermal bath during the production process of the DM N 1 through freeze-in mechanism) and Fig. 6 (with parameters are chosen as the bottom-right panel of Fig. 12 and assuming a negligible small θ). Therefore one can expect the kinematic threshold of the decay/inverse decay and the scattering processes that contribute to the DM production can be different from the tradi-tional FIMP DM calculations. More specifically, the mass splitting between m χ,A,H,H ± (T ) and m N 1 is dynamical and can be vastly different from the zero temperature case in a long time duration, which may lead to a enormous difference of DM relic abundance between thermal modified calculation and the traditional calculation without taking into account the thermal mass. For the case of χ producing through the freeze-in mechanism, the λ sφ,sη should also be negligible, with m χ (T ) ≈ m χ .

DM phenomenology
For the DM relic density calculation, the accumulated relic abundance of the previous temperature duration will be taken as the initial abundance of the next temperature duration. The DM relic abundance accumulated in the the symmetry phase can be significant depends on the reheating temperature which is taken as the initial condition of the calculation of FIMP where the entropy normalized number density is assumed to be null. For the DM number density calculation with phase transition, we use the thermal corrected mass to replace the physical mass at zero temperature in the decay and scattering processes.

FIMP N 1 DM
In the following, we study the case of m χ,H,H ± ,A > m N 1 without mass degeneracy. The relevant study can be found in Ref. [40], here the difference is that we have additional decay of χ → N ν. As in Ref. [40], the 2 ↔ 2 process is highly suppressed in the situation with phase transition pattern O → A → C, and we therefore focus on the decay/inverse decay of χ → N ν dominated FIMP DM production. We first study the freeze-in contribution to N 1 production, its production will be dominated by the decays of the scalars (H, A, H ± , χ) while they are in equilibrium with the thermal bath. The , can be computed by solving the following Boltzmann equation [14] dY where s is the entropy density of the Universe, H(T ) is the expansion rate of the Universe at a given temperature, with X = H, A, H ± , χ and being a SM lepton. In this equation, K 1 (x) is the Bessel function of the second kind, and g X is the number of internal degrees of freedom of particle X. Specifically, g H,A,H + ,H − ,χ = 1. The decay rates that enter into Γ (X → N 1 ) are calculated as Firstly, the requirements that the N 1 does not reach thermal equilibrium, i.e., Γ(X → N 1 ν) < H(T ), set the upper limits on the Yukawa couplings h 1α , and the lower limits on which is set by the Big-Bang nucleosynthesis (BBN) requirements [41]. At the BBN epoch, the finite temperature effects are negligible. We note that the increase of DM mass m N 1 means a larger parameter spaces allowed by the non-equibrium conditions and the BBN, as shown in Fig. 7. As for the Yukawa coupling of f 1 , one needs to worry about the disturbation to the BBN or CMB by the decay of dark scalars to N 1 , because all these dark scalars are in equilibrium before T f o ∼ m A,H,H ± /26 and these particles' number density quickly drops to be almost zero after T f o . process once the kinematical threshold is allowed. The Boltzman equation for this scenario is here the p γH ± = (s−m 2 H ± )/2 √ s, and the cross section of σ γH ± → N is given in Appendix C. In Fig. 8 is obtained by solving the Eq. 4.6 (all the thermal masses adopted here are the same as Fig. 6). The departure of the "with PT" yield Y N 1 (x) from the "without PT" yield is noticeable, this is because the kinematical thresholds for the decay and scattering process are dynamical during the temperature evolution. In particular, the thermal mass of A, H ± mostly smaller than the corresponding zero temperature masses within the temperature epoch under study. The two plateau behaviors of the solid curves reflect if the thermal masses reach the threshold of the DM production processes during the two stage phase transition.

Freeze-in N 1 augmented by late decay of FIMP χ
In this section, we study the case that DM N 1 is partially generated by freeze-in mechanism, and partially generated by late decay of FIMP χ, in the phase transition pattern O → B → C. When the yukawa couplings f 1α ∼ h 1α O(1) and λ sφ,sη , µ sof t 1, both χ(in this case mostly S component) and N 1 never reach thermal equilibrium(both are FIMPs). Then the coupled Boltzmann equations are The decay widths are given as follows, (4.10) The χ particle mostly comes from the inert sectors decay as given by Eqs. (4.8,4.9), the decouple conditions is given by Γ(H ± /A → χW ± /Z) < H(T ). The left panel of Fig. 9 depicts that the combination of non-equilibrium condition and BBN bounds the mixing angle 10 −11 < θ < 10 −6 depending on the temperature. To obtain the correct magnitude of the neutrino masses, as aforementioned in Sec. 2, the mixing angle can not be infinitely small. For the θ allowed by the non-equilibrium conditions and the BBN, the Yukawa couplings of f 2,3 and h 2,3 and the heavy neutrino masses are restricted, as indicated by Eq. 2.6. So we illustrate the constraints in the right panel of Fig. 9, the mixing angle θ and the heavy neutrinos mass are restricted into a narrow regions for typical values of the Yukawa couplings f α2,3 and h 2,3β .   We consider the χ particle mostly comes from the decay of H ± , then the late decay χ → N 1 ν dominantly contributes to the relic density of N 1 . Before that, the contribution of the annihilation process dominates the production of N 1 , as depicted in the right panel of Fig. 8. In the left panel of Fig. 10 we show the "without PT" case without taking into account the thermal mass, the decay H ± → χW ± is always active during the production of χ particle. For the case with the thermal effect labeled as "with PT" in the right two panels of Fig. 10, the thermal masses m H,A,H ± (T ) are adopted, and the m χ (T ) ≈ m χ since here we need to assume negligible small λ sη,sφ , θ to ensure χ does not enter into the equilibrium. The kinematical threshold of H ± → χW ± can only open at the pretty early stage at high temperature, which leads to a sharp increase of Y χ at the beginning. For the cases with the thermal effects, the correct DM relic density highly relay on the reheating temperature due to the modification of the kinematical threshold by the thermal masses, we plot the middle and the right panels of Fig. 10 to illustrate this situation. In comparison with the case of "without PT", the late decay of χ contributions dominate the DM relic density of N 1 , the much smaller f 1 and h 1 are chosen here to make sure the correct relic density Ωh 2 = 0.12 [32] can be obtained. The magnitude of h 1 is relatively larger than f 1 to ensure the χ particle can fully decay before BBN.

Gravitational waves from the SFOEWPT
The gravitational wave signals produced during the phase transitions process can be characterized by two parameters of α and β at the phase transition temperature T (bubble nucleation temperature). The first important parameter α is the latent heat normalized by the radiation energy, given by α = /ρ rad with ρ rad = π 2 g T 4 * /30. The latent heat includes the difference of the vacuum energy between the false and true vacuum and the entropy variation ∆s (see Appendix D) at the phase transition temperature, given by The second crucial parameter β reflects the duration of the phase transition, and characterize the peak frequency of the GW spectrum. Under the assumptions of adiabatic expansion of the universe, one has where S 3 is the three dimensional Euclidean action for the critical bubble (nucleating bubble). With the two parameter at hand, one can calculate the gravitional wave signals produced by the SFOEWPT, which includes three dominant contributions: bubble collisions, sound waves and Magnetohydrodynamic turbulence (MHD) in the plasma [54,55]. The total energy spectrum is given by 3) The first important contribution is the bubble collision, estimating with the envelop approximation [56][57][58], which is given by [59] with the bubble wall velocity v b and the efficient factor κ being functions of the crucial parameter α [60], 3α/2 1 + 0.715α , and the peak frequency located at Hz . (5.5) The other two main contributions are the sound waves and the MHD, which are given by with κ v ≈ α(0.73 + 0.083 √ α + α) −1 and κ turb ≈ 0.1κ v describing the fraction of latent heat transformed into the bulk motion of the fluid for sound waves and MHD. Today's Hubble parameter is The peak frequencies of sound waves and MHD locate at Hz , (5.9) Hz . For the one-step SFOEWPT pattern O → C, one can consider the energy barrier for the phase transition almost coming from the singlet part or the inert doublet part, since the mixing of the two scalars is assumed to be negligible. For the inert doublet dominant case, we refer to Fig. 11, which is analogy to the study of Ref. [48]. While the singlet dominant case requires large Higgs portal coupling λ sφ , see Ref. [36]. In Fig. 12, we plot the GW signals to be probed by the ALIA, eLISA, BBO, DECIGO and ultimate-DECIGO experiments, which are generated from the two-step SFOEWPT with patterns O → A → C and O → B → C. The top plots of Fig. 12 indicate that the increase of λ 2,L,sη leads to the shift of the frequency to a smaller value. The parameter λ 2,L,sη chosen for the top-right panel ensures that the inert scalars and the singlet scalar live in the thermal bath, which is the benchmark model corresponding to the freeze-in production of N 1 with O → A → C scenario, as studied in left panel of Fig. 8. The related studies on GWs in this scenario can be also found in Refs. [61][62][63][64][65] 7 . The bottom plots show that the peak of the spectrum/frequency shifts to the right by increasing λ s,sη,sφ for the O → B → C pattern. The freeze-in production of N 1 without late decay of χ (the right panel of Fig. 8) could be obtained with the parameter setup shown in the bottom-right panel of Fig. 12, except employing a negligible small θ that will not affect the GW signals production. Considering the freeze-in production of χ, as studied previously in Sec. 4.2 (see Fig. 10), the negligibly small λ hs,sφ and θ are required which also do not affect the phase transition dynamics. This can be realized due to the pattern O → B → C is the inert doublet-like two-step phase transition. The study of the GWs from the inert-doublet two-step phase transition is absent in literature, we fill the blank in this scenario.

Comments on collider searches
One of the signals predicted by the SFOEWPT is the deviation of the triple Higgs couplings, which might be able to probed at colliders, such as the lepton colliders [68,69] and high energy hadronic colliders [30]. We emphasize that this kind of search mostly apply to the SFOEWPT through the one-step pattern with a sizable dimensional-six operators after heavy particles being integrated out. For the phase transition patterns O → A → C and O → B → C, there is no sizable dimensional-six operators could be obtained since there is no contribution from the tree-level. In both cases, the dimensional-six operators can be obtained at one-loop level, see Ref. [47,70] for the O → A → C pattern, and Ref. [48] for the O → B → C pattern. As argued in Ref. [28], both the two scenarios might be tested at LHC with the off-shell Z-pair channel [71], and at the future lepton colliders [44,47,65] as well as 100 TeV pp collider [36]. As studied in Ref. [44], the displaced vertex signature might be able to tell the difference between the scalar sectors of this model and the inert doublet model. We checked the charged Higgs contribution to the Higgs diphoton signal following Ref. [44], the deviation from the SM prediction is negligible for the FIMP benchmark scenarios under study.
For the FIMP DM N 1 , the decay-lengths of χ, H, A, H ± particles are far beyond the ability of LHC. For the scenario of m N 3 > m H ± ,A,H,χ > m N 2 > m N 1 , as studied in Ref. [40], one can expect the relic abundance of N 1 DM coming from the late decay of N 2 after the freeze out and before BBN epoch, the decays of the heavier singlet fermions N 2 → N 1¯ mediated by the scalars (ignore the mixing of H and χ) should happen with Therefore we obtain parameter regions allowed by BBN, which is shown in Fig. 13.The freeze-in mechanism of DM N 1 production is not valid in this scenario due to the largess of the Yukawa couplings.  Ref. [49] studied the LHC dilepton pair signal in the FIMP DM scenario of the scotogenic model proposed by E. Ma [42], which can be applied to our study except the late decay which calls for much heavier m N 2,3 , see Fig. 9. In the model under study, the dilepton decay for m H ± ,A,H,χ > m N 3 > m N 2 > m N 1 is given by where we have an additional decay channel of N 3 → N 2 α β . The decay length is shown in Fig. 14, which means the channel can be tested at the LHC.

Conclusion and discussion
This work demonstrate that the two-step SFOEWPT modifies the thermal history of the FIMP DM production through modifying the kinematical thresholds of the decay and scattering processes. The phase transition effects may result in a larger or smaller DM relic density in comparison with the traditional calculation. We study: (1) the FIMP DM production via the the 2 → 2 scattering processes and/or the decay of the particles from the thermal bath; (2) and the DM production with the late decay of a neutral scalar, which is also produced through the freeze-in mechanism. The hidden scalars that induce the two-step SFOEWPT would lead to detectable signals at colliders and the gravitational wave signals to be probed on the future space-based interferometer. The realization of the second scenario here requires a tinny mixing angle between the inert CP-even neutral scalar and the singlet scalar to accommodate massive neutrinos of the scotogenic model. The SFOEWPT modified FIMP DM is general and apply to the context of DM models with more than one hidden sectors. For the explanation of the baryon asymmetry of the Universe utilizing the electroweak baryogenesis mechanism [72], an extra CP-violation operators involving the dark sector are necessary [63,64,73,74]. For the leptogenesis realization in the scotogenic model and it's connection with the WIMPs DM, we refer to Refs. [75][76][77][78].

A.1 Parameter relations
In the basis of real neutral scalars, (S, H 0 , A), the mass matrix for the dark scalars is The masses and interactions of the scalar sector are parameterized by the scalar quartic couplings and µ Φ,η,S,soft , which are in terms of the physical masses in the global EW vacuum with λ L,S = (λ 3 + λ 4 ± λ 5 )/2.

A.2 Perturbative unitarity
All the scattering matrices of the model are [44] 8πS Y =2,σ=1 =    Then the eigenvalues Λ X Y σi of the above scattering matrices (where i = ± or 1, 2, 3) can be calculated as 13) and the Λ even 00 1,2,3 correspond to the three roots of the polynomial equation (A.14)

A.5 LEP bounds
Considering the CP-even neutral Higgs H is a mixture of doublet and singlet, the LEP bounds on the scalar masses of this model are the same as in IDM. The lower limits on the scalar masses comes from the precise measurements of the W and Z widths, Utilizing the neutralino search at LEP II, the mass splitting ∆m HA should be either smaller than 8 GeV or greater than 100 GeV for m A < 80 GeV [80]. The production of the charged Higgs pairs H + H − at the LEPII give rise to [81], The corresponding parameter limits from the neutrino mass are plotted in Fig. 15, a larger mixing angle θ is preferred for a lower magnitude of f α2,3 h 2,3β . Since we assume the N 1 never reachs equilibrium with the thermal bath, the requirement of the neutrino mass of order ∼ 0.1 eV leads to a ballpark estimate of the size of the multiplication of Yukawa couplings f α2,3 h 2,3β and the soft terms µ sof t .

D Entropy injected by EWPT
In this section, we recall knowledge of entropy deviation induced by EWPT [84]. Assuming that reheating happens quickly relative to the expansion rate, the energy density ρ of the universe does not change during reheating. And there are only small amount of reheating by the release latent heat of the transition. We do not expect the phase coexistence stage for SFOEWPT. The injection of entropy from by supercooling process of EWPT for the two step pattern is evaluated by [34], with the finite temperature potential V included only the thermal mass correction evaluated at the nucleation temperatures, which has been normalized by the SM entropy in the radiation dominate universe, By taking g s = 100 there, it was estimated that ∆s/s EW is around percent level around the phase transition temperature T n for benchmarks of O → A(B) → C. As in the Standard model, (see Ref. [85]), one finds a negligible dilution factor, namely: with the entropy difference between the high temperature and low temperature entropy: ∆s = s + − s − .