A two-component vector WIMP -- fermion FIMP dark matter model with an extended seesaw mechanism

We consider an extension of the Standard Model that explains the neutrino masses and has a rich dark matter phenomenology. The model has two dark matter candidates, a vector WIMP and a fermion FIMP, and the sum of their relic densities matches the total dark matter abundance. We extensively study the dark matter production mechanisms and its connection with the neutrino sector, together with various bounds from present and future experiments. The extra scalar field in the model may induce a first-order phase transition in the early Universe. We study the production of stochastic gravitational waves associated with the first-order phase transition. We show that the phase transition can be strong, and thus the model may satisfy one of the necessary conditions for a successful electroweak baryogenesis. Detectability of the phase transition-associated gravitational waves is also discussed.


Introduction
The Standard Model (SM) of particle physics has proved extremely successful in the past decades with the experiments matching its predictions and the Higgs boson discovery being the final piece to complete it. Nonetheless, astrophysical and cosmological evidence have posed questions that are not explained by the SM and are still open problems to this date. It is well established by neutrino oscillation data (NOD) that the neutrinos have a non-zero mass while in the SM they are massless. A mechanism is therefore needed to symmetry breaking, the dimension-5 operators also give IR contributions, dominated by the Higgs decay. The dominant contribution will depend on values of the Higgs mass, the reheating temperature T R , and the scale of new physics Λ.
The visible sector described by the SM is composed of a complex arrangement of particles and gauge groups. Likewise, we could expect the similar complexity to arise in the dark sector. There is no experimental indication that the DM sector is composed of a single field. Since both the freeze-in and freeze-out mechanisms are viable production mechanisms, both the WIMP and FIMP could have been active in the early Universe, producing parts of the total DM relic density Ω Tot h 2 = 0.120 ± 0.001 as observed by the Planck experiment [5]. Although the simplest setup would be the case where there are two DM candidates both of which contribute to the total DM relic abundance, one may consider a more general multicomponent DM scenarios. Recent studies on the multi-component DM scenarios include Refs. .
In this paper, we consider a beyond the SM (BSM) scenario that addresses the aforementioned problems, exploring its viability and the possible experimental signatures. We introduce two sets of three-generation extra neutrinos N i L and S i L where the first two generations of neutrinos are used in the extended seesaw mechanism to explain the light neutrino masses while the third generation will be part of the dark sector. In the mass basis, through mixing, N 3 m and S 3 m become FIMP-type particles, and considering S 3 m to be the lighter one, it may become a viable DM candidate. We shall explore possible connections between the neutrino parameters and the DM relic density. The second DM candidate is the vector gauge boson W D associated with an extra dark U (1) D gauge symmetry. We study constraints from the lepton flavour violation (LFV) data to the mixing angles and to the DM production via the neutrinos sector.
The rest of the paper is organised as follows. In Section 2, we present the model under consideration in detail and explain the generation of the neutrino masses as well as the DM candidates. We then scrutinise the LFV bounds on the neutrino sectors in Section 3. In Section 4, we study the two-component DM scenarios and present the result together with various detection bounds. In Section 5, possibilities of having a FOPT are explored. We also discuss the detectability of the stochastic GW signals associated with the FOPTs. In doing so, we present three benchmark points (BPs) that realise the neutrino masses, the correct DM relic density, a strong FOPT, and a detectable GW signal at the same time. Finally, we discuss potential collider searches in Section 6 before we conclude in Section 7.

Model
We consider the following Lagrangian where L SM is the SM Lagrangian with the SM Higgs field φ h , L N is the Lagrangian associated with the additional singlet neutrinos which take part in the neutrino masses, and L DM corresponds to the DM Lagrangian. The fourth term describes the kinetic term for the extra U (1) D Higgs φ D , and the covariant derivative is given by D µ = ∂ µ − ig D W Dµ with g D being the U (1) D gauge coupling and W D the vector boson associated with the extra U (1) D gauge symmetry. F αβ D is the field strength of the vector boson W D , B αβ is the field strength tensor associated with the hypercharge U (1) Y gauge group, and the gauge kinetic mixing between these two field strength tensors is parametrised by ζ. Finally, the last term represents the scalar potential which is given by The neutrino sector is described by We have considered the Yukawa term for S i to be negligible compared to the one for N i , following the standard extended double seesaw model [8,9]. The Yukawa terms with the dark Higgs φ D are forbidden by symmetries. The parameters M ij S , M ij R , and µ ij are constants with mass-dimension one, while y ij are dimensionless coupling constants that compose the Dirac mass matrix M D that we will use later. We have considered that S 3 L and N 3 L are decoupled from the rest of the particle spectra by assuming that they are Z 2 -odd while the rest of the particles are Z 2 -even. Such a discrimination ensures that the lightest particle may be treated as a FIMP-type DM candidate. Productions of S 3 L and N 3 L are through dimension-5 operators which get naturally suppressed when the scale of new physics Λ is large, ensuring feeble interactions with the rest of the particle spectra; in the present work, we consider Λ ≥ 10 14 GeV. In the neutrino sector, the effect of such dimension-5 operators is negligible. The Lagrangian associated with S 3 L and N 3 L is thus given by Table 1. Particle contents and their corresponding charges under different gauge groups and discrete symmetry. The index i is for three flavours, running from 1 to 3 whereas the index j runs from 1 to 2.
Finally, the term proportional to the coupling ζ denotes the gauge kinetic mixing term between the SM U (1) Y gauge boson and the U (1) D gauge boson. It has been shown that small values of the parameter ζ are favoured from the viewpoint of the muon g−2 [122,123]; see also, e.g., Ref. [124] for various experimental constraints on ζ. In this work, we shall ignore the gauge kinetic mixing term to ensure that the U (1) D gauge boson W D becomes a stable WIMP DM candidate. One may alternatively impose an upper bound of ζ 10 −20 by requiring that the lifetime of W D is larger than the age of the Universe; see Appendix A for details. 2 Additionally, to consider W D as the WIMP DM, which is one of the main motivations of the present work, we consider all the particles to be neutral in U (1) D except the singlet scalar φ D which is necessary for obtaining the W D mass. Introduction of U (1) D charges to any other fields would make W D unstable. Table 1 summarises the particle contents of the model under consideration and their charges. The presence of an extra scalar that interact with the φ h induces a mixing between the two. In unitary gauge, the expressions of φ h and φ D , after the spontaneous breaking of the gauge symmetry, are given by with the mass matrix Here, v (v D ) denotes the VEV of the SM (dark) Higgs φ h (φ D ). Diagonalisation of the mass matrix leads to the mass eigenstates where θ is the mixing angle, given by (2.8) The mass eigenvalues are We consider the case where the dark Higgs is lighter than the SM Higgs. In other words, M H 2 denotes the mass of the dark Higgs, and M H 1 matches the SM Higgs mass. One may express the scalar quartic couplings in terms of the mixing angle and masses of the physical Higgses; we present the expressions in Appendix B. When the dark Higgs acquires a VEV, the U (1) D gauge boson W D gets the mass of M W D = g D v D .

