Upgrading sterile neutrino dark matter to FImP using scale invariance

In this article we propose a class of extremely light feebly interacting massive particle, FImPs. They are combination of feebly interacting massive particle with scale invariance, by which DM stability, mass origin and relic density are inherently related. In the scale invariant version of the Standard Model (SM) with three right-handed neutrinos (ν\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu $$\end{document}SISM), the lightest N1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_1$$\end{document} realizes the FImP scenario. In this example scalar singlets, which are intrinsic to the ν\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu $$\end{document}SISM, generate mass and relic density for this FImP simultaneously. Moreover, they are badly needed for electroweak symmetry spontaneously breaking. Interestingly, a 7.1 keV N1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_1$$\end{document} with correct relic density, which can explain the recent 3.55 keV X-ray line, lies in the bulk parameter space of our model.


Introduction and motivation
The existence of dark matter (DM) is commonly believed in; nevertheless, its particle properties are still in the dark. The hypothesis that DM is a weakly interacting massive particle (WIMP) prevails, mainly due to the so-called WIMP miracle, which says that if DM annihilates involving weak scale mass and strength, it will obtain the correct order of magnitude for relic density DM h 2 ∼ 0.1. Despite a lot of experimental efforts made (including (in)direct DM detections and also collider searches), one still has not detected any trace of DM with confirmation. These null results [1] begin to challenge the WIMP paradigm. Now it is the right time to reexamine the theoretical basis of WIMP DM and explore other DM frameworks, or more concretely, models. We would like to follow the ensuing basic questions about DM.
• First, why is it there? A DM candidate that is predicted rather than introduced would be more attractive. A good case in point is the lightest sparticle in SUSY and the a e-mail: zhaofengkang@gmail.com axion from the Peccei-Quinn models aiming at solving the strong CP problem. Another example is the core of this paper, the lightest right-handed neutrino (RHN) N 1 , or sterile neutrino [2]. Originally RHNs are introduced to explain the nonzero neutrino masses, but they are neutral and thus have potential to provide a DM candidate. • Second, why is it stable (or at least sufficiently long-lived on a cosmic time scale)? Again, we advocate a mechanism that naturally guarantees DM stability rather than imposing a protective symmetry like Z 2 by hand. For instance, such a Z 2 may be an accidental symmetry due to the field content along with the space-time and gauge symmetries of the model, see an example based on scale invariance (SI) [3,4], 1 a symmetry which may address the hierarchy problem . Sterile neutrino DM offers another line for stability: lightness along with feebly interactions. No exact symmetry, by hand or accidental, is invoked here, and thus DM is expected to decay, say, into the very light SM species like neutrino and photons. But the decay rate is greatly suppressed by powers of DM mass and couplings, so DM can be sufficiently long-lived. • Last but not the least, why is the DM relic density DM h 2 around 0.1? The WIMP scenario presents the WIMP miracle. 2 We find that a similar numerical coincidence may also arise in the framework of feebly interacting massive particle (FIMP) [32,33]. Moreover, it is tied to the mass origin and stability of DM.
Let us neglect the first question for the time being. The answer to the second and third questions may point to a framework which combines FIMP with SI. The FIMP DM X (e.g., a Majorana fermion like N 1 ) feebly interacts with other fields, including the Higgs-like field S, which generates its mass via the term S X 2 , so it will be very light as long as the vacuum expectation value (VEV) of S is not hierarchically above TeV. Lightness plus feeble interactions are appropriate to give accidentally long-lived DM. Additionally, that mass source term can simultaneously freeze-in FIMP with correct relic density. Hereafter that kind of FIMP is dubbed FImP with "m" indicating its lightness due to SI. Needless to say, FImP cannot be seen by the current DM detectors, except for these sensitive to X-ray. We mention that, with regard to experimental tests, FImP, which tends to produce X-ray, is more interesting than FIMP, which usually stays in the dark.
The FImP example N 1 from the νSISM is of special interest, because it gives answers to all those three questions. Interestingly, there are some suggestive hints which favor such DM in the X-ray line at energy 3.55 keV, which was reported recently and can be explained by the 7.1 keV N 1 with activesterile neutrino mixing angle 0.8 × 10 −5 [34,35]. Despite the difficulty existing in the conventional ways, such a N 1 can easily get a correct relic density via freeze-in [36][37][38]. In this paper we will concentrate on the N 1 example, and another interesting example based on scalar FImP is considered in Appendix C.
In a later stage of this paper, we found that the idea of introducing a scalar singlet to freeze-in N 1 has already been explored by several groups [39][40][41][42]. Even then, there are several obvious differences between our paper and theirs. First of all, only in this paper the singlet is a built-in rather than artificial ingredient. Second, only here the singlet plays a crucial role in EWSB, while in others it is just a spectator. EWSB will be one of the centers of our article. Third, only here N 1 could enjoy the merits of FImP. Last but not least, we actually will need two scalar singlets for the sake of successful EWSB with acceptable Higgs phenomenologies, and this will make a difference in freeze-in.
The paper is organized as follows. In Sect. 2 we introduce the νSISM, studying EWSB, the mass spectrum, and freezein dynamics. We present the numerical results in Sect. 3, and the next section includes a conclusion and a discussion. Finally, some supplementary materials are in the appendix.

