Upgrading Sterile Neutrino Dark Matter to FI$m$P Using Scale Invariance

In this article we propose a class of extremely light feebly interacting massive particle, FI$m$Ps. 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 ($\nu$SISM), the lightest $N_1$ realizes the FI$m$P scenario. In this example scalar singlets, which are intrinsic to the $\nu$SISM, generate mass and relic density for this FI$m$P simultaneously. Moreover, they are badly needed for electroweak symmetry spontaneously breaking. Interestingly, a 7.1 keV $N_1$ with correct relic density, that can explain the recent 3.55 keV $X-$ray line, lies in the bulk parameter space of our model.


I. INTRODUCTION AND MOTIVATION
Developments of particle physics are both theoretically and experimentally oriented. Therefore, when we are blindly exploring new physics beyond the standard model (SM), we should try our best to keep eyes on guidelines from both sides. Thus far, the hierarchy problem and existence of dark matter (DM), nonzero neutrino masses respectively provide the strongest guidelines from each side. A framework that could address them unifiedly, out of question, must be of special interest. But usually they are not inherently related. Supersymmetry (SUSY) is a famous example. It can elegantly solve the hierarchy problem and at the same time present a DM candidate, but it is not a theory of neutrinos. In this article we will follow another line, i.e., the (classical) scale invariance (SI) that is potential to be a protecting symmetry on the weak scale [1] and inspires a number of studies [3][4][5][6][7][8][9][10]. Along this line, we will find that those three aspects can be closely related.
We start from some thoughts about DM, from which we will spell out how they are related. Despite of a lot of experimental efforts, people still have not acquired much confirmative information about the particle natural of DM. Since the DM (in)direct detections are mainly based on the weakly interacting massive particle (WIMP) DM hypothesis, the null results [2] begin to challenge this paradigm. Now it is at the right time to reexamine the theoretical basis of WIMP DM. Or more ambitious, we may should figure out what kind of DM is theoretically preferred and accordingly suggest new focus of DM detections. We would like to follow the ensuing basic questions about DM.
• Firstly, why is it there? Obviously, a DM candidate that is predicted rather than introduced would be more attractive. A good case in point is the lightest sparticle (LSP) in SUSY and axion from the Pecci-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 [11]. RHNs are originally introduced to explain the nonzero neutrino masses, but they are neutral and thus potential to provide a DM candidate.
• Secondly, why is it stable (or at least sufficiently stable)? 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 SI [10] 1 . In this sense even the LSP, the most promising WIMP DM candidate, is criticized for its demand for the so-called R−parity. By contrast, the sterile neutrino N 1 , given a mass ∼ O(keV), can have lifetime at cosmic time scale even in the absence of an exact protective symmetry. Note that in this way of understanding DM stability, DM is expected to be decaying, and DM lightness (far below the weak scale) plays the key role in suppressing the decay rate.
• Last but not the least, why the DM relic density Ω DM h 2 ∼ 0.1? In the WIMP sceanrio, the answer is illustrative: DM annihilates involving weak scale mass and weak scale strength, which is known as the WIMP miracle and is just the reason that makes WIMP gain advantage over others. This article discoveries an even stronger miracle for the feebly interacting massive particle (FIMP) [12] with mass around the keV scale 2 . Sterile neutrino N 1 is a very good case in point and we will focus on it throughout the paper, but the discussion can be easily generalized to other particles (the scalar FIMP case can be found in Appendix C).
The FImP ("m" instead of "M" is used to indicate that the massive particle is as light as the keV scale) miracle heavily depends on SI, hence tied with the solution to the hierarchy problem. Simply, SI forbids the direct Majorana mass terms for N 1 and therefore scalar singlets S a with nonzero vacuum expected values (VEVs) are indispensable to give mass to N 1 , via terms λ a S a N 2 1 /2. On the other hand, S a are also badly needed to realize the hidden Coleman-Weinberg (CW) mechanism [13] that could radiatively break SI by nonzero S a and then induces electroweak symmetry spontaneously breaking (EWSB). It is found that S a ∼ O(TeV) are natural values to accommodate a SM-like Higgs boson around 126 GeV. Thus a keV N 1 means the coupling constants λ a ∼ O(10 −9 ). Such a small strength in turn is just able to produce relic density of N 1 at the correct order, through the freeze-in mechanism [12]. As one can see, the FImP miracle by virtue of SI is even stronger than the WIMP miracle, because in it (the keV) DM mass origin surprisingly implies the correct order of relic density.
After these thoughts, we draw the conclusion: N 1 in the scale invariant SM with RHNs (νSISM) is provided with all the theoretical merits, predicted, naturally stable, and moreover endowed with a relic density miracle thanks to SI. In the νSISM, DM and neutrino are unified. Moreover, the particle features of DM like mass and interactions, to some degree, are specified by SI. Therefore, it is an example to show the inherent relations between the three aspects referred at the beginning of the paper.
In the late 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 [14][15][16], as an alternative way to achieve correct relic density for N 1 (note that the formulation in this paper is a little different to theirs). But in our νSISM, as a consequence of SI, such a singlet is a built-in instead of added element. It further gives rise to the FImP miracle, which can be realized only by SI. Interestingly, such a FImP miracle may already leave hint in X−ray line observations. Recently a X−ray line at energy 3.55 keV is reported [17], that can be explained by a sterile neutrino DM with 7.1 keV and mixing angle with active neutrino around 0.8 × 10 −5 [17]. Such a sterile neutrino is problematic in acquiring correct relic density via the conventional ways, but is ready via freeze-in [20][21][22].
The paper is organized as follows. In Section II we introduce the νSISM, studying EWSB, the mass spectrum and freeze-in dynamics. We present the numerical results in section III, and its next section includes conclusion and discussion. Finally, some supplementary materials are casted in the appendix.

