Thermally modified sterile neutrino portal dark matter and gravitational waves from phase transition: the freeze-in case

We consider the thermal effects into the evaluation of the dark matter production process. With the assistance of the right handed neutrinos, the freeze-in massive particle dark matter production history can be modified by the two-step phase transitions. The kinematic of decay/inverse decay or annihilation processes can be affected by the finite temperature effects as the Universe cools down. The history of the symmetry respected by the model can be revealed by the DM relic abundance evolution processes. The strong first order electroweak phase transition generated gravitational waves can be probed. The number of extra scalars for the Hierarchy problem can be probed through the Higgs off-shell searches at the LHC.


Introduction
The baryon asymmetry of the Universe (BAU) and the dark matter (DM) are two fundamental physics problems that urge the particle physicists and cosmology physicists to propose variants of intelligent ideas and models. The electroweak baryogenesis mechanism has been studied extensively to solve the BAU problem due to the detectable signals of the strongly first order electroweak phase transition (SFOEWPT) at high energy colliders [1]. A SFOEWPT can generate a detectable gravitational wave signal with a typical peak frequency O(10 −3 − 10 −1 )Hz. Extensive studies of the gravitational waves generated by the SFOEWPT in various new physics models are inspired especially after the first discovery of the merging black holes with the gravitational waves signals detected by LIGO [2].
The existence of the DM is supported by substantial astrophysical and cosmological observations. With the accumulation of the direct and indirect detection experimental data of the DM, the Weakly Interacting Massive Particle (WIMP) paradigm confronts an increasing pressure through the interactions between the dark and the visible sectors which are detectable at colliders. The WIMP DM produced at super low temperature around -1 -

The model
In this paper, we utilize the model similar to refs. [14,20,22,23]. It contains a Z 2 -odd majorana fermion χ and a real-scalar boson φ. The SM-fields are all even under the Z 2 . We also introduce a Z 2 -even sterile neutrino. This sterile neutrino can be majorana or pseudo-Dirac. The sterile neutrino-Higgs-lepton boublet couplings in the pseudo-Dirac case can be in a much wider range than the couplings in the majorana case, and the freeze-in process is not much disturbed by details of the sterile-neutrino sector. In this paper, we focus on 1 For the WIMP DM modified by the phase transition, we refer to ref. [13]. 2 For the FIMP DM without thermal effects in the similar model we refer to refs. [18][19][20]. the pseudo-Dirac case. The general Lagrangian is given by where χ C = χ is the four Dirac four-spinor, l i with i = 1, 2, 3 or i = e, µ, τ are the SM left-handed lepton doublets.N D = N 1 iσ 2 N * 2 is a Dirac four-spinor, and the N 1 and N 2 are the Weyl-components. m D is the Dirac mass terms, and µ 1,2 are the Majorana mass terms, m χ,φ,N are the χ, φ, N mass terms, and y χ,χ5,N i , λ φ,φH are the coupling constants. We rotate to a basis that the y χ,χ5,N i , λ φ,φH and m χ,φ,N are real numbers. For simplicity, we omit the y χ5 which breaks the CP symmetry. In this paper, we adopt the convention in the left-handed lepton and Higgs doublets where G + , G 0 are the goldstone bosons, h is the SM Higgs boson, and v = 246 GeV. Note that A · B = A i (iσ 2 ij )B j , where σ 2 is the second Pauli-matrix. Although in this paper, we only introduce one pseudo-Dirac sterile-neutrino in our calculations, and this seems to be different from some common see-saw models, we can always rotate the mass basis in the m N,D ∝ I case so that only one sterile neutrino interact with the dark matter, see appendix A. Therefore, our calculations are still available. In general, the "pseudo-Dirac" particle with the nonzero µ 1,2 and y N Ci terms actually split into two nearly-degenerate majorana components. However, these terms are usually rather small in both the linear-and inverse-seesaw models, and their effects in the early universe are usually negligible. Therefore, for convenience, we just set all of them to be zero during the calculation processes.

Electroweak phase transition dynamics
In this work, we study the effects of a two-stage phase transition on the DM production. The first-stage of the phase transition is a second-order phase transition from the symmetric phase (both Z 2 symmetry and the EW symmetry are all preserved) to the phase where the hidden scalar φ got VEV (where the Z 2 symmetry breaks). The second-stage of the phase transition is a first-order phase transition from the Z 2 symmetry broken and EW symmetry preserved phase to the Z 2 symmetry restored and EW symmetry broken phase. The phase transition can be characterized by the thermal evolution of the thermal potential. The tree level scalar potential at zero temperature is given by,

