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.


I. 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 T f o ∼ m DM /26 where the DM particles decouple from the thermal plasma. An alternative approach is the feebly interacting massive particle (FIMP) production mechanism [3], where the DM particles never reach thermal equilibrium with the SM bath in the whole cosmology history due to the small interaction rate with the Standard Model particles. The DM particle freeze in occurs at a relatively higher temperature around T f i ∼ m DM /x f i with x f i ∼ O (1 − 5) in the FIMP paradigm. Considering the SFOEWPT occurs around T n ∼ O(10 ∼ 10 3 ) GeV which is indeed around the T f i when the m DM ∼ O(10 ∼ 10 3 ) GeV, we can expect that the production of the FIMP DM would be significantly modified by the SFOEWPT. This is due to the kinematical threshold that can be altered by the thermal correction to the particles that take part in the decay/inverse decay or annihilation processes contributing to the DM abundance production.
Recently, the Ref. [4,5] studied the FIMP DM scenario affected by the thermal masses and the phase transition effects. In this work, we present a novel multi-step FIMP production mechanism of the DM abundance after taking into account the SFOEWPT after the reheating of the Universe. In the model to be computed, a pseudo-Dirac sterile neutrino is introduced. In the literature, such kind of sterile neutrinos are utilized to give rise to the neutrino masses through the linear or inverse seesaw mechanisms, with the linear seesaw can be differentiated from the inverse seesaw with the feasibility of leptogenesis [6] and lepton number violation search at colliders. Here, we note that the behaviors of the dark matter are not significantly affected by the types of the sterile neutrinos (Dirac or Majorana). See Ref. [7][8][9] for recent works on collider searches. In the dark sector, we introduce a hidden singlet scalar and a hidden fermion. Both these particles can serve as the DM particle depending on the mass spectrum [64]. For completeness, we study the complete set of Boltzmann equations including the thermal effects. The non-thermal production of the FIMP fermionic DM is an excellent benchmark where the gravitational wave signals of SFOEWPT can be reached by the projected GW detectors. The mixing between the active neutrino and the sterile neutrinos may be beyond the colliders search sensitivity. While, the SFOEWPT signals can be reached at LHC, see Ref. [10].

II. THE MODEL
In this paper, we utilize the model similar to Ref. [6,[11][12][13]. 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 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.
 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 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.

III. 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, 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 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 (h, φ) = which constrains the parameters as, The electroweak phase transition can be studied in a gauge invariant approach [14] after taking into account the thermal corrections, the thermal potential used to estimate the vacuum structures at finite temperature is given by, where 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. 5, Eq. 6, and Eq. 9 together bounds the parameter spaces of λ φ , λ φh and m φ , as shown in Fig. 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, FIG. 1: Two-step EWPT select conditions allowed parameter spaces. The gray region is excluded by the Higgs invisible decay bounds.
In the Electroweak vacuum, the thermal corrected mass are given by, for Goldstones, and . 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 Sec. IV. 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, 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 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, This condition can be converted to S 3 (T n )/T n = 4 ln(T n /H n ) ≈ 140 − 150 [15]. Where 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 Ref. [16][17][18]. 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. [19].
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, 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 Fig. 9. With the entropy difference between the high temperature and low temperature entropy: The Fig 9 indicates one obtain a negligible dilution factor as in the SM case, see Ref. [20], and the two-step phase transition pattern being studied in Ref. [21].

IV. 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 Ref. [22,23]). 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 firststage 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.
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 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. [24], 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 m N > m h (T ), 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 [25], with the kinematic function 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 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 Sec. III 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 (41), 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 (41). 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 Eq. (32, 41), we use the ready-made function [26,27] embeded in the micrOMEGAs [28] for computing the stiff equations Eqn. (32,41). With the model file being prepared by FeynRules [29], after implemented the thermal masses and VEVs of φ and h as a function of temperature as in Sec. III, we calculate all the thermal σv (s) and the thermal decay widths of the particles using the CalcHEP [30] 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 Sec. III. The degree of freedoms g * and g * S implemented in the micrOMEGAs are adopted in order to calculate M pl , and entropy s = 2π 2 45 g * S T 3 . Here the Planck scale is M pl = 1.22 × 10 19 GeV.

A. 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. 33 and Eq. 34. Which results in the bounds on the y χ as shown in Fig. 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.

B. 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 Fig. 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 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. [31] 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 Fig. 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 [32,33], LUX [34] and Panda X-II [35] for our interesting DM mass range, see also Ref. [36]. 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. • When z 1, as the temperature is extremely high, v φ (T ) = 0 and φ acquires a significant thermal mass through the Eqn. (8). 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 ) < • 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 Fig. 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 there becomes a second dip at the left panels in the Fig. 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 Fig. 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 onshell 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 Fig. 7 appears due to the first-order 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.  Fig. 7)), we need a larger coupling constant y D compared with the Fig. 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 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.  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 Fig. 8), such an influence can be highly reduced by the small φ abundance. In the Fig. 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 Fig. 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 Fig. 8).

V. 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.

A. 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 )), 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 [37,38]. The total energy spectrum of the three sources is given by, The first source of the gravitational waves from the bubble collision estimated using the envelop approximation [39][40][41] is [42], 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 α [43], v w 1/ √ 3 + α 2 + 2α/3 1 + α , κ 0.715α + 4 27 3α/2 1 + 0.715α , and the peak frequency f env is, Hz .
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 ; (54) the Hubble parameter at present is given by In Fig. 9, we show the possibility to probe the benchmarks of the Fig. 7 and Fig. 8 with the GW detectors. The sensitivities for the proposed space based interferometers: LISA [44] with two design configurations in notation NiAjMkLl [37,45], BBO, DECIGO (Ultimate-DECIGO) [46] and ALIA [47] are shown with different color shaded regions. For the GW signals produced by the SFOEWPT around the benchmark of the Fig. 7 is show in the top panels of Fig. 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 Fig. 9, we plot the GW signals from the SFOEWPT around the benchmarks of Fig. 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 Ref. [48][49][50].

B. 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 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 Fig. 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.

C. Collider interaply
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 [51][52][53], 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, We show the value of δm 2 h /Λ 2 in the panel of N φ and λ hφ , see Fig. 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 [54]. 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 [10]. For the search of this scenario at future linear colliders including the CEPC, ILC, and Fcc-ee we refer to Ref. [54,55]. The study of Ref. [55] 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.
For the typical mass of the RHN being adopted in this work, we have the decay length being given by [56], 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. [56]. The proposed SHiP experiment is also unable to probe the scenario due to its proposed mass range [57].

VI. 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 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. pseudo-Dirac 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. (1). m Dij = y 3×3 where y 3×3 N (C) are also the corresponding 3×3 extensions. Seesaw models require that µ 3×3 I,1,2 elements, the mixings between the light and sterile neutrinos are not sensitive to many of them except 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.
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 (B2). 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 a pair of pseudo-Dirac neutrinos for the simplicity of the calculations. The rotate matrix is perturbatively calculated to be up to the first order of v φ (T ). Here χ w , N 1,2 are the new rotated Dirac spinors, and these can still be combined into a Majorana and a Dirac spinor Then, for the small mixing limit, analogy to Ref. [5], we have new decay channel of and the dominant annihilation is the φ mediate s-channel process, with cross section being given by, σ(H † H → χχ) = 4(y χ sin θ) 2 (λ hφ v φ (T )) 2 (s − 4m 2 χ ) 3/2 8πs(m 2 φ (T ) − s) 2 s − 4m 2 h (T ) .

(B7)
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.