Generation of neutrino masses with the extended seesaw mechanism
We consider the first two generations of the additional fermions, namely N 1 L , N 2 L , S 1 L , and S 2 L , to take part in the neutrino mass generation. The third generation is decoupled from the visible sector which is achieved by making them Z 2 -odd. Therefore, the neutrino mass matrix can be expressed as Here, M D is the 2 × 3 Dirac mass matrix, where m ij D = y ij v/ √ 2 and the superscript R (I) stands for the real (imaginary) part. On the other hand, M R and M S in Eq. (2.10) are 2 × 2 matrices which we choose to take as follows: (2.12) Finally, we choose µ as a symmetric matrix. It is in general complex, and we parametrise it as (2.13) To realise the non-zero neutrino masses in the extended seesaw framework, we consider the following hierarchy amongst the elements of M D , M R , and M S mass matrices [9]: (2.14) With these assumptions, we can diagonalise the mass matrix shown in Eq. (2.10) and obtain the following set of mass matrices [9]: Once we diagonalise the 3 × 3 matrix m ν , we get the masses of the active neutrinos. The other two matrices give the masses of the sterile neutrinos. After the diagonalisation, we find the relation between the flavour basis (ν L S L N L ) T and mass basis (ν m S m N m ) T as and where U , which is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix [126,127], W S , and W N diagonalise m ν , m S , and m N , respectively [9].

FIMP dark matter candidate
The remaining singlet neutrinos, N 3 L and S 3 L , comprise a 2 × 2 mass matrix, and the lighter one is a good DM candidate. The mass matrix for the DM sector takes the following form: where the elements are given by . (2.20) In the limit M 33 R M 33 SN , we diagonalise the mass matrix to obtain the eigenvalues expressed as The relation between the mass eigenstates and the flavour eigenstates are given by where δ ≡ M 33 SN /M 33 R . In our study, S 3 m is the lighter one, becoming a good DM candidate, and N 3 m is the next-to-lightest stable particle (NLSP). We shall drop the superscript '3' and use S m and N m to denote the lighter third-generation mass eigenstate, which becomes the FIMP DM, and the heavier third-generation mass eigenstate, which is the NLSP, respectively, hereinafter.
In Section 4, we shall explore the DM phenomenology in detail, focusing on the parameter space where Λ 10 14 GeV and α = α = ξ = ξ = κ = κ = O(1). In this case, both the DM candidate S m and the NLSP N m are produced out-of-equilibrium in the early Universe through the freeze-in processes of the type H i (+H i ) → FIMP + FIMP, where i = 1, 2. Thus, the effective couplings are in the ballpark of the FIMP-type DM. After the production, the NLSP decays into the lighter eigenstate; see Section 4 for details.

WIMP dark matter candidate
The vector gauge boson W D associated with the extra dark U (1) D is, on the other hand, a good WIMP DM candidate in our model. The vector boson W D , being a WIMP, is produced via the standard freeze-out processes SM + SM ↔ W D + W D which keep the WIMP in thermal equilibrium with the SM thermal bath.
As the model features both the FIMP DM, S m , and the WIMP DM, W D , a twocomponent DM scenario naturally arises in our model. We will present the detailed analysis in Section 4.