JHEP12(2018)006
The desired vacuum structure for the two-step phase transition can be obtained easily with the local vacuum localized (h, φ) = (0, v φ ) and the global Electroweak vacuum localized at (h, φ) = (v, 0), Which give raise to µ 2 = λv 2 and v φ = −µ 2 φ /λ φ . Then, the mass of the hidden scalar φ is given by: Firstly, the local vacuum at (h, φ) = (0, v φ ) should be higher than the global one at which constrains the parameters as, The electroweak phase transition can be studied in a gauge invariant approach [24] after taking into account the thermal corrections, the thermal potential used to estimate the vacuum structures at finite temperature is given by, with the contribution of fermions χ negligible since we focus on the freeze-in DM scenarios.
To ensure the second-stage strong first order phase transition to occur, one needs the vacuum in the direction of φ (at temperature T φ ) to appear earlier than the one in the direction of h(at temperature T h ) during the phase transition, i.e., T φ > T h , Then, the eq. (3.3), eq. (3.4), and eq. (3.7) together bounds the parameter spaces of λ φ , λ φh and m φ , as shown in figure 1. To evaluate the DM at different temperatures, the field dependent mass at zero and finite temperatures are crucial, the zero temperature field dependent mass are, In the Electroweak vacuum, the thermal corrected mass are given by, (3.14) for Goldstones, and for m 2 h (T ) < m 2 φ (T ) and m 2 h (T ) > m 2 φ (T ). The VEVs of the Higgs h and scalar φ in the h-vacuum and the φ-vacuum as a function of temperature are given by, Which leads to the variation of particle masses that enter into the DM production process, see the section 4. As the temperature cools down, the scalar masses are essential for the evaluation of the evolution of the DM number density. Prior to the second stage first order phase transition, all particles are set in the φ-vacuum(wherein the Z 2 symmetry is broken and the SM EW symmetry is preserved), and there is no Electroweak symmetry breaking, the scalar field dependent thermal corrected mass are,

JHEP12(2018)006
Assume an instantaneous transition from the electroweak symmetric phase (φ-vacuum) into the electroweak symmetry breaking phase (h-vacuum, wherein the Z 2 symmetry of φ is restored and the SM EW symmetry is broken), i.e., a instantaneous supercooling process, we have the scalar field dependent thermal corrected mass being, To evaluate the assumption, we calculate the bounce solution to find the nucleation temperature of the bubble T n where the phase transition from the φ-vacuum to h-vacuum occurs. This makes the assumption established for the DM number density evaluation, the final DM number density in the φ-vacuum will be taken as the initial number density in the h-vacuum. Therefore, we directly cast the DM number density in the φ-vacuum to the h-vacuum around the second stage first order phase transition. When the Universe cools down to the nucleation temperature T n , which is pretty close to the critical temperature, the phase transition from φ-vacuum to h-vacuum proceeds. At this temperature, the bubble nucleation rate per unit volume per within the Hubble horizon H −1 reaches unity, one can obtain the bounce configuration of the two-fields h and φ that connects the EW broken vacuum (h-vacuum, the true vacuum) and the Z 2 broken vacuum (φ-vacuum, the false vacuum) through solving the equation of motion for h and φ, with the boundary conditions of where φ b represents the h n and φ n . In particular, we adopt the method developed in refs. [26][27][28]. Before the study of the DM in detail we firstly study the entropy induced by the EWPT. We recall the knowledge of entropy deviation induced by the EWPT following ref. [29]. 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 amounts of reheating by the release latent heat of the transition. We do not expect the phase coexistence stage for SFOEWPT.
The injection of the entropy from supercooling process of EWPT for the two step pattern is evaluated by, with the finite temperature potential V included only the thermal mass corrections evaluated at the nucleation temperatures, which has been normalized by the SM entropy in the radiation dominate universe, (3.28) By taking g s = 100, it was estimated that ∆s/s EW is about percent level around the phase transition temperature T n for both benchmarks of figure 9. With the entropy difference between the high temperature and low temperature entropy: ∆s = s + − s − , one have The figure 9 indicates one obtain a negligible dilution factor as in the SM case, see ref. [30], and the two-step phase transition pattern being studied in ref. [31].