The scale invariant version of νSM (νSISM)
As stated in the Introduction, the νSISM predicts the FImP DM N 1 , which takes advantage in addressing the basic properties of DM such as mass and relic density origins and stability. In this section we will first introduce the model and then study its two main phenomenologies, SI spontaneously breaking and the N 1 freeze-in production. The scalar singlets will play crucial roles in both aspects.
We begin with the model setup. Asides from the SM field content, the νSISM introduces three RHNs, N i , to produce realistic neutrino masses and mixings. However, in order to generate Majorana masses for N i , scalar singlets, S a (a = 1, 2 . . . n), with non-vanishing VEVs are indispensable. Interestingly, they are also badly needed to implement radiative SI spontaneously breaking. It is found that the minimal case with n = 1 fails to accommodate the current Higgs data [43,44], so we consider n = 2 real scalars or a complex singlet S = (J + iσ )/ √ 2. These minimal degrees of freedom are subject to the SU (3) C × SU (2) L × U (1) Y gauge symmetries and the Poincaré and classical scale invariance space-time symmetries. Then the most general interacting Lagrangian reads To reduce parameters, we also impose CP-invariance on the scalar potential, which forces λ 4,5,6 to be real. But this symmetry is not physically necessary (we allow a complex Yukawa coupling to break it explicitly). For later use the expansion of Eq. (1) in real degrees of freedom is given in Appendix A.
A comment about the novelty of the νSISM is in order. As a matter of fact, it is not entirely new and similar versions have been investigated in Refs. [43][44][45]. However, our physical motivations and the arguments of the model presented here are quite different from theirs. In addition to that, our focused parameter space will be totally different.