Neutrino masses and lepton flavour violation bounds
The eigenvalues of the mass matrix m ν represent the masses of the active neutrinos as we discussed in Section 2.1. Differences of their mass-squared will give us the solar mass difference ∆m 2 12 and the atmospheric mass difference ∆m 2 31 . On the other hand, the elements of the PMNS matrix U give us the oscillation angles θ 12 , θ 13 , and θ 23 . In this work, we consider the recent bounds on the oscillation parameters [ Additionally, we consider the bound on the sum of the non-decaying active neutrino masses from cosmology, i.e., i m ν i < 0.23 eV [4,5]. Since we have additional sterile neutrinos, we also take into account the LFV processes. The most stringent bounds on the LFV processes come from µ → eγ, µ → eee, and µ-to-e conversion CR(µ − Ti → e − Ti). The recent bounds are given by Br(µ → eγ) < 4.2 × 10 −13 [128], Br(µ → eēe) < 1 × 10 −12 [129], and CR(µ − Ti → e − Ti) < 6.1×10 −13 [130]. In determining the branching ratios for µ → eγ, µ → eēe, and µ − e conversion rate, we follow Refs. [131,132]. In order to satisfy all the aforementioned constraints, the elements shown in Eq. (2.11) and Eq. (2.12) cannot take arbitrary values. We thus vary the model parameters as below and obtain allowed parameter spaces by imposing the aforementioned constraints: where α = {e, µ, τ } and i, j = 1, 2, and we have chosen The rest of the model parameters, which affect the DM relic density directly but do not take part in the neutrino mass, are fixed as These fixed values are inspired by the DM studies as well as the FOPTs, as we will discuss later in the paper.
In the left panel (LP) and right panel (RP) of Fig. 1, allowed parameter spaces are shown in the m e1 D -m µ1 D and µ R 11 -m e1 D planes, respectively. The cyan points are obtained after imposing the NOD constraints. The blue points are obtained when we additionally impose the LFV bounds. From the LP of Fig. 1, one may see a sharp correlation between m e1 D and m µ1 D . This is because both of them actively contribute to the neutrino mass, i.e., they are the leading contributions in two different elements of the neutrino mass matrix m ν . Since we have taken the elements of M S to be equal, for a large value of M S , we need a small value for m e1 D and m µ1 D , and similarly, for a small value of M S , we need a large value for m e1 D and m µ1 D . Moreover, the LFV processes are mediated by the gauge bosons (W ± , Z), so those processes mainly depend on the active-sterile mixing terms, namely M D /M S and M D /M R . Therefore, when we apply the LFV bounds, higher values of m µ1 D get ruled out for each value of m e1 D . On the other hand, in the RP of Fig. 1, we see an anti-correlation between m e1 D and µ R 11 , which is mainly due to the neutrino mass relation. Furthermore, elements of the matrix µ do not actively contribute to the LFV processes. Therefore, there is practically no shrink in the µ R 11 -m e1 D plane after applying the LFV bounds. In   of the active-sterile mixing. Therefore, in the LP of Fig. 2, we see that both parameters cannot take higher values simultaneously due to the LFV bounds. The same conclusion is also observed for the RP. Depending on the strength of the active-sterile mixing, we may detect the sterile neutrinos in many ongoing and future experiments which we shall discuss later in Fig. 4 and Fig. 5. In Fig. 3, we have shown the scatter plots in the M D /M S -µ plane after imposing the NOD and LFV bounds. In the LP of Fig. 3, blue, magenta, green, and cyan points respectively correspond to m e1 In the RP, we have blue, magenta, green, and cyan points for , respectively. One interesting point to note here is that there is a strong correlation amongst the blue, magenta, and green points in both the LP and RP, while we observe no relation amongst the cyan points. The points that exhibit the strong correlation strictly follow the relation (M D /M S ) 2 µ < 10 −11 GeV, which is the mass of the active neutrinos. Parameters denoted by the cyan points do not affect the neutrino mass directly; they either come with the multiplication of other terms or are absent in the neutrino mass matrix. Thus, in the end, their combinational effect never exceeds the light active neutrino mass.
In the LP of Fig. 4, we have shown the allowed parameter space in terms of the Yukawa coupling y e1 and the DM relic density that is coming solely from the neutrino sector. S m , which is a FIMP DM candidate as we discussed in Sec. 2, may be produced via annihilations of active neutrinos and extra heavy neutrinos, mediated by the Higgses, as Fig. 4 indicates that this contribution is subdominant. One may understand the general behaviour as follows. When M N 1 is smaller than 500 GeV, we have a linear relation between y e1 and the DM relic density coming from the active and heavy neutrinos annihilations. It reflects the fact that Ω ν−Cont.
DM h 2 ∝ y 2 e1 . When M N 1 is larger than 1000 GeV, the contribution to the DM relic density is small as the mass is close to the chosen reheating temperature of T R = 3 TeV; thus, a suppression occurs. We observe that, for the chosen range of DM h 2 (left) and m S1 -|V eS1 | 2 (right) planes after imposing the NOD constraints. Bounds coming from various ongoing (solid lines) and future (dashed lines) experiments are overlaid; see text for more details. Figure 5. Allowed parameter spaces in the m S1 -|V µS1 | 2 (left) and m S1 -|V τ S1 | 2 (right) planes after imposing the NOD constraints. Bounds coming from various ongoing (solid lines) and future (dashed lines) experiments are overlaid; see text for more details. parameter values (3.2), the contribution of the active and extra heavy neutrinos to the total DM relic density is at most ∼ 3%. The RP of Fig. 4 depicts the allowed region in the active-sterile mixing angle associated with electron |V eS 1 | 2 and sterile neutrino mass m S 1 plane after imposing the NOD bounds. The solid lines represent the present bounds which come from CHARM [133,134], BELLE [135], and DelPhi [136], depending on the mass of the sterile neutrino. DelPhi demands the allowed range |V eS 1 | 2 < 10 −4 for the sterile neutrino mass up to 100 GeV, whereas CHARM puts a bound on the active-sterile mixing angle |V eS 1 | 2 10 −7 for the sterile neutrino mass less than 2 GeV. There are various proposed experiments, including DUNE [137,138], SHiP [139], MATHUSLA [140], LNV-Meson [141], FCC-ee [142,143], CMS [144], and LHCb [144,145], which have the sensitivity reaching up to |V eS 1 | 2 ∼ 10 −10 for the sterile neutrino mass up to 100 GeV. Figure 5 shows the allowed region in the active-sterile mixing associated with the muon |V µS 1 | 2 (LP) as well as tauon |V τ S 1 | 2 (RP) and the sterile neutrino mass m S 1 planes, after imposing the NOD bounds. In the LP, the recent bounds put by NuTeV [146] and FMMF [147] already rule out the sterile neutrino mass up to 2 GeV for the active-sterile neutrino mixing larger than 10 −7 . Various future experiments such as DUNE [137,138], SHiP [139], MATH-USLA [140], FASER2 [148], LNV-Meson [141], AL3X [149], FCC-ee [142,143], CMS [144], and LHCb [144,145] are also presented by dashed lines which will probe the active-sterile mixing, |V µS 1 | 2 , up to 10 −10 for sterile neutrino mass as large as 100 GeV. On the other hand, from the RP of Fig. 5, we see that the DELPHI experiment [136] already rules out |V τ S 1 | 2 > 3×10 −5 for the sterile neutrino mass up to 100 GeV. The future experiments shall cover the active-sterile mixing up to |V τ S 1 | 2 ∼ 10 −10 for the mass range up to 100 GeV. For all the active-sterile mixing and sterile neutrino mass planes, there exits a bound coming from the Big Bang Nucleosynthesis (BBN) as well, if the sterile neutrino decays after the BBN. However, the BBN bound is weak for the parameter space we have considered.

Dark matter phenomenology
We now discuss the production and detection prospects of the DM candidates in our model. 3 Our model features both the WIMP and FIMP DM candidates as we discussed in Section 2. The dark gauge boson W D plays the WIMP role, and the lighter singlet neutrino of the third generation S m becomes the FIMP DM; the NLSP, N m , will eventually decay to the FIMP DM S m . Thus, a two-component DM scenario naturally arises. The WIMP part ensures the potential detectability in future, whereas the FIMP DM will be difficult to probe by the direct, indirect, or collider detection techniques.
As we shall discuss in Section 5, a low-mass BSM dark Higgs is favoured from the FOPT point of view [153]. Thus, in this section, we mainly focus on the range 1 − 200 GeV for the dark Higgs. Furthermore, to avoid any potential problems with collider searches due to the low mass of the dark Higgs, we consider the mixing angle θ in Eq. (2.8) to be small, focusing on | sin θ| < 0.1. In doing so, we may easily evade the Higgs signal strength bounds [154,155].
There are mainly five constraints that we have taken into account for the discussion of DM phenomenology: i) relic density, ii) direct detection bounds, iii) indirect detection bounds, iv) Higgs invisible decay, and v) Higgs signal strength bound. We explain each category below before presenting the results. • DM relic density: We consider the bound on the total amount of DM relic density coming from the Planck experiment [4,5]. Specifically, the following 3σ bound is used, unless stated otherwise: Here, Ω W D h 2 (Ω Sm h 2 ) denotes the WIMP (FIMP) DM relic density.
• Direct detection: In our model, DM can have the elastic scattering with a nucleon N as depicted in Fig. 6. The analytical estimate for such a process takes the form [156], is the reduced mass, with M N being the nucleon mass, v is the electroweak VEV, Z is the atomic number, A is the atomic weight, and f α (α = p, n) can be expressed as .026(0.020), and f p,n Ts = 0.043 [157]. As we take the WIMP DM mass to be in the range of 1 − 100 GeV, the DM may be detected by different experiments. A part of the parameter space in the spin-independent direct detection (SIDD) cross-section, σ SI , and DM mass, M W D , plane is already ruled out by LUX-ZEPLIN-5.5T [158], PandaX-4T [159,160], and Xenon-1T [19] for the 10 − 100 GeV DM mass range. On the other hand, the mass range below 10 GeV will be explored by experiments such as DarkSide-50 [161,162], XENON-1T (M) [163], CDMSlite [164], and CRESST-III [165,166]. Our SIDD cross-section is a few orders of magnitude below the current bound.
• Indirect Detection: The WIMP DM can also be detected by observing the annihilation products, namely bb, ττ , µμ, and eē. When the WIMP DM mass is above the b-quark mass, the bound from the bb final state dominates. Fermi-LAT + MAGIC Segue 1 [21] puts the stringent bound on the σv bb -M W D plane. On the other hand, when the DM mass is smaller than the b-quark mass, the DM annihilates to ττ , µμ, and eē dominantly. The bounds come from the study of FERMI-LAT [167,168], CMB [168], and AMS [168,169]. We shall discuss the details of the indirect detection bound when we present our resultant plots.
• Invisible decay: When the DM mass is below half of the SM Higgs mass, there is a possibility that the SM Higgs will have an invisible decay, Γ inv H 1 . Thus, one needs to make sure that the invisible decay is always smaller than the current bound [170], The decay width of the SM Higgs to the WIMP DM in the present case takes the form, The allowed parameter region in the M W D -σ SI plane after imposing the invisible decay constraints shall be presented in Section 4.2 with the SM Higgs decay width Γ H 1 = 4.156 MeV.
• Higgs signal strength: The Higgs signal strength can be estimated by measuring its production and decay ratio with the SM values. It can be defined as where µ H 1 = σ H 1 /σ SM H 1 is the ratio of the Higgs production in the new model and the SM, and µ f = B f /B SM f is the ratio of the branchings of the Higgs to a channel f . The current bound onμ after a combined analysis is given by [154] µ = 1.17 ± 0.10 . (4.7) Assuming that the Higgs boson has the same kind of branchings as the SM case, we find thatμ ∼ cos 2 θ. By taking the 3σ range, we obtain that sin θ < 0.36. Since we consider a small mixing angle, namely sin θ < 0.1, we thus always satisfy the bound from the Higgs signal strength.