Thermal effects modified dark matter production
The calculations of the relic abundance of the thermally modified dark matter are based on the Boltzmann equations (We derive the equations according to refs. [32,33]). For completeness, we consider both the case of m φ > m χ + m N (D) and m φ < m χ + m N (D) . Keeping in mind that, with the temperature cooling down, the phase goes through a first-stage second-order phase transition and then a second-stage first-order phase transition with the vacuum tunneling from the φ-vacuum to the h-vacuum. We assume the second stage phase transition from φ-to h-vacuum to occur instantaneously without large reheating.

JHEP12(2018)006
The Boltzmann equations in the Z 2 symmetry phase and the h-vacuum are given by where the Y X = n X s is the actual number of the constituent X per-comoving-volume, and the Y Xeq = n Xeq s is the equilibrium number of the constituent X per-comoving-volume, n X(eq) is the (equilibrium) number density of the constituent X, s is the entropy density, z = mχ T , and T is the temperature, H is the Hubble constant. Here, we note that the reheating temperature is supposed to be above the φ-vacuum to h-vacuum phase transition temperature, which is taken to be the bubble nucleation temperature as been estimated in figure 9 of section 5.1.

JHEP12(2018)006
Firstly, we define the thermal decay width ofΓ φ→χN as with the thermal mass m φ (T ) calculated as in the previous section and If one have m χ > m φ (T ) + m N temporary during the phase transition process, we need need to replace the Γ(φ → χN ) by the decay width of the χ, i.e., Γ(χ → N φ), being given by For completeness, we consider also the decay (or inverse-decay processes) N → h ±0 l( or h ±0 → lN ). After the phase transition of φ-to h-vacuum, we need to consider the processes N ↔ W ± /Z/h also. As in ref. [34], we calculate the N → h ±0 l or h ±0 → lN with all the four states of the Higgs doublets having the Higgs boson mass m h (T ). If As the temperature cools down to T m h (= 125 GeV), the magnitude of theΓ N is highly suppressed by a factor of e −2m h +m N T . In order to let the right-handed neutrino decay, we usẽ being calculated at the zero temperature. We further note that in this work, the thermal history of the sterile neutrino does not significantly affect the dark matter production. In this work, we restrict our interest to the mass of N being smaller than 80 GeV. In this case, these two body decays are kinematically forbidden. Meanwhile, three body decays is dominant by the process involving the exchange of a virtual W [35], (4.8) lnx}. An overall factor of 2 appears in the bracket for the N being Majorana case. This factor does not appear in the Dirac case due to the interference of the two components decay. As for the pseudo-Dirac particles, the splitting of the two components will remove the interference, leaving a factor between 1 and 2. In this paper, we adopt the factor 2, but practically, the detailed value does not affect the qualitative, and most of the quantitative results. Secondly, the thermally averaged cross section times velocity σv AB→CD is given by (4.9) where δ CD = 1(0) if C and D are identical(different) particles, g A and g B are the degrees of freedoms of particle A and B, K 1 is a Bessel function, σ AB→CD (s ) is the cross section of the process AB → CD with the total energy in the center of mass frame is √ s . We note that all thermal effects being explored in section 3 has been implemented here.
On the other hand, considering the reheating temperature to above the temperate of the first-stage second-order phase transition, there should be other terms arising from the breaking of the Z 2 symmetry, The number "4" appearing in the σv hχ↔hN indicates that in the φ-vacuum, all the charged and neutral four Higgs bosons are identical. We should note that an on-shell t/u-channel particle appears in many of the processes in the (4.10), due to the on-shell decay of the initial particles. Therefore, when m χ > m N + m φ or m φ > m N + m χ , we do not include all the terms in the (4.10). And that all terms involving the hφφ vertices should be take into account since the EW symmetry is respected in the φ−vacuum.
To solve the Boltzmann equations of eqs. (4.1), (4.10), we use the ready-made function [36,37] embeded in the micrOMEGAs [38] for computing the stiff equations eqs. (4.1), (4.10). With the model file being prepared by FeynRules [39], after implemented the thermal masses and VEVs of φ and h as a function of temperature as in  section 3, we calculate all the thermal σv (s) and the thermal decay widths of the particles using the CalcHEP [40] embeded in the micrOMEGAs. Since the first step phase transition to the φ-vacuum is second order, one can expect negligible entropy injection, the strongly first order phase transition of φ-vacuum to h-vacuum injected entropy is estimated to be around precent level in section 3. The degree of freedoms g * and g * S implemented in the micrOMEGAs are adopted in order to calculate the Hubble constant H = 1.66 M pl , and entropy s = 2π 2 45 g * S T 3 . Here the Planck scale is M pl = 1.22 × 10 19 GeV.