Flat direction and tree level spectrum
We proceed to investigate the vacuum of the scalar potential given in Eq. (1). In the absence of a symmetry, it is supposed that both components of S develop VEVs. So, along with the Higgs doublet VEV, we should work in the threedimensional field space. The treatment of such a situation becomes more complicated than the single field case, and the Gildener-Weinberg [46] instead of the Coleman-Weinberg (CW) [47] approach is appropriate. Following this approach, we should work out the flat directions of the scalar potential, namely flat = ϕ n with ϕ and n, respectively, the modular and directional vector of flat in the multi-dimensional field space. The existence of such a kind of directions is a consequence of SI at tree level, which leads to the modularindependent minimum lines of the potential, i.e., the lines satisfying the equation dV tree /dϕ = 0. The flat directions will be lifted by radiative corrections that introduce sources of SI violation and then fix the value of ϕ. We recall that these operations are employed at the renormalization scale Q, at which all couplings are inputted. Later, it will be chosen at the SI spontaneously breaking scale.
We have simplified the potential by imposing CPinvariance, and this simplification allows us to obtain analytical expressions of the flat direction (and mass spectrum as well). Such a strategy is also adopted in Ref. [45], but it further imposes a Z 4 discrete symmetry to eliminate the λ 4 -and λ 5 -terms. The flat direction is the solution of the following three tadpole equations: where we have introduced two VEV ratios x ≡ σ/J and y ≡ J/ h. As a convention, we always choose J as the largest scale and thus y > 1 (we will turn back to it later). The above equations are derived from the potential expressed in terms of real components; see Eq. (A2) and the definitions of coupling constants Eq. (A3). Immediately, from these equations one can see that in the Z 4 -symmetric case, where λ 4 = λ 5 = 0 and thus λ J = λ σ , λ h J = λ hσ , only solutions x = ±1 are possible. In other words, Z 4 forces the flat direction to lie along the special direction σ = J . But more generic flat directions can be admitted in the absence of such a symmetry. The resulting Higgs physics is different. In our later numerical demonstration, the VEV ratios x and y will be taken as fixed inputs, since they have direct physical meaning. Note that y can be expressed in terms of x and other quartic coupling constants, In all cases of interest we shall find that y is at least moderately larger than 1, and thus there should be a mild hierarchy, i.e., λ h J , λ hσ λ h . Otherwise we have to arrange a cancellation between the two terms in the numerator. It is convenient to express the flat direction in terms of x, y, and h as follows: By definition, we have the relations v h,J,σ = n h,J,σ ϕ. Later we will find that R(x, y) (or n h ) has a clear interpretation.
The tree level spectrum along the above flat direction consists of three Higgs states, the Goldstone boson (GSB) of SI spontaneously breaking P and two massive states H 1,2 . We should calculate radiative corrections of the tree level potential due to these states. The mixings between these scalars are largely determined at tree level. Concretely, the three states are obtained from the mass squared matrix M 2 , which in the basis (h, J, σ ) takes the form of We have expressed λ J and λ σ in favor of others via the tadpole conditions Eqs. (3) and (4), respectively. One can explicitly check that DetM 2 = 0, hence implying a GSB.
The two massive Higgs bosons can be analyzed in a simple way. The Higgs sector must present a quite SM-like Higgs boson h SM , which must be dominated by the Higgs doublet, merely containing small fractions of singlets. It means that the Higgs sector can be split into a doublet sector and a singlet sector. In this simplified case, the singlet sector breaks SI through large singlet VEVs, which are then directly mediated to the doublet sector through the mixing term S 2 |H | 2 , generating the negative Higgs mass term for EWSB. 3 Now, the mass squared of h SM can be approximated by So as usual λ h ≈ 0.75 is almost fixed by the 125 GeV Higgs boson, up to a small correction from the mixing effect. With this result and taking into account the presence of a GSB, which is dominated by singlets, it is not difficult to deduce the mass squared of the heavier Higgs boson H 2 : In the ensuing discussion we will see that, to guarantee the presence of h SM , at least one singlet should develop a larger VEV than that of the doublet, v. As mentioned before in our notation we take v J = J v or y = v J /v 1. Then the mass of H 2 is well approximated by Or one can express it in terms of v ϕ , m H 2 = x λ σ 3 + λ hσ x 2 y 2 v ϕ , by solving λ J σ from Eq. (4). In most cases it can be further approximated to be of the very simple form x √ λ σ /3v ϕ .