Dark matter production
Let us temporarily consider a regime where only the coupling κ is active and the mixing between S 3 L and N 3 L is negligible, in which case S m S 3 L . Let us also consider the case where the mixing between the Higgses is small, i.e., cos θ 1. In this case, the FIMP DM S m is produced dominantly by the Higgs scattering process, and the analytical solution for the yield is given by [28] where T R is the reheating temperature, T end 1 MeV is the temperature after which we may safely assume that no DM production occurs, and the entropy S and the Hubble parameter H are given by with M Pl = 1.22 × 10 19 GeV, and g S and g ρ being respectively the entropy and the energy density degrees of freedom of the Universe; we take g S,ρ ∼ 100. In achieving Eq. (4.8), it is assumed that the masses of the associated particles in the production may be neglected compared to the temperature at which DM production happens which, for the process under consideration, is the reheating temperature. Thus, the production depends on the highest temperature, obtaining the UV freeze-in contribution. Consideration of masses of the associated particles has therefore a negligible effect in the DM production. Moreover, as in the IR freeze-in case, it is assumed that one may safely ignore the back-reaction of the DM in the Boltzmann equations since the number density is always smaller than the equilibrium number density. With these assumptions, the final relic density for the T R -dominated regime is given by where S 0 /ρ c 2.74×10 8 is the ratio of the entropy today and the critical energy density. We have checked that the result from numerical analyses performed by using micrOMEGAs [152] for the T R -dominated regime matches well with the analytical expression (4.8). Figure 7 shows the relic density of the FIMP DM, Ω Sm h 2 , in the T R − Λ plane, using the analytical result. The black dashed line indicates the correct relic abundance. We see that low values of the reheating temperature are preferred. For the rest of the work, we will therefore concentrate on the low reheating temperature. In particular, we shall choose T R = 3 TeV throughout this section.
With the knowledge obtained above, we now re-introduce all the couplings and numerically evolve the full Boltzmann equations using micrOMEGAs to obtain the DM relic densities. The relevant Boltzmann equations are where z ≡ M H 1 /T , σv AB→CD is the thermally-averaged cross section times velocity, Γ A→BCD is the thermally-averaged decay width, and θ(x) is the Heaviside step function. We present the relevant Feynman diagrams in Appendix C. In Fig. 8, the DM production by freeze-out and freeze-in mechanisms are shown. The green double-dot-dashed line corresponds to the WIMP DM production by the freeze-out mechanism. It freezes out at T M W D /20 which corresponds to z 2500. The cyan dashed line represents the production of the NLSP N m , which later decays to the FIMP DM S m at z 3500. The NLSP is produced in the early Universe at T 3000 GeV, i.e., z 0.03 through 2 → 2 processes present in our model. The purple dot-dashed line indicates the FIMP DM production via the freeze-in mechanism. At its initial production, we see a sharp rise at z = 0.03 which represents the production by the 2 → 2 processes like the NLSP case. There exists a second rise in the production shortly after z = 1 which is due to the decay of the SM-like Higgs, H 1 . Finally, a third rise happens at z 3500 when the NLSP decays to the FIMP DM. The sum of the WIMP and FIMP DM relic densities is depicted by the black solid line, which coincides with the Planck measurement of total DM relic density Ω Tot h 2 = 0.12 today which is represented by the grey solid line. We have chosen the parameter values in such a way that the WIMP and FIMP DM contribute equally, namely Ω Sm h 2 ≈ Ω W D h 2 ≈ Ω Tot h 2 /2.
Dependence of the DM relic densities on the FIMP DM mass is shown in the LP of Fig. 9. One may see that the variation of the FIMP DM mass does not alter the WIMP DM relic density, which is depicted by double-dot-dashed lines. The dashed lines correspond to the NLSP (N m ) relic densities. The decay length of NLSP is not affected by the DM mass unless we choose M Nm M Sm + M H 1,2 . The dot-dashed lines below Ω DM h 2 10 −1 are the FIMP DM evolutions. We see that, for M Sm = 1 and 20 GeV, there is a slight rise in the DM density which corresponds to the FIMP DM production from the SM Higgs decay at around z = 1. This rise is, however, negligible for the M Sm = 50 GeV case due to the phase space suppression from the SM Higgs decay. The second rise at z ∼ 3500 happens when NLSP decays to the FIMP. Total DM relic density, which is the sum of the WIMP and FIMP relic densities, is represented by the solid lines. We see that the total DM relic density mainly follows the WIMP DM relic density. The RP of Fig. 9, shows the dependence of the DM relic densities on the NLSP mass. The NLSP masses are all above 100 GeV, so its production happens through 2 → 2 processes. The NLSP relic density varies linearly with its mass, and its contribution to the DM relic density is given by (M Nm /M Sm )Ω Nm h 2 . This is similar to the SuperWIMP mechanism [171] and associated with the conservation of the comoving number densities between two out-of-equilibrium species. For example, for a process A → B + · · · , if A and B are out of equilibrium, then one has Y A = Y B , and thus,  their production strength is inversely proportional to Λ. We see that the green lines, which correspond to Λ = 5.5 × 10 13 GeV, have larger relic densities than the blue and pink lines, which respectively correspond to Λ = 5.5 × 10 14 GeV and Λ = 5.5 × 10 15 GeV. The RP of Fig. 10 shows the variation of the DM relic densities for three different values of the U (1) D gauge coupling g D . The green lines are for g D = 10 −3 , the blue lines are for g D = 3.1×10 −4 , and the pink lines are for g D = 10 −5 . We see that the g D = 10 −5 case, which is depicted by the pink lines, has the largest WIMP DM relic density. It is because a smaller value of g D reduces the WIMP DM annihilation cross-section (W D + W D → SM + SM) which affects inversely the WIMP DM relic density. On the other hand, the FIMP DM relic density is controlled by the strength of the dark Higgs VEV v D = M W D /g D . The VEV v D is linearly proportional to the coupling strength responsible for the FIMP DM production through the BSM Higgs decay. One interesting thing we may note is that g D does not affect the NLSP production as its production is governed by the 2 → 2 processes. However, the NLSP decay depends on the VEV v D ; the NLSP decays faster with the increment of v D or decrement of g D . Moreover, FIMP DM production from the decays of the SM Higgs has an effect only when v D is large enough; otherwise, the v D -associated part in H 1 → S m + S m production is suppressed.
In Fig. 11, we show the dependence of the DM relic densities on three different values of M H 2 (LP) and M W D (RP). From the LP, we see that there is no effect of M H 2 on the FIMP DM production, while the effect on the WIMP DM production is significant. It is the case since the WIMP DM relic density is mainly controlled by how far we are from the resonance region of the second Higgs H 2 . On the other hand, from the RP of Fig. 11, we see that changing the WIMP DM mass affects both the WIMP and FIMP DM productions. The effect on the WIMP DM is due to the fact that, with the change of M W D , we are moving away from the resonance region of the second Higgs H 2 , and thus we have more production of the WIMP DM. Additionally, the NLSP decay is proportional to the VEV v D = M W D /g D . Therefore, by increasing the value of M W D , the value of v D increases as well, which triggers an early decay of NLSP.