On the decouple conditions bounds on the parameter spaces
For the FIMP production of χ, that would never reach equilibrium, which can be used to estimate the bounds on the couplings through the decouple condition Γ < H. For the decay or inverse induced FIMP, we just needs Γ χ < H, with decay width being given by eq. (4.2) and eq. (4.3). Which results in the bounds on the y χ as shown in figure 3. The figure illustrates that for FIMP production process, the typical freeze in temperature T f i (∼ m χ /z with z ∼ 1) will be smaller than the phase transition temperature T n . Therefore, one can expect the phase transition process will affect the FIMP process.

Dark matter
Firstly, we note that the sterile neutrino can become in thermal equilibrium with the thermal bath through its decay/inverse decay shortly after the reheating process for the typical y N > 10 −6 . And we note that the values of the y N does not have a strong impact on the thermal history of the dark matter. However, since there are only two couplings y χ and y N which connect the sterile neutrino with the SM sector, and usually y χ y N , the thermal history of the sterile neutrino is mainly decided by y N . We show in the figure 4 the thermal evolutions of y N = 10 −5 , 10 −6 , 10 −7 respectively. We can see that 10 −6 is some critical value, below which the sterile neutrino evolution will be similar to a "freeze-in" process, and its decay will be delayed significantly.
Since the SFOEWPT prefer a moderate λ hφ , with which the φ particle can be produced by freeze out mechanism, we study χ freeze in scenario to reveal the EWPT effects. Before the study of DM production details, we explore the mass threshold related to the following DM analysis as the temperature cools down, see figure 5. Before the Z 2 respected by φ and EW symmetry is broken, we denote the symmetric phase in the figure, the m φ (T ) is dominated by the thermal corrections at temperatures higher than the second order phase transition temperature of Z 2 (we denote T φ ) where there is no VEV for h or φ field. In this symmetric phase, one have m φ (T ) = m sys φ (T ) and v h,φ = 0. During the temperature of T φ and the strongly first order phase transition from the Z 2 broken EW symmetry phase (φ−vacuum) to the Z 2 preserved EW broken phase (T = 0) = 246GeV. Our study shows that, a larger quartic coupling of λ hφ and a larger m φ will leads to a smaller T φ , means a later happening of the first-stage second-order phase transition.

φ-DM
In this section, we consider m φ +m N < m χ , in this case we have φ as dark matter candidate. When one consider the φ produced from the FIMP mechanism as in ref. [41] the SFOEWPT fails. The relic abundance of φ can be obtained from the late decay of χ → φN after χ  is generated, while the freeze out mechanism can also contribute significantly when the φ couples with the SM Higgs moderately for a successfull SFOEWPT. One can simplify the evaluation of the relic density of the φ from the two contributions, one from freeze-out and the other from the late decay of χ. with Note that the productions of the dark matter can be separated into different stages, which is similar, but much simpler than the stages in the next subsection. For our benchmark point, since the m χ is much larger, after the first-stage second-order phase-transition, the m φ never exceeds the m χ −m N . Therefore the structure in the figure 6 are much simpler. Finally, when the z 100, most of the χ decays to the φ. For the reasons to be described below, the direct detection experiments do not favor the φ-dark matter case. Therefore, we do not illustrate the detailed steps in this subsection, and leave the descriptions in the χ-dark matter case below.
Generally, a small fraction of DM relic abundance from the freeze-out mechanism calls for a larger λ hφ that is needed for a SFOEWPT, while a larger λ hφ can easily been excluded by the direct detection experiments of Xenon 1T [42,43], LUX [44] and Panda X-II [45] for our interesting DM mass range, see also ref. [46]. On the other hand, a larger fraction of DM relic abundance from the freeze-out requires a relatively small λ hφ which fails the SFOEWPT. This ambiguity almost rules out the possibility of the φ as dark matter.