II. THE SCALE INVARIANT VERSION OF νSM
As stated in the Introduction, the νSISM demonstrates the paradigm that can inherently address the most guiding aspects of new physics, the gauge hierarchy problem, DM and neutrino masses. Moreover, DM, the lightest RHN N 1 , takes theoretical advantages. In this section we will first introduce the model and then study its two main phenomenological aspects, SI spontaneously breaking and the N 1 freeze-in production. The singlets play crucial roles in both aspects and create the FImP miracle for N 1 . Obviously, this model is a surprisingly economic model, in which each member plays a very important role.
We begin with the model setup. Asides from the SM field content, the νSISM introduces three RHNs N i to produce neutrino masses and mixings. Besides, in order to generate Majorana masses for N i , singlet(s) S a are indispensable. At the same time, the singlets are capable of triggering SI and EW symmetries breaking. It is found that the minimal case which introduces only one real singlet fails to accommodate the current Higgs data [9], so two real singlets or a complex singlet S = (J + iσ)/ √ 2 will be introduced. For such minimal degrees of freedom subject to SU (3) C × SU (2) L × U (1) Y gauge symmetries and the Poincare and classical scale invariance space-time symmetries, we can write down the most general interacting Lagrangian, We have imposed CP-invariance on the Higgs sector, but practically it is not necessary and is merely for the sake of reducing parameters. The most general potential expressed in both complex and real degrees of freedom is casted in Appendix. A. A comment about the novelty of the νSISM is in order. As a matter of fact, it is not brand new and its similar version has been investigated in Ref. [8,9]. However, our physical motivations and arguments of the model are quite different to theirs. Additionally, our focused parameter space will be totally different.
A. Scale invariance spontaneously breaking

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, both components of S are supposed to develop VEVs, which is important to implement the FImP miracle. So, along with the Higgs doublet VEV, we should work in the three-dimensional field space. Treatment with such a situation becomes more complicated than the single field case, and the Gilderner-Weinberg approach [23] is suitable for it. Following this approach, we should work out the flat directions of the scalar potential, namely Φ f lat = ϕ n with ϕ and n respectively the modular and directional vector of Φ f lat in the multi-dimension (three in our case) field space. The existence of such kind of directions is a consequence of SI at tree level, which implies the modular-independent minimal lines of the potential, i.e., they satisfy the equation dV tree /dϕ = 0. The flat directions will be lifted by radiative corrections, that introduce sources of SI violation and then fixes the value of ϕ. We remind that these operations are employed at a scale Q, at which all couplings are inputed.
For the model under consideration, we have simplified the potential by imposing CP-invariance. This simplification allows us to obtain simple analytical expressions of the flat directions (and mass spectrum as well). Such a strategy is also adopted in Ref. [8], but it further imposes Z 4 discrete symmetry to eliminate the λ 4 − and λ 5 −terms. Interestingly, we will find that this Z 4 will force the flat direction to lie along the special direction σ = J, however, more generic flat directions can be admitted without Z 4 . And they show difference in Higgs physics. Without loss of generality, the flat directions are solutions of the following three tadpole equations: where we have introduced two VEV ratios x ≡ σ/J and as well y. The equations are derived from potential in terms of real components, see Eq. (A2) and the definitions of coupling constants Eq. (A3). Immediately, from the above equations one can see that for the Z 4 −symmetric case, where λ 4 = λ 5 = 0, only solutions x = ±1 are allowed. In our treatment, y will be considered as an intermediate variable and can be expressed in terms of x and other quartic coupling constants, Thus, to have y ≫ 1, the case interested in this article, the quadratic coupling constants between the doublet and singlets (λ hJ and λ hσ ) should be much smaller than λ h . Otherwise we have to arrange a subtle cancelation between them. It is convenient to express the flat directions in terms of x, etc., as the following , n J = y R(x, y) , n σ = xy R(x, y) .
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 . They furnish the basis of calculating radiative corrections on the tree level potential. The mixings between them are also largely determined at tree level. Concretely, they are obtained from the mass squared matrix M 2 , which in the basis (h, J, σ) takes a form of We have expressed λ J and λ σ in favor of others via the tadpole conditions Eq. (3) and Eq. (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 is dominated by the doublet and its mixings with other two singlets are sufficiently small. 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 via large singlet VEVs, which are then directly mediated to the doublet sector through the mixing term S 2 |H| 2 3 . Now, the mass squared of h SM can be approximated by So λ h is almost fixed by the 125 GeV Higgs boson to λ 1 ≈ 0.75, 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 difficulty to deduce 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. Without loss of generality, we will use the convention v J = J ≫ v, or equivalently y ≫ 1. Then, mass of H 2 is well approximated to be To get the second expression we have solved λ Jσ from Eq. (4). This expression will be useful and it can be further approximated to be x λ σ /3v ϕ in most cases.

CW potential and pseudo-GSB mass
With the massive spectrum at hand, we discuss SI spontaneously breaking in the one loop CW effective potential. In general, the effective potential V CW can be written as [8] 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 for real scalar/vector boson/Dirac fermion respectively. In writing these expressions, we have collectively assumed the relation m = gv ϕ , which is true in the SI models. It is worthy of note that for the field that gains mass via coupling to the scalar with a smaller VEV, e.g., h, the resulting effective coupling g is suppressed by the friction n h ≪ 1. As a result, the failure of the scale invariant SM owing to the heavy top quark is saved by the presence of a much larger singlet VEV. Within the perturbative region, namely both A and B are much smaller than 1, the VEV of ϕ is pined down by solving dV CW /dϕ = 0. From it we get the relation Alternatively, we can choose Q = v ϕ , which eliminates the potential large logarithmic terms, and then get a relation among the couplings at scale Q, i.e., A + 2B = 0. It is nothing but dimensional transmutation [13]. Perturbativity of the effective potential requires A, B < 1, which is satisfied for g s 1. Actually, in the case of both A and B receiving contributions mainly from a single field dependent mass term, such as m H2 in this paper, we can approximately determine the effective coupling constant or m s = √ e Q. Eq. (15) tells that, if we choose Q ≫ v ϕ , then g s ≫ 1, which violates purturbativity at Q. In the opposite, one may want to choose Q ≪ v ϕ to get a weak coupling g s ≪ 1. But it will result in an intolerablely light PGSB with mass suppressed by Q/v ϕ (or g s ). Therefore, g s ∼ O(1) is preferred. This fact is useful to find out the favored region of x. Recall that v σ = n σ v ϕ with n σ ≈ x, thus we have g H2 ≈ x(λ σ /3) 1/2 . Then to make g s ≡ g H2 ∼ 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 be 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 towards its mass, we would like to first figure out its doublet fraction, that is labeled as F Ph and useful in constraining/discovering this PGSB at colliders. After an explicit calculation, one finds that this faction is nothing but just the doublet projection in the tree level flat direction, see Eq. (6): up to a overall sign. F Ph is solely determined by two 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. Obviously, to sufficiently suppress it, at least one of the singlet VEVs should be properly larger than v, i.e., y ≫ 1 and/or xy ≫ 1. Following the previous convention, we take y ≫ 1 and then we have F Ph ≈ 1/y √ 1 + x 2 . The classical SI is violated by the quantum effect which thus generates a mass for the tree level GSB. In general, it is given by To ensure that the extreme from dV CW /dϕ = 0 is indeed a minimum, m 2 P or B must be positive. In the above expression, top quarks are potential to drive B < 0 but it is stopped by H 2 . The stability condition then is It actually is requiring λ σ > 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, one needs a large λ σ ∼ O(10) to compensate the relatively larger suppression by n σ in g H2 (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, only the case x 2 y 2 ≫ 1 will be of interest in this article. In the H 2 −dominance limit, after using Eq. (15) the PGSB mass is related to m H2 in a simple form: As an estimation, we write m P = 75.7 × (y/5.0)(λ σ /2.0) GeV with x = 0.8 fixed. Increasing y not only helps to suppress the doublet fraction of P but also enhances the mass of P, thus making it safe under the LEP exclusion. The price is pushing H 2 into the TeV region, hard to be detected at the near future colliders. But P is still promising. We leave more quantitive analysis in the coming section and in the ensuing subsection we enter into the discussions about dark matter.

B. FImP Miracle for the Sterile Neutrino
Before heading towards the freeze-in production of N 1 , we brief the 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 it never enters the thermal plasma, given no other couplings. N 1 can be non-thermally produced via non-resonant sterile-active neutrino oscillation, known as the Dodelson-Widrow (DW) mechanism [25]. However, it has been ruled out, see Fig. 2. The resonant production [26] and thermal production mechanisms [27] are still allowed. However, both mechanisms have some theoretical defects since they need additional physics. The former requires an anomalously large lepton asymmetry, and the latter requires a large entropy release. Model extensions are then necessary 5 .
In the νSISM new interaction of N 1 with the singlet S (Schematically only one singlet is considered for the time being), i.e., λ n SN 2 1 with S ∼ O(TeV), naturally appears. For the keV scale N 1 , which is good for stability of N 1 (and of cosmological interest), the strength of the coupling λ n ∼ 10 −8 . Hence this new vertex like the Yukawa coupling to active neutrino will not 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 , with negligible inverse decay. Note that generically the annihilation processes like ff →N 1 N 1 contribute to the freeze-in sub-dominantly [31] , 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. [29]). Generically, for freeze-in proceeding via two-body decay P →XX, the final yield is formulated to be [12,[30][31][32] Y X (∞) ≈ 45 g P 1.66π 4 g S * g ρ * Γ(7/2)Γ(5/2) 16 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 internal degrees of freedom of P . For multi mother particles contributing to freeze-in X, there is a summation over P . Specified to the schematic example for freeze-in N 1 , its relic density depends on the unknown parameters as proportional to λ 3 n S . Then, for a TeV scale S , it is not difficult to find that the FImP miracle arises.
Of course, analysis in the realistic model becomes more complicated than the toy model from two aspects. Firstly, there two two singlets J and σ coupling to N 1 , both having VEVs. And 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. Secondly, it is well known that a single family of RHN fails to accommodate neutrino phenomenologies, and thus a nontrivial flavor structure should be took into account. It, 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. So, as a warm up we first consider the case with only one RHN, N 1 . The relevant Lagrangian to the freeze-in processes is derived as the following, with S α = (J, σ) and M N1 = η Sα v Sα = (x η J + η σ )v σ . We have written S α = F SαHa H a with H a = (P, H 1 , H 2 ). In the interested region with x 1 and y ≫ 1, the H 1 (= h SM ) component can be neglected. While other fractions are approximately parameterized by one mixing angle θ between the singlets: F JP = F σH2 = cos θ and −F JH2 = F σP = sin θ with tan θ = x and θ ∈ [0, π/4]. It is illustrative to write η Ha = f Ha M N1 /v J with Barring cancellation between the two contributions to M N1 , one get the naive estimations η σ,J ∼ M N1 /v σ,J . In particular, if only one singlet couples to N 1 , there will be no cancelation and then for x ≃ 1 it is expected f Hα ∼ 1, the referred value for f Hα hence. Now substituting the decay width of H a → N 1 N 1 into Eq. (21) where m DM = M N1 . The FImP miracle manifests in the above equation. But there is already a mild tension between the Higgs sector and dark sector. It becomes more serious in the lighter N 1 region, e.g., in the potential warm dark matter region with M N1 ∼ 1 keV. The dark sector wants the singlets' VEV to lie significantly below the TeV scale, which however is disfavored by the Higgs phenomenology. So, we may have to endure a substantial cancelation between the two singlets coupling to N 1 , so as to make at least one f Ha ∼ 10. As a deduction, the case with only one singlet coupling to RHN fails in freezing-in a quite light N 1 , because it always gives f Ha ∼ 1. Now let us detail how incorporating the flavor structure for RHNs makes a big difference. Concretely, it opens a novel scenario for freeze-in. To see it, we will 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 [28]. Now the Yukawa couplings turn out to be η Sα,ij S α N i N j /2 + c.c. with η Sα,ij = η 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 ij N j , arriving at the Lagrangian Owing to multi singlets, the mass matrix and the Yukawa coupling matrix of RHNs can not be diagonalized simultaneously. The interactions of N 1 is more involved than the previous case, which is only a limit of negligible mixing effect. To under this limit better and for latter convenience, we turn back to the original mass matrix M N (for simplicity 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.
As one can see, η H2,12 ∼ η Sα,12 ∼ 10 b /v Sα . Therefore, it is b − 6 orders of magnitude larger than η Ha,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, η H2,12 can be easily around 10 −8 . According to the estimation in Eq. (24), it is able to freeze-in a quite light (∼keV) N 1 via decay H 2 → N 1 N 2 . Moreover, even in the limit only one singlet coupling to RHNs, we do not need cancelation as the case discussed before. We stress again that the presence of multi singlet scalars is a key to preserve a significant mixing effect η Ha,12 = 0. In summary, correct relic density of N 1 , no matter light or heavy, can be readily achieved by means of freeze-in.

III. 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 do not involve the original parameters, since from the previous analysis it is known that, for M N1 spanning over a fairly wide region, the correct relic density can always be achieved. Consequently, the dark matter phenomenologies actually involve two parameters, M N1 and the mixing angle between N 1 and active neutrinos sin θ 1 . Noticeably, the recently reported X−ray line at 3.55 keV can be readily explained by the decaying N 1 . While it is problematic for N 1 produced by other mechanisms.
The analysis of the Higgs sector includes one more 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 it 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 BFB requires that all of its sub-matrices Λ n×n with n = 1, 2, 3 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 are requiring the non-diagonal terms to be not very large. Note that these terms also lead to the mixing terms in the Higgs mass matrix, thus the above conditions are consistent with the Higgs phenomenologies. Aside from the renomalization scale Q, altogether the Higgs sector contains six real parameters, λ h,σ,J and λ hσ , λ hJ , λ 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 the VEV ratios have clear physical implications, thus they are taken to be inputs and fixed in our numerical analysis. Eventually, 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. To demonstrate the spectrum on a plane, λ σ has benn chosen schematically. For each chosen λ σ (the dashed line), it is accompanied by two other lines m H2 (the thick line) and m P (the dotted line) corresponding to the spectrum. Note that the λ σ taken in these three lines are slightly different, but it is enough for illustration. For comparison, in Fig. 1 we take y = 5 (left panel) and y = 3 (right panel). Some observations are available: • In the singlet-doublet decoupling limit, m hSM in expectation is merely sensitive to the diagonal quartic coupling constant λ h . Similarly, masses of H 2 and P are sensitive to the diagonal quartic coupling constant in the singlet sector λ σ but insensitive to the singlet-doublet mixing coupling constant λ hσ .
• For all of the BEFB conditions listed in Eq. (30), practically the last one suffices. It excludes the regime outside of the green-shadowed area. We call the intersection between the Higgs mass curvature and λ σ −line, the black circle, as a solution which gives a set of parameters (λ σ , λ hσ , λ Jσ , ...). Now it is seen that BEFB, for the larger y(= 5), is able to rule out the smaller values for λ σ , but for the smaller y(= 3) this power tends to be lost. For illustration, for y = 5 we show three typical solutions with λ σ = 0.5, 2.0 and 4.0 respectively. Given m hSM = 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 that makes m P < 114.4 GeV also subjects to the LEP constraint, which restricts (m P , F Ph ). In turn, λ σ and y will be bounded. We do not intend to scan the whole parameter space since the discussing on the SM Higgs doublet mixing with a singlet scalar has been done in many previous works, say, the most relevant one [9]. But as long as λ σ is around 1 and y is relatively large, P can easily avoid the LEP bound. Generically, only for a not very large y and λ σ , there is a good prospect to discover at least one of the two extra Higgs bosons. If in the future we do really hunt two new Higgs bosons with hierarchical masses, the νSISM will be one of the top candidates to account for them.
Now we turn our attention to the phenomenologies of dark matter N 1 . We show that the most elegant scenario, i.e., N 1 with DW mechanism, has been ruled out already. As stated in the introduction, the light FIMP usually decays because it does not require an exact protecting symmetry. Here N 1 via its mixing with the active neutrinos decays into a neutrino plus photon at two body level, with decay width [11] where θ α1 is the mixing angle between N 1 and the active neutrino flavor ν α . Observations of X−ray line constrain on M N1 and α sin 2 θ α1 ≡ sin 2 θ 1 [33]. Since the decay rate is proportional to M 5 N1 , thus for a heavier N 1 , its mixing angle is bounded to be very small. Consequently, relic density of N 1 , in the DW production [11], is insufficient: The region M N1 3 − 4 keV is 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 λ f s . According to the latest result, it gives a lower bound on the thermal N 1 , M N1 8 keV [34]. Therefore, the entire region has been excluded, see Fig. 2.
The DM production via freeze-in changes the situation dramatically. Firstly, the X−ray bound can be avoided because 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 with active neutrinos two heavy RHNs can still account for neutrino masses and mixing 6 . Secondly, the Ly-α bound relaxes significantly due to 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 DM scenario p T ≈ 2.8 T [14]. 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 Ha ∼ 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 ), that substantially cools N 1 down. Eventually, the average momentum becomes The Ly-α bound accordingly weakens and now it gives M N1 >1 keV, a fairly loose bound. But requiring N 1 not to be hot can yield a stronger lower bound. Concretely, in terms of Eq. (B4) the free streaming of N 1 is estimated to be Mpc. (34) where √ t nr is the time at which RHN becomes nonrelativistic. The derivation of this expression can be found in Appendix B. Note that different to the freeze-in scenario through a frozen-in scalar [16,21], here the length of free streaming is independent on mass of the decaying scalar boson. To ensure that N 1 does not turn out to be hot, it is required λ f s 0.1 Mpc and in turn the lower bound M N1 3 keV. FIG. 2: In this plot, the shadowed region is excluded by the X−ray observations, and the vertical line at MN 1 = 3 keV indicates the lower bound on MN 1 from the free streaming limit. On the line with lepton asymmetry L = 0, correct relic density of N1 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. [28]). The red star labelling the point that fits the 3.55 keV X−ray anomaly. It lies in the bulk parameter space of the freeze-in scenario.
Interestingly, a X−ray line at energy 3.55 keV is recently reported with 3σ significance evidence [17], through the observation of galaxy clusters and the Andromeda galaxy. Despite of the controversy [18], 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 [17]. Production mechanisms of correct relic density for such N 1 are of concern. The non-resonant production fails already. The resonant production mechanism may save it [19], but the latest paper Ref. [22] showed that taking into account the Ly-α bound, this production mechanism also fails. By contrast, the freeze-in mechanism that gives a colder RHN, either via a frozen-in [21,22] or a thermal scalar boson decay, successfully accommodates the required N 1 with correct relic density in the bulk parameter space, see Fig. 2.