Exploration of allowed parameter spaces
With the understandings we have acquired in Section 4.1 by studying the behaviours of the DM relic densities near the point (3.3), we attempt to obtain allowed parameter regions amongst the different parameters after imposing that the DM relic density satisfies the range 0.01 ≤ Ω DM h 2 ≤ 0.12. The lower limit of 0.01 is to ensure more allowed points. Note, however, that our conclusion remains unchanged even if the 3σ range shown in Eq. (4.1) is considered. We also discuss bounds on the WIMP DM parameters coming from the direct and indirect detections of DM. We perform parameter scans with the following parameter ranges: with the rest of the model parameters being fixed as those given in Eq.  due to the observation that the WIMP DM mass needs to be close to the resonance region, as shown in Fig. 11. The LP of Fig. 12 shows the allowed region in the M Sm -Λ plane where the colour represents the FIMP DM relic density. One may easily see that, for a fixed value of the FIMP DM mass, increasing the value of Λ makes the FIMP contribution to DM relic density decrease, as the FIMP DM production is inversely proportional to Λ. The lower limit in the Λ value comes from the maximum allowed range for DM relic density, Ω DM h 2 = 0.12, since the relic density is proportional to M Sm /Λ 2 . In the RP of Fig. 12, we present the allowed range in the M W D /g D -Λ plane. One may again observe that, as we go to a higher value of Λ, we have a smaller FIMP DM contribution. The VEV v D = M W D /g D linearly contributes to the FIMP DM relic density, and thus, for a higher value of v D , we need a higher value of Λ to get the correct DM relic density value; we notice this in particular in the region g D < 10 −3 and 10 < M W D [GeV] < 100. This correlation between v D and Λ is observed for higher values of v D , whereas we do not see such a correlation for lower values of v D . In both the LP and RP of Fig. 12, we see that there is no upper bound on Λ. This is because such higher values of Λ reduce the FIMP DM relic density, and the DM relic density bound can be satisfied from the contribution of the WIMP DM.
The LP of Fig. 13 shows the allowed range in the M W D -M H 2 /M W D plane after imposing 0.01 ≤ Ω DM h 2 ≤ 0.12. It is clearly shown in the figure that, to obtain the WIMP DM relic density below 0.12, we need to stay near the resonance region, i.e., M H 2 ∼ 2M W D . It is also clear that, when we are very close to the resonance region, we have a smaller WIMP contribution in the DM relic density, while a larger WIMP contribution is obtained as we depart from the resonance region. Moreover, we see that, for M W D 62.5 GeV, the dark Higgs mass M H 2 may take any value. This is due to the fact that the dominating contribution comes from the SM Higgs resonance. The RP of Fig. 13 shows the allowed How close one needs to be to the resonance region is studied in Fig. 14 where the WIMP DM relic density is shown in terms of the mixing angle sin θ and r ≡ 2M W D /M H 2 . The parameter r quantifies the closeness to the resonance region, and r = 1 corresponds to the exact resonance point. We find that r typically takes a value between 0.92 and 0.98, depending on the value of the mixing angle, if we ask for the WIMP DM component to be a significant part of the total relic density. We observe that the window for a WIMP DM relic density of at least 10% of the total DM relic density is narrower for smaller values of the mixing angle. At a fixed value of r, the relic density decreases as the value of sin θ increases. One may understand this as follows: The process keeping the WIMP DM in thermal equilibrium is DM + DM ↔ SM + SM, and it is mediated by H 2 . We thus find that the cross section is proportional to cos 2 θ sin 2 θ. Higher values of r mean being closer to the resonant point where the cross section increases. Therefore, the relic density decreases following the standard behaviour, Ω WIMP h 2 ∼ 1/ σv . Once we depart too much from the resonance region, we may overproduce the WIMP DM and overclose the Universe. Figure 15 shows the allowed region in the M W D -(Ω W D /Ω Tot )σ SI (LP) and M W D -(Ω W D /Ω Tot ) σv bb (RP) planes, together with various direct and indirect detection bounds that are depicted by solid lines. Note that we have rescaled the y-axes by the amount of the WIMP DM relic density compared to the total DM in the Universe Ω Tot h 2 = 0.12. The LP of Fig. 15 may be easily understood with the direct detection expression given by  Eq. (4.2), which states that σ SI is proportional to g 2 D . One may estimate the percentage of the WIMP DM relic that each sample point corresponds to with the help of the RP of Fig. 13. Comparing the LP of Fig. 15 and RP of Fig. 13, one can easily see that lower values of g D correspond to lower values of σ SI (∝ g 2 D ) and higher values of the WIMP DM relic density as the density is inversely proportional to g 2 D . A sharp dip at M W D 62.5 GeV happens because of the mutual cancellation between the SM Higgs-and the BSM Higgs-mediated processes as one may see from Eq. 4.2. A part of the M W D > 7 GeV region is already ruled out by the different direct detection experiments such as XENON-1T [19], PandaX-4T [159,160], and LUX-ZEPLIN-5.5T [158]. The region of DM mass below 7 GeV will be explored by DarkSide-50 [161,162], XENON-1T(M) [163], CDMSlite [164], and CRESST-III [165,166]. The black solid line corresponds to the bound from the Higgs invisible decay which is obtained by staying near the dark Higgs resonance region, i.e., M H 2 ∼ 2M W D , so that the WIMP DM never becomes over-abundant. The region above the black solid line is already ruled out by the current bound on the branching of the Higgs invisible decay mode. We note that our model predicts much lower values for σ SI compared to the aforementioned bounds. From the RP of Fig. 15, we see that there is a dip in σv bb for the WIMP DM mass below 5 GeV. This is due to the fact that, for this range, the channel W D W D → bb is not active. The region of M W D 10 GeV is constrained by the Fermi-LAT + MAGIC Segue 1 data [21]. We observe that most of the parameter space which contributes dominantly to the DM relic is already ruled out by the indirect detection bound. We have also checked the present bounds on the DM annihilation to µ + µ − and τ + τ − , and we present the results in Appendix D.