χ as DM
In this section, we study both the case of m φ > m χ + m N and m φ < m χ + m N , with χ being DM candidate. Eqs. (4.1), (4.10) tells us that both the 1 ↔ 2 and the 2 ↔ 2 channels contribute to the freeze-in processes. Practically, the χχ → N N contributions are highly suppressed by the extremely small coupling constants on their extra vertices compared with the 1 ↔ 2 processes. The dominant 2 → 2 process appears to be χφ → N h, which makes a comparable contribution with the 1 ← 2 processes in the h-vacuum, since its extra vertex is proportional to the λ hφ v, which can be relatively large. The χ freeze-in processes can be affected due to the kinematical threshold can be changed by the thermal effects. We use two typical benchmarks to show the thermal effects in the freeze-in process. The difference is if the 1 ↔ 2 processes is active or not when one do not take into account the thermal effects. We show the two scenarios in top and bottom panels of figure 7. We show the traditional calculations results for comparison when the thermal effect is ignored, denoted as "No Phase Transition" in the figure. The thermal effects modified scenario are denoted as "With Phase Transition". We first explain the physical picture of the thermal effects modified scenario shown in the top-left panel (where the 1 ← 2 processes is always active when the thermal effects are not considered, say the "No Phase Transition" case.) in details: • When z 1, as the temperature is extremely high, v φ (T ) = 0 and φ acquires a significant thermal mass through the eq. (3.6). This induces m φ (T ) > m χ + m N , and the main contribution to the freeze-in processes is the φ ↔ χ + N channel.
• The mass of the φ decreases as the temperature drops. Once |m χ − m N | < m φ (T ) < m χ + m N , the 1 ↔ 2 processes are closed. Much smaller 2 ↔ 2 channels dominates the freeze-in processes. This causes the first dip at the left panels in the figure 7.
• The mass of the φ continues to decrease until m φ (T ) < m χ − m N , 1 ↔ 2 channels reopen, however become χ ↔ φ + N for this time. Therefore, dY χ /dz arises again as can be learned from the figure 7. Note that after this period, a second-order phase transition happens and v φ (T ) becomes nonzero. After this phase transition, m φ (T ) will rise again as v φ (T ) arises.
• As the m φ (T ) arises and reaches the |m χ − m N | < m φ (T ) < m χ + m N range again, there becomes a second dip at the left panels in the figure 7. 2 ↔ 2 channels dominate the freeze-in processes again, however due to the absent of the Z 2 symmetry on this stage, the mixing between the φ and the SM-Higgs boson introduces much larger SM coupling constants in the 2 ↔ 2 diagrams. Therefore, this dip becomes much shallower than the first one.
• When the m φ (T ) arises above the m χ + m N again, the φ ↔ χ + N opens again and the dY χ /dz recovers. However, in the figure 7, there is a third fake dip before the final recovery. This is because we have dropped all of the 2 ↔ 2 channels with an on-shell t/u above the φ ↔ χ + N , which is not a very good approximation around the thresholds. Fortunately, this happens within an extremely small period of time, so the final relic abundance results will not be seriously disturbed.
When the temperature drops to T n , a first-order phase transition happens.
• A final discontinuity/dip at the left panels of the figure 7 appears due to the firstorder phase transition, where one have |m χ − m N | < m φ (T ) < m χ + m N and the 1 ↔ 2 process are closed, and again a sub-dominate process φχ ↔ N h take a role here. After the m φ raise above m χ + m N a tiny increase of dY χ /dz shows up due to the 1 → 2 process is active.
• As can be imagined, with the decrease of the number density of the φ, as shown in the top-right panel of figure 7, one have a smoothly drop of dY χ /dz. After φ freeze out, one have a decay of φ → χN , which leads to the tiny smooth uplift of dY χ /dz.