CW potential and pseudo-GSB mass
With the massive spectrum at hand, we now discuss SI spontaneous breaking in the one loop CW effective potential V CW . In general, it can be written as [45] The numerical factors A and B are functions of the dimensionless couplings only. In the MS scheme, they are explicitly given by with s/V / f denoting real scalar/vector boson/Dirac fermion, respectively. In writing these expressions, we have collectively assumed the relation m = gv ϕ , which is true in the models with SI. It is worthy of note that in the multidimensional field space, if the field that gains mass mainly via coupling to the scalar with a smaller VEV, e.g., h in this context, the resulting effective coupling g will be suppressed by the fraction n h 1. In this way, we can avoid the failure of the scale invariant SM, which is caused by the large negative contribution to B from the heavy top quark. We will come back to this point later.
Within the perturbative region, namely both A and B much smaller than 1, the VEV of ϕ is pinned down by solving the equation dV CW /dϕ = 0. From this we get the relation If we choose Q = v ϕ , which eliminates the potential large logarithmic terms, we get a relation among the couplings at scale Q, i.e., A + 2B = 0. It is nothing but the phenomenon of dimensional transmutation [47]. Perturbativity of the effective potential is satisfied for g s 1. If both A and B receive contributions mainly from a single field dependent mass term, such as m H 2 in this paper, we can approximately determine the effective coupling constant to be or m s = √ e Q. Equation (15) tells that, if we choose Q v ϕ , then g s 1, and thus it violates perturbativity at Q. In the opposite case, one may want to choose Q v ϕ to get a weak coupling g s 1. But it will result in an intolerably light PGSB with mass suppressed by Q/v ϕ (or g s ). Therefore, (11)); thus we have g H 2 ≈ x(λ σ /3) 1/2 . Then to make g s ≡ g H 2 ∼ 1 we need λ σ ∼ 3/x 2 , but it blows up even for x 0.3. So the vacuum figuration with x 1 is favored, which means that H 2 tends to a strong mixture of σ and J . Now we turn our attention to the pseudo GSB (PGSB), P, which is massless at tree level but gets a mass from the effective potential. Before heading toward its mass, we would like to first figure out its doublet fraction F Ph , which is useful in constraining/discovering this PGSB at colliders. After an explicit calculation, one finds that this fraction is nothing but just the doublet projection in the tree level flat direction, see Eq. (6), up to an overall sign. F Ph is determined solely by the VEV ratios x and y. A small F Ph is necessary not only to hide a fairly light PGSB at LEP but also to keep h SM overwhelmingly dominated by the doublet. To sufficiently suppress it, at least one of the singlet VEVs should be significantly larger than v, i.e., y 1 and/or x y 1. Following the previous convention, we take y 1 and then we have Note that in this convention x ≤ 1, thus x y 1 also means y 1. The classical SI is violated by quantum effect which thus generates a mass for the tree level massless GSB. In general, it is given by To ensure that the extremum from dV CW /dϕ = 0 is indeed a minimum, m 2 P or B must be positive. In the above expression, top quarks have the potential to drive B < 0 but it is stopped by H 2 . The stability condition is .
It is actually required that λ σ > 6.3/y 2 x 2 (1 + x 2 ). No surprise, when x < 1/y 1 (or x 2 y 2 1), namely one of the singlet having VEV below the weak scale, that one needs a large λ σ ∼ O (10) to compensate for the relatively larger suppression by n σ in g H 2 (compared to that by n h in g t ). This is not an appealing situation if we want to keep the model perturbative up to a very high scale. Moreover, the resulting spectrum, in particular P, is fairly light and thus may have already been excluded by the present experiments like LEP and LHC. 4 Therefore, we will focus only on the case x 2 y 2 1, consistent with the analysis below Eq. (15), which concludes that x should be near 1. In the H 2 -dominance limit, after using Eq. (15), the PGSB mass is expressed in the following form: As an estimation, we write m P = 75.7×(y/5.0)(λ σ /2.0) GeV with x = 0.8 fixed. Increasing y helps to not only suppress the doublet fraction of P but also enhance the mass of P, thus making it safe under the LEP exclusion. The price is pushing H 2 into the TeV region, thus making it hard to detect at near future colliders. But P is still promising. We leave more quantitative analysis in the coming section and in the ensuing subsection we enter into a discussion of dark matter, the FImP N 1 .

Freeze-in sterile neutrino
Before heading toward the freeze-in production of sterile neutrino dark matter N 1 , we briefly discuss its conventional production mechanisms. The tiny active neutrino mass leads to a naive upper bound on the Yukawa coupling constant given in Eq. (1) (the bound may be spoiled somehow in the presence of flavor structure), That feeble coupling means that N 1 never enters the thermal plasma, given no other interactions. It can be nonthermally produced via non-resonant sterile-active neutrino oscillation, known as the Dodelson-Widrow (DW) mechanism [48]. But it has been ruled out (we will give reasons later). The resonant production [49] and thermal production mechanisms [50] are still allowed. However, both suffer some theoretical defects since they require big modifications. The former requires an anomalously large lepton asymmetry, and the latter requires a large entropy release. Model extensions are then unavoidable. 5 In the νSISM, new interactions of N 1 with the singlets appear naturally. These new interactions originally are supposed to generate very light mass for N 1 (lightness is necessary for stability), but at the same time they surprisingly can account for the correct relic density of N 1 via the freeze-in mechanism.