First-order phase transitions and associated gravitational waves
The extra dark U (1) D Higgs field not only gives a mass to the WIMP DM W D , but it also changes the vacuum evolution. We study the evolution of the vacuum state and the dynamics of the phase transition in this section. We first compute the one-loop finitetemperature effective potential, where V (0) is the tree-level scalar potential, V CW is the one-loop Coleman-Weinberg potential, and V (1) T is the finite-temperature correction. In terms of the background fieldsH and H D of the SM Higgs and the U (1) D dark Higgs, the tree-level scalar potential is given by The Coleman-Weinberg potential can generically be written as where the + (−) sign is for bosons (fermions), n i is the number of degrees of freedom of the species i, M i is the field-dependent mass, the constants c i are 1/2 for transverse gauge bosons and 2/3 for the rest, andμ is the renormalisation scale. The expressions for the field-dependent masses M i and n i are summarised in Appendix E. Depending on the choice of the renormalisation scaleμ, the effective potential V eff changes, and hence one may arrive at different results. This renormalisation scale dependence has been explored in e.g. Refs. [172,173] 4 . Together with the gauge dependence issue, we do not attempt to address the issue of the renormalisation scale dependence as it goes beyond the scope of the current work. We thus ignore the one-loop Coleman-Weinberg corrections in the followings by assuming that the renormalisation scaleμ is chosen in such a way that the Coleman-Weinberg corrections are minimised. The finite-temperature correction is given by [179] V (1) where I + is for bosons and I − is for fermions. The re-summed ring diagrams are taken into account by replacing the field-dependent masses as [180] where Π i (T ) are the thermal masses [181]. We present the thermal mass expressions in Appendix F. Up to the leading O(T ) order, the effective potential is thus given by where E SM ≡ (2g 3 2 + g 2 1 + g 2 2 3 )/(32π) and '· · · ' include sub-leading, negligible terms. Note that we have assumed that the dark gauge coupling g D is small enough not to affect the leading-order terms.
In the presence of the extra Higgs field, FOPTs may arise. FOPTs with a dark U (1) D have been studied in e.g. Refs. [112,[182][183][184]. For studies of the phase transition with an extra scalar field, see, e.g., Refs. [153,172,178,[185][186][187][188][189]. In particular, in Ref. [153] where the studied scalar potential has the same form as Eq. (5.7), it was shown both analytically and numerically that the FOPTs could be strong, characterised by v c /T c 1. Here, T c is the critical temperature at which the potential minima become degenerate, and v c is the VEV of the SM Higgs field at T c . It indicates that one of the Sakharov conditions for successful electroweak baryogenesis can be fulfilled [190]. Since the scalar field space is now two-dimensional due to the extra Higgs field, one may achieve either one-step or two-step phase transitions. Noting that the VEV of the U (1) D Higgs is non-zero at zero , the second step breaks the electroweak symmetry, giving [153] v One may see that strongly FOPTs, v c /T c 1, can be achieved when the dark Higgs is lighter than the SM Higgs. The one-step phase transition shows a similar behaviour [153]. Utilising the publicly available tool CosmoTransitions [191], we numerically compute v c /T c for a wide range of the parameter space and present the result in Fig. 16. The explored parameter range is as follows: We note that these four input parameters are the only relevant model parameters. The other model parameters can be derived from the above input parameters. The x-axis of Fig. 16 is defined as λ m ≡ λ h − λ 2 hD /(4λ D ). We observe that strong FOPTs could be achieved for small values of λ m , which is in good agreement with both the analytical estimate (5.8) and the results of Ref. [153].
FOPTs may produce observable stochastic GWs [95]. There are three main contributions to the GWs from the FOPT: bubble wall collisions Ω col h 2 , sound wave in plasma Ω sw h 2 , and the magneto-hybrodynamic turbulence Ω turb h 2 . The total GWs are then Ω GW h 2 Ω col h 2 + Ω sw h 2 + Ω turb h 2 . The GWs coming from the bubble wall collisions are given by [108] Ω col h 2 = 1.67 while the turbulence contribution is [108] Ω turb h 2 = 3.35 × 10 −4 H * β where h * = 1.65 × 10 −5 Hz T * 100GeV g * 100 1 6 . (5.12) Finally, the sound-wave contribution to the GW signal can be expressed as [188,[192][193][194] Ω sw h 2 = 4.80 × 10 −6 min 1, (5.13) For the expressions for f col , f sw , f turb , v w , κ, κ v , and κ trub , see Appendix G. Three parameters that play the key roles in the GW signal are where S E = S 3 /T is the Euclidean action of a bubble with S 3 being the three-dimensional action, ρ vac the released energy density during the phase transition, and ρ * rad = g * π 2 T 4 * /30, with g * being the number of effective degrees of freedom at T = T * . We take T * to be the nucleation temperature T n . We employ CosmoTransitions [191] to numerically compute the three key parameters, α, β/H * , and the nucleation temperature T n . In Fig. 17, we present the FOPT-associated GW signals for three BPs together with the sensitivity curves of future space-based GW experiments such as LISA, DECIGO, and BBO. The three BPs, that account for not only the neutrino masses and the correct DM relic density, but also the strong FOPTs, are summarised in Table 2. One may notice that the three BPs have different DM compositions. In the case of the first BP (BP1), both the WIMP and FIMP contribute equally to the total DM relic density, while the BP2 (BP3) is mostly composed of the FIMP (WIMP) DM. We see from   [195] for the data for LISA, BBO, and DECIGO, and Ref. [196] for the Ultimate-DECIGO.
BPs  Table 2. Three BPs. The first four columns represent the input model parameters, the fifth, sixth, and the seventh columns are GW-related quantities, the eighth column shows the strength of the FOPT. The last two columns denote the WIMP and FIMP contributions to the total DM relic density Ω Tot h 2 = 0.12; for the first BP, both the WIMP and FIMP equally contribute to the total DM relic density, while the second (third) BP is mostly composed of FIMP (WIMP) DM. In all the three cases, the LFV bounds are satisfied, and the neutrino masses can successfully be generated. The GW signals corresponding to the three BPs are shown in Fig. 17.

Collider Searches
The present work deals with the WIMP and FIMP-type DMs. Due to the feeble interaction of the FIMP, it is difficult to probe it at collider experiments. We can, however, focus on general search strategies for BSM particles in the context of the present work. In particular, we may study the production of the second Higgs H 2 at the pp or e + e − colliders and look for its subsequent decay. Suitably adjusting the WIMP DM mass allows the second Higgs to decay mainly to the WIMP DM, and we may look for the missing energy with mono-jet or di-jet signals in the final state. Otherwise, if H 2 does not dominantly decay to the WIMP DM, then it will decay to the SM particles such as the SM Higgs. See, e.g., Refs. [45,197].
The relevant signal channels for our current study at the pp collider would be where j corresponds to the initial or final state jets, l is associated with the SM lepton, and E T is the transverse missing energy. Similar to the SM Higgs searches, we can also investigate e + e − → ZH 2 → nj + pl (n, p → 1) at the e + e − collider. The exact values of the integers, n, m, and p, depend on the production cross section and dominance of the associated backgrounds. Moreover, exploring the singlet fermions (S L , N L ) of our model at different colliders is an interesting direction; see, e.g., Ref. [197]. Further comments require a full-fledged collider study which is out of the scope of the current work, and we leave it for future study.

Conclusion
We have studied an extension of the Standard Model that accounts for the dark matter and the smallness of the neutrino masses under the extended seesaw framework. In our model, two sets of three-generation neutrinos are introduced; the first two generations provide the light neutrinos with a mass, and the third-generation neutrinos become FIMP-like particles. Amongst these FIMP-like particles, the heavier one eventually decays into the lighter one, and thus, we have the lighter third-generation neutrino as the FIMP dark matter candidate. Our model also contains a WIMP dark matter candidate, namely the dark U (1) D gauge boson. Thus, a two-component WIMP-FIMP dark matter scenario naturally arises in our model. We have explored allowed parameter spaces by using the lepton flavour violating bounds as well as the neutrino oscillation data. Much of the parameter spaces are already tightly constrained, but we have shown that there are viable parameter regions which are free from the constraints. Prospects of various future experiments have been discussed as well. Interestingly, the contribution to the FIMP dark matter relic density coming from neutrinos scattering is found to be up to a 3% of the total relic density for the range of the model parameters considered in our study. We have also discussed the dependence of the relic density on the model parameters. Utilising publicly available tools, we have performed extensive numerical parameter scans in order to study the evolutions of the dark matter candidates. Parameter spaces compatible with the bounds from (in-)direct detection and collider searches are presented. In particular, we have showed regions where a twocomponent dark matter scenario is realised and testable by future (in-)direct experiments.
The dark U (1) D Higgs field plays a major role in the FIMP and WIMP dark matter productions. In addition, the extra scalar field also changes the evolution of the vacuum state in the scalar sector, making a first-order phase transition possible. We have demonstrated that the strength of the electroweak first-order phase transition, quantified by the quantity v c /T c , where T c is the critical temperature and v c is the SM Higgs vacuum expectation value at T c , may become larger than unity for small values of the dark U (1) D Higgs mass. Therefore, one of the essential ingredients for a successful electroweak baryogenesis is achieved in our model. We have also studied stochastic gravitational waves associated with the first-order phase transitions and showed that the gravitational wave signals are strong enough to be detectable by future experiments such as BBO and DECIGO.
Three benchmark points, that explicitly demonstrate the capability of i) having a correct dark matter relic density, ii) generating the non-zero neutrino masses with the extended seesaw mechanism, iii) achieving a strongly first-order phase transition, and iv) emitting stochastic gravitational waves detectable by future experiments, are presented. Thus, the model studied in this work has an exciting potential detectability not only with future (in-)direct detection experiments and collider searches, but also with future gravitational wave experiments.
Acknowledgments J.K. would like to thank Yikun Wang for useful discussions on the phase transition and the use of CosmoTransitions. S.K. would like to acknowledge Geneviève Bélanger for the help related with micrOMEGAs. The work of F.C. is supported by the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 860881-HIDDeN. This work used the Scientific Compute Cluster at GWDG, the joint data center of Max Planck Society for the Advancement of Science (MPG) and University of Göttingen.