JHEP12(2018)006
In the top-middle panel of figure 7, we plot the corresponding thermal abundance evolution of the DM χ, where one can find the thermal effects modified case denoted as "With Phase Transition" reveals the phase transition history of the model, including the second order phase transition of φ and the strongly first order phase transition of h. The results of "With Phase Transition" and "No Phase Transition" are close to each other, which is because that the phase transitions of second order and first order occurs around z ∼ m χ /T n ∼ 0.5 and the final and the largest increase of the Ω χ h 2 occurs after the phase transitions.
We show the evolution history in the bottom panel of figure 7 for the case of m χ +m N > m φ , and m χ < m φ . In this scenario, the φ ↔ χ + N process is kinematically prohibited in the h-vacuum shortly after the first-order phase transition. Although φ → χ + ν, where ν indicates a active neutrino, becomes the dominant channel, this is suppressed by the squared mixing angle of light and sterile neutrinos θ 2 (see eq. (A.2)), and we do not include its negligible contributions in our calculations. All the thermal steps are basically similar to the top panel, except that shortly after the last sharp leap in the left panel, where the first-order phase-transition happens, the dY χ /dz rapidly drops after m φ (T ) becomes larger than m χ − m N . Then, dY χ /dz is dominated by the φχ ↔ N h process and smoothly drops down. Notice that since the φ ↔ χ + N processes are absent shortly after the first-order phase-transition (as can be found in figure 5 when one have m χ + m N = 130 GeV (the parameter set of the bottom panel of figure 7) instead of 70 GeV (the parameter set of the top panel of figure 7)), we need a larger coupling constant y D compared with the figure 7 for a correct dark matter relic abundance. Therefore, one can expect a larger discrepancy between the DM relic density calculated with and without the thermal effects. We stress that, different from the first benchmark shown in the top panel, • The first step second order phase transition of φ and the first order phase transition of h all happens around the z ∼ m χ /T n ∼ 1.
• A large increase of Ω χ h 2 happens before the second order phase transition of φ in comparison with the top panels case, mostly due to a larger y D and z ∼ m χ /T n ∼ 1.
• A significant accumulation of the thermal abundance of Ω χ h 2 happens before the decrease of the number density of φ, which can be found through comparison of the bottom-middle and bottom-right panel plots.
• Another difference of the two benchmark scenario can be found through the comparison of the top-left panel and the bottom-left panel, no uplift of the dY χ /dz shows up in the bottom panel's benchmark due to the φ ↔ χ + N process is not active after the φ freeze-out for the case of "With Phase Transition" (indeed the φ ↔ χ + N process only active before the second-stage first-order phase transition) and never active for the case of "No Phase Transition".
We note that for the m N + m χ > m φ case, the suppressed χ + ν decay width will be very likely to exceed the BBN starting time scale ( 1s). In this case, the following decay of the freeze-out φ can potentially disturb the BBN and even CMB by its decay and injecting particles into the plasma during these epochs. However, if m φ > m h , and the freeze-out abundance Ω φ h 2 can be significantly reduced to be 10 −5 (as shown in figure 8), such an influence can be highly reduced by the small φ abundance. In the figure 8, the dY χ /dχ had completely disappeared shortly after z>10. This is because we have omitted the φ → χ + ν channels in our calculations, which depend on θ 2 and does not affect the final results of Ω χ h 2 . The same reason that result in the "With Phase Transition" behavior of Ω χ (see bottom panels of figure 7) leads to the large discrepancy between the DM relic densities (Ω χ h 2 ) of "With Phase Transition" and "No Phase Transition" scenarios (see the middle plot of figure 8).

Comments on the SFOEWPT and the FIMP DM
A hint from the previous benchmarks is that there can be one tight connection between the phase transition and the produce of the FIMP DM. In this section, we demonstrate the connection and the possibility to search the feature with GW signals and at colliders.

Gravitational wave signals
One of the crucial parameter for the Gravitational wave is the strength of the phase transition, the parameter α. Which is the energy budget of SFOEWPT normalized by the radiative energy, being defined by, where the radiation energy of the bath or the plasma background ρ R is given by The ∆ρ is the released latent heat ( vacuum energy density or energy budget of SFOEWPT) from the phase transition to the energy density of the radiation bath or the plasma background. This is given by the difference of the energy density between the false (here it is φ vacuum, ρ(φ n , T )) and the true vacuum (the h-vacuum or EW broken vacuum, ρ(v n , T )),