A toy analysis
As a toy analysis, only one singlet is considered for the time being. For the keV scale N 1 , again for stability (and for cosmological considerations [51]), the strength of the coupling is extremely small λ n ∼ 10 −8 . Hence this new vertex cannot thermalize N 1 either. But the magnitude of λ n is just at the correct order to admit freeze-in production of N 1 . In this scenario, N 1 is produced via the slow decay process S → N 1 N 1 , which has a negligible inverse decay rate. Note that generically the annihilation processes like ff →N 1 N 1 contribute to the freeze-in process sub-dominantly [52,53], due to the reason, among others, that they are suppressed by extra small couplings. The freeze-in process lasts until the mother particle decouples from the plasma and quickly decays away. The peak of production is around the mass scale of the mother particle S. In other words, it is UV-insensitive (for UV-sensitive freeze-in, please see Ref. [54]). For freeze-in proceeding via two-body decay P →X X, the final yield of X is formulated to be [32,33,52,53,55,56] with g S,ρ * , respectively, the effective numbers of degrees of freedom for the entropy and energy densities at T m P , the mass of the mother particle. g P is the number of internal degrees of freedom of P. For multiple mother particles contributing to freeze-in X , there is a summation over P. Specified to this schematic example for freeze-in N 1 , the relic density depends on the unknown parameters as proportional to λ 3 n S . Then, for a TeV scale S , it is found that a keV scale FImP allows for correct relic density. This kind of feeble interaction admitting a correct relic density of extremely light DM is somewhat reminiscent of a weak interaction admitting the correct relic density of the weak scale DM. Here the TeV scale, a sign of new physics, is also involved, but it is in the form of a VEV rather than the DM mass scale itself.
The analysis in the realistic model becomes more complicated in two respects. First, there are two singlets J and σ coupling to N 1 , both having VEVs. Moreover, two physical Higgs bosons P and H 2 (actually three but the contribution from the SM-like Higgs boson is suppressed by small mixing) contribute to the freeze-in process. Second, it is well known that a single family of RHN fails to accommodate neutrino phenomenologies, and thus a nontrivial flavor structure should be taken into account. This, along with the multi singlets, is going to make a big difference. We find that there are two distinguishable scenarios that can successfully freeze-in N 1 , and one of them is just reduced to the single RHN case.