A WIMP DM decay width through kinetic mixing
In the presence of the mixing between the WIMP DM and the U (1) Y gauge boson, the DM may decay to, e.g., electrons, through the coupling between the DM and SM fermions [123]. For simplicity, we consider the decay of the DM to electrons. The decay width is then given by where g W D ee = 3eζ/(4 cos θ w ), e = √ 4πα, α is the fine-structure constant, θ w is the weak angle, and ζ is the gauge kinetic mixing parameter introduced in (2.1). Considering the DM mass of 10 GeV, and requiring the lifetime of the DM to be is larger than the age of the universe, we get an upper bound on the gauge kinetic mixing parameter as ζ < 10 −20 . When the decay of the DM to the SM fermions is open, the γ-ray observation may become relevant [125]. In this case, the DM lifetime should be greater than 10 29 s [125], which puts an even stronger bound of ζ < 10 −26 .  Figure 18. Feynman diagrams for LFV processes.

B Quartic couplings
The scalar quartic couplings may be written in terms of the mixing angle and masses of the physical Higgses as follows: C Feynman Diagrams Figure 18 shows the diagrams which contribute to the processes µ → eγ, µ → eeē, and µto-e conversion. We have considered these diagrams for the discussion of the LFV bounds. The Feynman diagrams relevant for our DM analysis are shown in Fig. 19.