IV. CONCLUSION AND DISCUSSION
The νSISM offers a good example to address the three leading hints for new physics beyond the SM in an unified framework. In this model, the right-handed neutrinos introduced to explain neutrino masses provides an accidental DM candidate, the lightest RHN N 1 with mass a few keVs. Owing to scale invariance, it is necessary to incorporate extra singlets that develops VEVs to generate Majorana masses for the RHNs. Such VEVs, viewing from the radiative breaking SI, along with the current Higgs data, are supposed to lie at the TeV scale. Then the keV scale N 1 means feeble Yukawa couplings between N 1 and the singlets. However, such feeble interactions are able to non-thermally produce N 1 via the freeze-in mechanism. It not merely provides an alternative way to alleviate the relic density problem of N 1 . As a matter of fact, here the dynamical DM mass generation and DM production share the common dynamics, in which the correct order of relic density of the FIMP is implied by the keV FIMP mass. This numerical coincidence is the FImP miracle. The recent X−ray line at 3.55 keV, which can be elegantly explained by a 7.1 keV N 1 , may be a smoking gun of the FImP miracle.
Two open questions deserve further attention. Firstly, 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 its coupling to active neutrinos can not decouple. In that case, the resulting X−ray line will be closely correlated with neutrino phenomenologies. Secondly, it is of special interest to explore sterile neutrino DM in the scale invariant B − L models where the RHNs have more natural physical origin, i.e., they are required by anomaly cancelation. But in such kind of models RHNs carry B − L charge and consequently they are thermalized, except for in the decoupling new gauge dynamics limit, e.g., the new gauge coupling is vanishingly small or the massive gauge boson is very heavy such that it decouples before the reheating temperature. We leave such interesting possibilities for further investigations.

V. ACKNOWLEDGEMENTS
We would like to thank Xiaoyong Chu and P. Ko for helpful discussions.