Consider multi-singlets and flavor structure
As a warm up, we consider the case with only one RHN N 1 but two singlets. The Lagrangian relevant to freeze-in is derived as follows: with S α = (J, σ) and M N 1 = η S α v S α = (x η J + η σ )v σ . 6 We have written S α = F S α H a H a with H a = (P, H 1 , H 2 ).
In the region of interest with x 1 and y 1, the H 1 (= h SM ) component can be neglected. Other fractions are approximately parameterized by one mixing angle θ between the singlets: F J P = F σ H 2 = cos θ and −F J H 2 = F σ P = sin θ with tan θ = x and θ ∈ [0, π/4]. It is illustrative to Barring cancellation between the two contributions to M N 1 , one gets the naive estimations η σ,J ∼ M N 1 /v σ,J . In particular, if only one singlet couples to N 1 , there will be no cancellation and then for x 1 it is expected that f H α ∼ 1, the reference value for f H α hence. Now substituting the decay 6 The SM Higgs doublet also contributes to the mass of N 1 via the dimension-five operator |H | 2 N 2 1 / , which is obtained after integrating out the singlet S with VEV v s and mass m S . Roughly, with m DM = M N 1 .
From the above equation we find that there is a mild tension between the Higgs sector and the dark sector. The tension gets more serious as N 1 becomes lighter, e.g., as light as the potential warm dark matter with mass around 1 keV. For that light DM, Eq. (24) shows that the dark sector wants the singlets' VEV to lie significantly below the TeV scale, which, however, is disfavored by the Higgs phenomenology in terms of the previous discussions. Therefore, we may have to endure a substantial cancellation between the two singlets coupling to N 1 , so as to make at least one f H a ∼ 10. Immediately, we know that the case with only one singlet coupling to RHN fails in freezing-in a quite light N 1 , because it always gives f H a ∼ 1. Now let us detail how incorporating the flavor structure for RHNs makes a big difference. Actually, it opens up a novel scenario for freeze-in. To see it, we consider a simplified case, i.e., there is a large mass hierarchy between the DM candidate and the extra RHNs; for instance, the other two RHNs lying at the GeV scale inspired by baryogenesis [66]. The genetic Yukawa couplings are η S α ,i j S α N i N j /2 + c.c. with η S α ,i j = η S α , ji and i/j = 1, 2, 3. After writing S α → v S α + S α , we as usual can work in the mass eigenstates of RHNs through an unitary rotation N i → U i j N j and eventually arrive at the Lagrangian with Owing to multi singlets with VEVs, the mass matrix and the Yukawa coupling matrix of RHNs cannot be diagonalized simultaneously. Consequently the interaction of N 1 is more involved than that of the previous case, which is only the limit of a negligible mixing effect in the case considered here. To understand this limit better and for later convenience, we come back to the original mass matrix M N (for illustration only two RHNs are considered): with −3 a 0. To make the lighter eigenvalue naturally ∼ O(10 −6 ) GeV, b should be smaller than a − 3. The single RHN case discussed before corresponds to b a − 3, a condition to decouple the heavy from the light.
On the contrary, when b a − 3 the light-heavy RHN mixing effect is not negligible and can play an important role. Still we consider the hierarchical scenario, where the mixing angle is estimated to be sin θ N (M N ) 12 1. Two off-diagonal Yukawa coupling constants η H a ,12 , by naive estimation, are approximated to be After some exercise, it is not difficult to find that at leading order η P,12 = 0. But the one involving H 2 is not zero, and explicitly it is given by 22 ) x η σ,12 + η J, 12 x η σ,22 + η J,22 .
As one can see, η H 2 ,12 ∼ η S α ,12 ∼ 10 b /v S α . It is b − 6 orders of magnitude larger than η H a ,11 , which is expected to be ∼ 10 −9 for a keV RHN. Recall that b < a −3, thus as long as a is not less than −2, η H 2 ,12 can easily be around 10 −8 . In the light of the estimation in Eq. (24), this light-heavy mixing effect is able to freeze-in a quite light (∼keV) N 1 via decay H 2 → N 1 N 2 . We stress again that the presence of multi singlet scalars is key to preserving a significant mixing effect η H a ,12 = 0. In summary, in our paper the correct relic density of N 1 , no matter light or heavy, can be readily achieved by means of freeze-in.

The numerical analysis of the νSISM
In this section we investigate the numerical results of the Higgs and dark sectors phenomenologies, respectively. For the former, we are not interested in making a parameter space scan. Instead, we try to make the features of the Higgs boson spectrum visual. For the latter, we even are not concerned with the original parameters, since from the previous analysis it is known that the correct relic density of N 1 can always be achieved over a fairly wide region of M N 1 . So the dark matter phenomenologies actually involve only two parameters, M N 1 and sin θ 1 , the mixing angle between N 1 and active neutrinos. Noticeably, the recently reported X-ray line at 3.55 keV can be readily explained by the decaying N 1 (7.1 keV DM with freeze-in production to explain this line was also considered in Refs. [57,58]). It is problematic for N 1 produced by other mechanisms.