JHEP12(2018)006
Another crucial parameter of β characterizes the inverse time duration of the SFOEWPT and thus the GW spectrum peak frequency, with H n being the Hubble constant at the bubble nucleation temperature T n . The gravitational wave signals generated by the SFOEWPT mainly include three sources: bubble collisions, sound waves and Magnetohydrodynamic turbulence (MHD) in the plasma [47,48]. The total energy spectrum of the three sources is given by, (5.6) The first source of the gravitational waves from the bubble collision estimated using the envelop approximation [49][50][51] is [52], where the bubble wall velocity v w and the efficient factor κ that characterizing the fraction of latent heat deposited in a thin shell, are all functions of the parameter of α [53], and the peak frequency f env is, Hz . (5.8) The second and the third important sources of the sound waves and the MHD are, .
Here, the fraction of latent heat transformed into the bulk motion of the fluid for sound waves and MHD are given by κ v ≈ α(0.73 + 0.083 √ α + α) −1 and κ turb ≈ 0.1κ v ; the peak frequency of sound waves and MHD are, Hz , (5.10) Hz ; (5.11) the Hubble parameter at present is given by h * = 1.65 × 10 −2 mHz T * 100GeV g * 100   Figure 9. The GWs generated by the SFOEWPT from three main sources: sound waves(brown dotdashed line), collision(blue dotted line), turbulence(orange dashed line) and total contribution(cyan solid line). The colored regions depict the experimental sensitivities of eLISA(two configurations with notation NiAjMkLl), ALIA(gray), BBO(green), DECIGO(yellow) and Ultimate-DECIGO(purple) In figure 9, we show the possibility to probe the benchmarks of the figure 7 and figure 8 with the GW detectors. The sensitivities for the proposed space based interferometers: LISA [54] with two design configurations in notation NiAjMkLl [47,55], BBO, DECIGO (Ultimate-DECIGO) [56] and ALIA [57] are shown with different color shaded regions. For the GW signals produced by the SFOEWPT around the benchmark of the figure 7 is show in the top panels of figure 9. A larger scalar quartic coupling λ hφ can leads to a larger ∆ρ, and therefore the corresponding GW signal from the SFOEWPT is of high magnitude and can be probed easier. We note that a tiny increase of λ hφ drives a lower peak frequency. In the lower panel of the figure 9, we plot the GW signals from the SFOEWPT around the benchmarks of figure 8, the feature is that a smaller λ φ leads to a lower peak frequency and higher magnitude of the spectrum, which is because that a smaller λ φ can induce a larger ∆ρ and a lower β/H. The study of dynamics of phase transition and the produced GW signals with the same scalar sectors under Z 2 symmetry can be found in refs. [58][59][60].

On the GW signals and FIMP DM production
From the side of DM, one ingredient is the open/close of the 1 ↔ 2 process at zero temperature, which is determined by the mass threshold at the zero temperature. If it is not kinematic allowed at zero temperature, the thermal effects would make it active before or during the phase transition process, and therefore results in a significant difference thermal relic abundance in comparison with the traditional one without taking into account the -19 -JHEP12(2018)006 thermal effects. Another ingredient is that when there is a larger portal coupling, there can be a larger contribution to DM production from the annihilation cross section of 2 ↔ 2 process χφ ↔ N h.
Furthermore, a larger m φ ( or λ hφ ) leads to the second order phase transition to occur much latter (or earlier) that can decrease (or increase) the contribution of decay/inverse decay process to the FIMP DM productions, especially when the decay/inverse decay channel is not open at zero temperature as the benchmark scenario of the bottom panel in figure 7. This is due to the interval between the reheating temperature T R and T φ is decreased (increased), and therefore a decrease (increase) of the DM abundance's accumulation.

Collider interplay
In the scenarios being explored in this work, we do not expect to search for the χ field at colliders due to its extremely long decay length far beyond the scope of the detectors. The other two new particles beyond the SM are the dark scalar field φ and the sterile neutrino that are relevant for the neutrino mass generation.
With the gauge invariant approach by taking into account the tadpole contributions to the Higgs two-point green's function [61][62][63], the φ introduced here can be related with the fine-tuning of the Higgs, the quadratic corrections from the hidden scalar singlet is given by, Supposing there are N scalars of φ, then we have, (5.14) We show the value of δm 2 h /Λ 2 in the panel of N φ and λ hφ , see figure 10. In the situation, we can expect the FIMP DM production to be more efficient, and to produce the correct DM relic abundance one just need a y N φ ∼ y χ / √ N . More exactly, y N φ χ ∼ y χ / √ 18 in the parameter region where one can obtain the cancellation of δm 2 h to alleviate the Hierarchy problem where the GW signals can be probed (in this work we use the benchmarks with λ hφ ∼ 0.5) with the phase transition not to be affected by the number of scalars for the two-step pattern [64]. The GWs signatures being explored in this work mostly focus on the case of m φ > m h /2, and therefore, we expect the benchmarks can be probed by the off-shell Z-pair search at LHC [21]. For the search of this scenario at future linear colliders including the CEPC, ILC, and Fcc-ee we refer to refs. [64,65]. The study of ref. [65] illustrate that the Fcc-ee and 100 TeV pp collider are complementary and of great potential to probe this kind of phase transition, the numbers of N φ might be able to be determined. The study of ref. [66] shows that a 100 TeV collider can close most of the parameter space allowed in the SM extended by a scalar singlet when the singlet is allowed to mix with the Higgs but is not heavy enough to resonantly decay to two Higgs bosons.  For the typical mass of the RHN being adopted in this work, we have the decay length being given by [67], For the typical yukawa y N and mixing angle between the sterile neutrino and active neutrino (θ ij ), with the mixing angle being around θ 2 ij ∼ O(10 −11 − 10 −9 ) the decay length is estimated to be cτ N ∼ O(10) m, which is beyond the scope of LHC trilepton and lepton jet search performed in ref. [67]. The proposed SHiP experiment is also unable to probe the scenario due to its proposed mass range [68].