On the Higgs bosons
The analysis of the Higgs sector includes a theoretical constraint, i.e., the potential should be bounded from below (BFB) in order to make the minimum stable. We check this at tree level. It is more convenient to do so in the three-dimensional field space spanned by (h, J, σ ), where the potential can be written as V = X X T with X = (h 2 , J 2 , σ 2 ) the bilinear vector and the matrix of coupling constants being BFB requires that all of its sub-matrices n×n with n = 1, 2, 3 should satisfy the conditions Tr n×n > 0 and Det n×n > 0. Equivalently, the following conditions for the quartic couplings should be fulfilled: Basically, they require sufficiently small off-diagonal elements in . Note that they practically lead to mass mixings in the Higgs mass matrix; thus the above conditions are consistent with the requirement that there should be a quite SM-like Higgs boson.
Aside from the renormalization scale Q, altogether the Higgs sector contains six real parameters, λ h,σ,J and λ hσ , λ h J , λ J σ . Among them, λ h is almost fixed by the SM-like Higgs boson mass, and three can be expressed in terms of two VEV ratios x, y, and λ hσ , λ J σ . Since x and y have clear physical implications, they are taken to be inputs and fixed in our numerical samples. Finally, only three free parameters, λ σ , λ hσ , and λ J σ , are left. We study the Higgs spectrum varying with them and display the results in Fig. 1, on the λ hσ − λ J σ plane, with λ σ chosen schematically. For each chosen λ σ (dashed line), we show the corresponding mass spectrum through two lines, the thick line and dotted line for m H 2 and m P , respectively. For comparison, in Fig. 1 we show two choices of y: y = 5 (left panel) and y = 3 (right panel). Some observations are in order from the figure: • In the singlet-doublet decoupling limit, m h SM in the expectation is merely sensitive to the diagonal quartic BEFB is able to rule out the smaller values of λ σ for a larger y(= 5), but for the smaller y(= 3) this power tends to be lost. For illustration, in the y = 5 case we show three typical solutions with λ σ = 0.5, 2.0, and 4.0, respectively. One can see that given m h SM = 123-126 GeV, the first one has been excluded by BEFB, the second one is near exclusion, while the last one is safe. • The region giving m P < 114.4 GeV is subject to the LEP constraint and (m P , F Ph ) are restricted. In turn, λ σ and y are bounded. We do not intend to scan the whole parameter space since the studies of the SM Higgs doublet mixing with a singlet scalar have been done in many works, say, the most relevant one [43,44]. Generically, as long as λ σ is around 1 and y is relatively large, P can easily avoid the LEP bound; if and only if y and λ σ are of normal size, there is hope in the hunt for at least one of the two extra Higgs bosons. If in the future we do really hunt for two new Higgs bosons with hierarchical masses; the νSISM will be a good candidate to account for them.

On the FImP with a benchmark at 3.5 keV X-ray line
Before heading toward the freeze-in production, we explain why N 1 with the DW mechanism has been ruled out already. N 1 can decay into a neutrino plus photon, with decay width [2] N 1 →νγ where θ α1 is the mixing angle between N 1 and the active neutrino flavor ν α . The nonobservation of an X-ray line stringently constrains on M N 1 and α sin 2 θ α1 ≡ sin 2 θ 1 [59]. The width is proportional to M 5 N 1 , so θ 1 is restricted to be very small for a heavier N 1 . Consequently, the final yield of N 1 via the DW production is insufficient [2]: Concretely, the region M N 1 3-4 keV has been excluded. On the other hand, the Lyman-α (Lyα) forest gives a compensatory constraint on the lighter RHN which has a longer length of free streaming λ fs . The latest analysis yields M N 1 8 keV [60]. Therefore, the entire region of M N 1 has been excluded; see Fig. 2.
The DM production via freeze-in changes the situation dramatically. First, the X-ray bound can be avoided because In this plot, the shadowed region is excluded by the X-ray observations, and the vertical line at M N1 = 3 keV indicates the lower bound on M N1 from the free streaming limit. On the line with lepton asymmetry L = 0, a correct relic density of N 1 is achieved via the DW mechanism; while on the line with L = 1.24 × 10 −4 , it is via resonant production (we use data from Ref. [66]). The red star labels the point that fits the 3.55 keV X-ray anomaly. It lies in the bulk parameter space of the freeze-in scenario the freeze-in production mechanism has nothing to do with the sterile-active neutrino mixing. In principle, it can be arbitrarily small given three families of RHN, because even in the limit of decoupling N 1 from the active neutrinos the other two heavy RHNs can still account for the neutrino masses and mixing. 7 Second, the Ly-α bound relaxes significantly because of two reasons. One is that f ( p, t), the initial spectrum of N 1 , becomes slightly colder in the freeze-in scenario, where the average momentum of f ( p, t) is p T ≈ 2.45 T, while in the DW scenario p T ≈ 2.8 T [39]. The other one is due to a significant entropy dilution. Here N 1 is produced at the very early Universe t in , corresponding to the temperature T in m H a ∼ 0.1-1 TeV. From t in to the current time t 0 there is a large entropy release S ≡ g * (t in )/g * (t 0 ), which substantially cools N 1 down. Eventually, the average momentum becomes The Ly-α bound accordingly weakens and merely gives M N 1 >1 keV, a fairly loose bound. But to avoid a hot N 1 yields a stronger lower bound. In terms of Eq. (B4) the free streaming of N 1 is estimated to be 8 7 Such a limit amounts to the singlet fermonic FImP via a singlet scalar portal [56]. Obviously, in that limit one does not need to worry about the X -ray bound. 8 Different from the freeze-in scenario through a frozen-in scalar [37,42], here the length of free streaming is independent on mass of the decaying scalar boson. λ fs = S −1/3 √ a eq t eq √ t nr 5 + ln t eq t nr 0.038 × 100 g * Mpc, (34) where √ t nr is the time at which RHN becomes nonrelativistic. To ensure that N 1 does not turn out to be hot, one requires λ fs 0.1 Mpc and in turn M N 1 3 keV.
Interestingly, a X-ray line at energy 3.55 keV was recently reported with 3σ significance evidence [34,35], through the observation of galaxy clusters and the Andromeda galaxy. Despite the controversy [61][62][63][64], it is tempting to interpret it as a smoking gun of decaying sterile neutrino with mass about 7.1 keV and mixing angle sin 2 2θ 1 ≈ 7×10 −11 [34,35]. Production mechanisms of the correct relic density for such N 1 are relevant here. The non-resonant production fails already. The resonant production mechanism may work [65], but the latest work Ref. [38] showed that it also fails after taking into account the Ly-α bound. The freeze-in mechanism that gives a colder RHN, either via a frozen-in [37,38] or a thermal scalar boson decay, successfully accommodates that N 1 with correct relic density in the bulk parameter space; see Fig. 2.

Conclusion and discussion
We proposed the FImP framework for dark matter, a combination of FIMP with SI. It is consistent with the null results from all kinds of DM detections. Besides, it shows advantages in addressing basic questions about DM, stability, mass origin, and relic density generation, in an inherent way. In the golden example, the νSISM, the FImP candidate, the lightest RHN N 1 , furthermore is predicted instead of introduced artificially. We would like to stress that another attractive feature of the νSISM is its economy and self-consistence. Owing to scale invariance, it is necessary to incorporate extra singlets that develop VEVs around the TeV scale to generate Majorana masses for the RHNs. At the same time, they are capable of producing N 1 via freeze-in, addressing its relic density problem. It is not only so that these singlets are badly needed for SI spontaneously breaking itself.
Two open questions deserve further exploration. First, in this paper we actually work in the three RHNs scenario, so the X-ray bound is simply gone in the decoupling limit of N 1 . But it is tempting to work in a more predictive framework where N 1 plays a more active role in neutrino physics, e.g., only two RHNs are introduced and then N 1 's cannot decouple from active neutrinos. In that case, the resulting X-ray line will be closely correlated with neutrino phenomenologies. Second, it is of special interest to explore sterile neutrino DM in the scale invariant B − L models [4,[67][68][69][70] where the RHNs have a more natural physical origin, i.e., they are required by anomaly cancellation. But in such kinds of models RHNs carry B − L charge and thus they are thermalized, except in the limit of decoupling new gauge dynamics, e.g., the new gauge coupling is vanishingly small or the massive gauge boson is very heavy, in such a way that it decouples before reheating.
redshifts the momentum of dark matter by a factor S −1/3 as indicated above. It is not surprising that the above expression is the same as the one derived in Ref. [42], since the average momentums have the same scaling, behaving ∝ T . The difference is manifest in the expression in t nr .