Conclusions
In this paper, we implement the thermal history of symmetry change with the Universe cools down (phase transition) into the dark matter production history. With the assistance of right handed sterile neutrino that lives in thermal bath when the FIMP DM is producing, the dark matter production process can be multi-step due to the effects of the EWPT. More precisesly, the thermal history of the EWPT can change the kinematic threshold of DM production process because the thermal corrected mass and the VEVs of h and φ are all depends on the temperature. We studied the two-step EWPT impacts, all the phase transition history can be revealed by the relic density evolution with the temperature drops, including the second order phase transition of the Z 2 symmetry and the first order phase transition from the Z 2 symmetry broken phase with the EW symmetry to the EW symmetry broken phase with the Z 2 symmetry being respected. Here, both the decay length of the FIMP DM and sterile neutrino are beyond the scope of the present colliders.

JHEP12(2018)006
The interaction rate for the two-step phase transition, i.e., the quartic coupling between the Higgs and the extra scalar φ, can results in the gravitational wave signals to be probed by the future gravitational wave experiments. To alleviate the Hierarchy problem the number of φ should be around ∼ 18 with a y N φ χ ∼ y χ / √ 18 to yields the correct relic abundance through thermal corrected FIMP mechanism. To address the puzzle of BAU within the EWBG mechanism, an extra CP violation source from the high dimensional operators are necessary, see refs. [59,[69][70][71][72].
Although the reheating and inflation is beyond the scope of this work, we briefly comment on cosmological signatures of the FIMP DM. In the FIMP paradigm, the feeble interactions with the thermal bath is assumed to populate the DM, until the DM abundance freeze-in [3,5,73]. A general feature of the FIMP is the absence of DM thermalization with the SM, such that primordial perturbations in the DM density spectrum would not be washed out and may leaves an imprint in the CMB. Future probes of the primordial tensor perturbation spectrum would be able to constrain freeze-in models via the isocurvature perturbations observations [4]. For the asymmetry reheating scenario, ref. [9] shows that the inflaton mediations would actually bring all the fields (hidden and the SM) to thermal equilibrium. For the case where the DM is populated by reheating, the DM may also undergo dark freeze-out without contact with the SM bath [74][75][76]. Supposing there are interactions between inflaton and the DM sectors, the requirement of the correct DM abundance can set limits on the inflaton-DM couplings [6,9,73,77,78].

JHEP12(2018)006
sterile neutrino case, a general mass matrix is given by where all the sub-matrices are 3 × 3 mass matrices. m 3×3 N D , µ 3×3 1,2 are the 3 × 3 mass matrix extended from the corresponding terms in the eq. (2.1). m Dij = y 3×3 N ij v/ where i = e, µ, τ , and j = 1, 2, 3. The electroweak precision measurements constrain the θ ij parameters. The most stringent bounds originate from the FCNC processes, and thus constrain the off-diagonal elements in the m D . In the literature, people are interested in the cases when m N,D ∝ I, therefore the significant FCNC processes are avoided. In this case, all the sterile neutrinos have a unified mass and decay width.
B Two stage phase transition and mixing between χ and N For the case that φ get a VEV during the phase transition, one have the mixing between the sterile neutrino and the χ. In the Majorana sterile neutrino case, N χ = cos θ sin θ − sin θ cos θ with mixing angle being given by tan(2θ) = 2y χ v φ (T )/(m N − m χ ). In the pseudo-Dirac case, things are a little bit complicated. Written in the basis of two-component Weyl spinors, the mass complete mass matrix is given by where χ w is the Weyl component of the four-spinor χ. We do not need to fully diagonalize the (B.2). We only need to rotate out all the terms between the χ and N D components. Therefore, we can still treat the rotated sterile neutrino-like component of the fermions as -23 -

(B.7)
After the vacuum transit from the φ-vacuum to the true Electroweak h-vacuum, these channels are shut down, and the decay channel is replaced by φ → χN , and the annihilation channel are mostly from the t/u channel HH(φφ) → χχ. We comment that due to the mixing angle θ is very small, for the FIMP production of χ with a y χ ∼ 10 −12 , these contributions are negligible.
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.