Axion–Sterile Neutrino Dark Matter

Abstract

Extending the standard model with three right-handed neutrinos and a simple QCD axion sector can account for neutrino oscillations, dark matter and baryon asymmetry; at the same time, it solves the strong CP problem, stabilizes the electroweak vacuum and can implement critical Higgs inflation (satisfying all current observational bounds). We perform here a general analysis of dark matter (DM) in such a model, which we call the $a\nu$MSM. Although critical Higgs inflation features a (quasi) inflection point of the inflaton potential, we show that DM cannot receive a contribution from primordial black holes in the $a\nu$MSM. This leads to a multicomponent axion–sterile neutrino DM and allows us to relate the axion parameters, such as the axion decay constant, to the neutrino parameters. We include several DM production mechanisms: the axion production via misalignment and decay of topological defects as well as the sterile neutrino production through the resonant and non-resonant mechanisms and in the recently proposed CPT-symmetric universe.

1. Introduction

Despite the remarkable success of the standard model (SM), there is no question that it needs to be extended. The observational evidence for neutrino oscillations and DM is indeed enough to draw this conclusion.

A minimal phenomenological completion of the SM up to the Planck scale was presented in [1], where the SM was extended to include three right-handed neutrinos with a generic flavor structure and the extra fields of the simplest invisible QCD axion model, the KSVZ one [2,3]. The model of [1], which we refer to as the $a\nu$MSM, not only accounts for neutrino oscillations and DM, but it can also provide the observed amount of baryon asymmetry in the universe, stabilize the electroweak (EW) vacuum, realize Higgs inflation [4,5,6,7] and solve the strong CP problem through the Peccei–Quinn (PQ) symmetry1 [11,12] at the same time.

In Ref. [13], it was found that Higgs inflation can be realized in its critical version [14,15,16] within the $a\nu$MSM: critical Higgs inflation (CHI) occurs when the SM lies extremely close to the border between the absolute stability and metastability of the EW vacuum [17]. CHI is particularly interesting for two reasons. One is that it can occur with a moderate, $\mathcal{O}\left(10\right)$, non-minimal coupling ${\xi }_H$ between the Higgs and the Ricci scalar. Consequently, the scale of breaking of perturbative unitarity, which was noticed in [18,19,20,21,22], is pushed just below the Planck scale where new physics is required to UV complete gravity. Furthermore, in Ref. [23], it was shown that CHI, unlike standard Higgs inflation [24], does not suffer from fine tuning in the initial conditions before inflation. It is also interesting that one will be able test this inflationary scenario with future space-borne interferometers [25].

To date, DM in this model has been accounted for exclusively through the axion. However, the $a\nu$MSM is rich enough to contain other potential DM candidates. DM is one of the biggest mysteries in fundamental physics, it represents the majority of matter in our universe, but its nature is still unclear. Motivated by these and other facts (see below), here we perform a general analysis of DM in the $a\nu$MSM.

We now provide an outline of this paper, which includes a summary of the results and highlights the motivations and the original parts.

In Section 2, we briefly review the $a\nu$MSM. The gauge group, $\mathrm{SU}{\left(3\right)}_{\mathrm{c}}\times \mathrm{SU}{\left(2\right)}_{\mathit{L}}\times \mathrm{U}{\left(1\right)}_{\mathit{Y}}$, is the same as that of the SM, but the field content is extended to include the three right-handed neutrinos, a complex scalar (gauge singlets) and two Weyl fermions that are charged under the color gauge factor $\mathrm{SU}{\left(3\right)}_{\mathrm{c}}$ only. The gravitational sector includes non-minimal couplings of all scalars to gravity, which allow inflation to take place. Section 2 also includes a discussion of the generic observational bounds that are needed for our purposes (other than the bounds related to DM, which are then discussed in the following sections).

Section 3 focuses on the axion contribution to DM. As explained there, we include the contribution from both the misalignment mechanism [26,27,28] and the decay of topological defects [29,30,31,32,33,34,35,36], which have been computed for the KSVZ model in [37]. The latter contribution to the DM energy density has a dependence on the quartic coupling of the extra scalar, whose value in the relevant parameter space of the $a\nu$MSM is determined here explicitly.

Section 4 and Section 5 are dedicated to the contribution to DM due to the lightest sterile neutrino. This is a good warm dark matter candidate when its mass is around the keV. Three possible mechanisms are found. The first two are the non-resonant [38] and resonant [39] production mechanisms, which occur thanks to the mixing between this sterile neutrino and the active neutrinos of the SM (see Refs. [40,41,42] for reviews). The third one takes place in a recently proposed CPT-symmetric universe [43,44], where inflation and the above-mentioned mixing are not required. For all these mechanisms, we derive the contributions to the DM energy density and the observational bounds as functions of the DM fraction $X_s$ due to the lightest sterile neutrino. Some of these functions were already known in the literature, while others are extracted here, as discussed in those sections.

Since in all the sterile neutrino production mechanisms the masses of these neutral fermions are below the ∼${10}^{14}$ GeV scale, they necessarily have a negligible impact on the running and, consequently, the parameter space of the $a\nu$MSM with absolute EW vacuum stability is enlarged [1]. This is because, generically, a Yukawa coupling (that is proportional to the mass of a fermion) contributes negatively to the $\beta$-function of the Higgs quartic coupling, as explained at the end of Section 5. Furthermore, the presence of a sizeable sterile neutrino contribution to DM, as we will discuss explicitly, allows the reduction of the mass of the extra scalar for fixed values of its couplings and so stabilization of the EW vacuum more efficiently [1,13,45,46]. All the sterile neutrino production mechanisms, therefore, favor EW vacuum stability. This is another motivation for realizing a fraction of DM through sterile neutrinos in the $a\nu$MSM.

Yet another motivation for this work is the fact that the well-motivated presence of the axion also significantly enlarges the viable region of parameter space for sterile neutrino DM in the $a\nu$MSM, compared to the case2 $X_s=1$ (where this region is quite narrow [52,53]): all observational bounds become weaker when the sterile neutrino has to account for only a fraction $X_s<1$ of DM.

Another possible source of DM in the $a\nu$MSM could be due to primordial black holes (PBHs): CHI features a (quasi) inflection point in the inflaton potential, which has been proposed in [54,55,56,57,58] as a potential trigger for PBH DM production (see Ref. [59] for a review). However, in Section 6, we show that, although this feature is qualitatively present, the $a\nu$MSM is not quantitatively able to account for any fraction of DM in the form of PBHs.

Therefore, the $a\nu$MSM leads to an axion–sterile neutrino DM scenario, which allows us to relate the axion parameters such as the axion decay constant $f_a$ to the sterile neutrino parameters (the masses of these neutral particles and their mixing with the active neutrinos); this provides us with an interesting link between neutrino and axion physics. The allowed parameter space for this combined axion–sterile neutrino DM scenario is identified in Section 7 taking into account the previously discussed bounds.

Finally, in Section 8 we offer our conclusions.

2. The $a\nu$MSM and Generic Observational Bounds

We now give the details of the $a\nu$MSM that are needed for our purposes (see Refs. [1,13] for an introduction to this model). The SM is extended with three sterile neutrinos $N_i$ and the fields of the KSVZ axion model [2,3] (two Weyl fermions $q_1$, $q_2$ neutral under $\mathrm{SU}{\left(2\right)}_{\mathit{L}}\times \mathrm{U}{\left(1\right)}_{\mathit{Y}}$ and a complex scalar A) .

Correspondingly, the SM Lagrangian, ${\mathcal{L}}_{\mathrm{SM}}$, is extended by adding three terms,

$$\mathcal{L}={\mathcal{L}}_{\mathrm{SM}}+{\mathcal{L}}_N+{\mathcal{L}}_{\mathrm{axion}}+{\mathcal{L}}_{\mathrm{gravity}},$$

(1)

which we define in turn. ${\mathcal{L}}_N$ represents the N-dependent piece:

$$i{\overline{N}}_i/\partial N_i+\left(\frac{1}{2}{N}_i{M}_{ij}{N}_j+Y_{ij}{L}_{i}H{N}_j+\mathrm{h}.\mathrm{c}.\right).$$

(2)

We take the Majorana mass matrix M diagonal and real, $M=diag(M_1,M_2,M_3),$ without loss of generality, but the Yukawa matrix Y is generic. ${\mathcal{L}}_{\mathrm{axion}}$ is the KSVZ piece:

$${\mathcal{L}}_{\mathrm{axion}}=i\sum _{j=1}^2{\overline{q}}_j/D{q}_j+{|\partial A|}^2-(y{q}_{2}A{q}_1+h.c.)-\Delta V(H,A),$$

where $\Delta V(H,A)$ is the A-dependent piece of the classical potential

$$\Delta V(H,A)\equiv {\lambda }_A{\left(\right|A|}^2-f_a^2{/2)}^2+{\lambda }_{HA}{\left(\right|H|}^2-v^2{\left)\right(\left|A\right|}^2-f_a^2/2),$$

$v\simeq 174$ GeV is the EW breaking scale and $f_a$ is the axion decay constant. The Yukawa coupling y is chosen real and positive without loss of generality. Finally,

$${\mathcal{L}}_{\mathrm{gravity}}=-\left(\frac{{\overline{M}}_{\mathrm{Pl}}^2}{2}+{\xi }_H{\left(\right|H|}^2-v^2)+{\xi }_A{\left(\right|A|}^2-f_a^2/2)\right)R-\Lambda ,$$

(3)

where ${\overline{M}}_{\mathrm{Pl}}$ is the reduced Planck mass, R is the Ricci scalar, ${\xi }_H$ and ${\xi }_A$ are the non-minimal couplings of H and A to gravity and $\Lambda$ is the cosmological constant. In our model, the inflaton is identified with the Higgs; it is possible to do so with ${\xi }_H\sim \mathcal{O}\left(10\right)$, as discussed in Ref. [13], when we are close to the frontier between the stability and the metastability of the EW vacuum (critical Higgs inflation).

After EW symmetry breaking, the neutrinos acquire a Dirac mass matrix $m_D=vY$, which can be parameterized in terms of column vectors $m_{Di}$ ($i=1,2,3$), i.e., $m_D=\left(m_{D1},m_{D2},m_{D3}\right).$ The active neutrino masses $m_i$ ($i=1,2,3$) are obtained by diagonalizing the matrix

$$m_{\nu }=\frac{m_{D1}{m}_{D1}^T}{M_1}+\frac{m_{D2}{m}_{D2}^T}{M_2}+\frac{m_{D3}{m}_{D3}^T}{M_3}.$$

(4)

We then express Y in terms of the $M_i$ and $m_i$ as done in Refs. [1,13].

On the other hand, the PQ symmetry breaking induced by $\langle A\rangle =f_a/\sqrt{2}$ leads to the quark mass $M_q=y{f}_a/\sqrt{2}$ and the scalar squared mass

$$M_A^2=f_a^2\left(2{\lambda }_A+\mathcal{O}\left(\frac{v^2}{f_a^2}\right)\right).$$

(5)

Since $f_a\gtrsim {10}^8$ GeV (see Ref. [60] for a review), the $\mathcal{O}\left(v^2/f_a^2\right)$ term is very small and will be neglected.

Let us now discuss the other generic observational bounds that are relevant for our purposes3 (with the exception of the bounds related to DM, which will be discussed in the following sections). As far as the active neutrinos are concerned, we have several data from oscillation and non-oscillation experiments. For example, Refs. [61,62] presented some of the most recent determinations of $\Delta m_{21}^2$ and $\Delta m_{3l}^2$ (where $\Delta m_{ij}^2\equiv m_i^2-m_j^2$ and $\Delta m_{3l}^2\equiv \Delta m_{31}^2$ for normal ordering and $\Delta m_{3l}^2\equiv -\Delta m_{32}^2$ for inverted ordering), as well as of the active neutrino mixing angles and the CP phase in the Pontecorvo–Maki–Nakagawa–Sakata (PMNS) matrix. Here, we take the currently most precise values reported in [61,62] for normal ordering (which is currently preferred). Regarding the SM sector, we also have to fix the values of the relevant SM couplings at the EW scale, say, at the top mass $M_t\simeq 172.5$ GeV [63]. We take the values computed in [17], which expresses these quantities in terms of $M_t$, the Higgs mass $M_h\simeq 125.1$ GeV [64], the strong fine-structure constant renormalized at the Z mass, ${\alpha }_s\left(M_Z\right)\simeq 0.1184$ [65] and $M_W\simeq 80.379$ GeV [64] (see the quoted literature for the uncertainties on these quantities).

3. Axion Dark Matter

As we will discuss, axion DM is produced in the $a\nu$MSM by two mechanisms: the misalignment one [26,27,28] and the decay of topological defects [29,30,31,32,33,34,35,36,37]. In order to determine these contributions to DM, the topological susceptibility $\chi$ (given in terms of $m_a$ by $\chi =m_a^2{f}_a^2$) is needed. The axion mass $m_a$ and thus $\chi$ are complicated functions of the temperature T. Here, we use the precise calculations of $\chi$ provided by [66,67].

As discussed in [37], the energy density ${\rho }_a^{\mathrm{mis}}$ due to axions produced by the misalignment mechanism contributes a fraction ${\Omega }_a^{\mathrm{mis}}={\rho }_a^{\mathrm{mis}}/{\rho }_{\mathrm{cr}}$ given by4

$${\Omega }_a^{\mathrm{mis}}{h}^2=(0.12\pm 0.02){\left(\frac{f_a}{1.92\times {10}^{11}\mathrm{GeV}}\right)}^{1.165}.$$

(6)

Requiring that the axion energy density does not exceed the total DM energy density ${\rho }_{\mathrm{DM}}$ (and using ${\Omega }_{DM}\equiv {\rho }_{\mathrm{DM}}/{\rho }_{\mathrm{cr}}=(0.1186\pm 0.0020)/h^2$ [68]), we find the upper bound

$$f_a\lesssim 2\times {10}^{11}\mathrm{GeV}\left(\mathrm{from}\mathrm{misalignment}\right).$$

(7)

Higgs inflation features a high reheating temperature, $T_{\mathrm{RH}}\gtrsim {10}^{13}$ GeV, thanks to the sizable couplings between the Higgs and other SM particles [69,70]. Thus, $T_{\mathrm{RH}}\gg f_a$ and the PQ symmetry is restored after inflation in the $a\nu$MSM.

Therefore, here, axion DM is also produced through decays of topological defects, which leads to a contribution ${\rho }_a^{\mathrm{string}}\equiv {\Omega }_a^{\mathrm{string}}{\rho }_{\mathrm{cr}}$ to the energy density, which in our model is given by [37]

$${\Omega }_a^{\mathrm{string}}{h}^2=0.{37}_{-0.2}^{+0.3}{\left(\frac{f_a}{1.92\times {10}^{11}\mathrm{GeV}}\right)}^{1.165}\frac{ln\left(f_a{t}_{co}\sqrt{{\lambda }_A/\zeta }\right)}{50}.$$

(8)

In Equation (8), the time $t_{co}$ can be determined in terms of $m_a$ by

$$\frac{2\pi {ϵ}_a}{t_{co}}=m_a\left(t_{co}\right)$$

(9)

and the numerical simulations of [71] give ${ϵ}_a=4\pm 0.7$ and $\zeta =1\pm 0.5$.

Therefore, the time $t_{co}$ can be computed by using the precise calculations of $\chi$. However, since ${\Omega }_a^{\mathrm{string}}$ depends on $t_{co}$ only logarithmically, we can, as we explain now, simply estimate $t_{co}$ by using the dilute instanton gas approximation, which gives a power-law temperature dependence of $\chi$,

$$\chi \left(T\right)={\chi }_0{\left(\frac{T_{\mathrm{QCD}}}{T}\right)}^n,$$

(10)

where $T_{\mathrm{QCD}}$ is the temperature of the QCD confining phase transition ($T_{\mathrm{QCD}}\simeq 157$ MeV), $n\simeq 8.16$ and ${\chi }_0\simeq 0.0216{\mathrm{fm}}^{-4}\simeq {\left(75.6\mathrm{MeV}\right)}^4$ [67]. Here, our treatment starts to diverge from that of [37] because the quartic coupling ${\lambda }_A$ does not need to be tiny in our model (unlike in [37]). First, note that in the radiation dominated era the Friedmann equation can be written in the form

$$t^{-2}=\frac{2{\pi }^2}{45{\overline{M}}_{\mathrm{Pl}}^2}{g}_{*}\left(T\right)T^4,$$

(11)

where $g_{*}$ is the effective number of relativistic species. Using this result and Equations (9) and (10), one finds

$$T_{co}={\left(\frac{45{\overline{M}}_{\mathrm{Pl}}^2{\chi }_0{T}_{\mathrm{QCD}}^n}{8{\pi }^4{ϵ}_a^2{f}_a^2{g}_{*}}\right)}^{\frac{1}{4+n}}\simeq \mathrm{GeV}{\left(\frac{2\times {10}^{11}\mathrm{GeV}}{f_a}\right)}^{\frac{2}{4+n}}{g}_{*}^{-\frac{1}{4+n}},$$

(12)

so

$$T_{co}\simeq 0.8\mathrm{GeV}{\left(\frac{2\times {10}^{11}\mathrm{GeV}}{f_a}\right)}^{\frac{2}{4+n}},$$

(13)

where we used well-known determinations of $g_{*}$ in the SM (see, e.g., [67]) and the fact that the contributions of the extra particles beyond the SM to $g_{*}$ are negligible at those temperatures. As a check of this result, note that the power-law temperature dependence of $\chi$ already fits reasonably well the full lattice results from temperatures of the order of a few hundred MeV (see Figure 2 of [67]). Now, using again the Friedmann equation in (11), we have

$$t_{co}\simeq 4\times {10}^{-7}\mathrm{s}{\left(\frac{f_a}{2\times {10}^{11}\mathrm{GeV}}\right)}^{\frac{4}{4+n}}.$$

(14)

We can equivalently use Equation (13) or Equation (14) to estimate the argument of the logarithm in (8) because (9) and (10) tell us

$$f_a{t}_{co}\sqrt{{\lambda }_A/\zeta }=\frac{2\pi {ϵ}_a{f}_a^2\sqrt{{\lambda }_A}}{\sqrt{{\chi }_0\zeta }}{\left(\frac{T_{co}}{T_{\mathrm{QCD}}}\right)}^{\frac{n}{2}}.$$

(15)

Therefore, Equation (13) or Equation (14) allows us to estimate ${\Omega }_a^{\mathrm{string}}$ for each value of $f_a$ and ${\lambda }_A$. So, fixing ${\lambda }_A$, we obtain a maximal value of $f_a$, which we call here $f_a^{\max }$, from the requirement that the axion energy density does not exceed the total DM energy density ${\rho }_{\mathrm{DM}}$, namely

$${\Omega }_a^{\mathrm{mis}}+{\Omega }_a^{\mathrm{string}}\le {\Omega }_{\mathrm{DM}},$$

(16)

where ${\Omega }_{\mathrm{DM}}\equiv {\rho }_{\mathrm{DM}}/{\rho }_{\mathrm{cr}}$. If one considers, for example, ${\lambda }_A=0.1$ the argument of the logarithm in (8) becomes

$$f_a{t}_{co}\sqrt{{\lambda }_A/\zeta }\simeq 4\times {10}^{28}{\left(\frac{f_a}{2\times {10}^{11}\mathrm{GeV}}\right)}^{\frac{8+n}{4+n}}.$$

(17)

and so

$$f_a\le f_a^{\max }\simeq 5\times {10}^{10}\mathrm{GeV}(for{\lambda }_A=0.1),$$

(18)

which is significantly lower than the pure misalignment bound in (7). In Figure 1, we show how $f_a^{\max }$ depends on ${\lambda }_A$.

4. Sterile Neutrino Dark Matter

Another important source of DM in the $a\nu$MSM is the sterile neutrino ${\tilde{N}}_1$ with the smallest mass $m_s$. Generically, this is not exactly $N_1$ due to an active–sterile neutrino mixing. Indeed, in an inflationary5 universe, the production of this particle occurs through its mixing with the active neutrinos of the SM. Such mixing is described by three (generically complex) quantities ${\theta }_{\alpha 1}$, where $\alpha$ represents the flavor of the active neutrino ${\nu }_{\alpha }$ ($\alpha =e,\mu ,\tau$). The ${\theta }_{\alpha 1}$ are the elements of the matrix $\Theta \equiv m_D{M}^{-1}$. It is convenient to introduce a total mixing parameter $\theta$ defined by [42]

$${\theta }^2\equiv \sum _{\alpha =e,\mu ,\tau }{\left|{\theta }_{\alpha 1}\right|}^2.$$

(19)

The mixing with the active neutrinos leads to the production of sterile neutrinos in the early universe in two ways. One is provided by the oscillations between active and sterile states. Furthermore, sterile neutrinos are produced in scatterings. Though this production mechanism is “thermal” in the sense that the ${\tilde{N}}_1$ are produced in scatterings in a thermal plasma, generically these particles are not in thermal equilibrium because of their tiny couplings. As we will discuss in the following subsections, the energy density ${\rho }_s$ of the lightest sterile neutrino can generically account for a non-negligible fraction $X_s$ of the total DM abundance, i.e., ${\Omega }_s=X_s{\Omega }_{\mathrm{DM}}$, where ${\Omega }_s={\rho }_s/{\rho }_{\mathrm{cr}}$ and $0\le X_s\le 1$.

In order to identify the allowed regions of the parameter space, it is necessary to have the observational bounds for an arbitrary value of $X_s$. Of course, the bounds will be generically weaker for $X_s<1$ than for $X_s=1$, but we want to know how they change varying $X_s$.

First, a fermionic DM candidate is subject to a phase-space lower bound on its mass (also known as the Tremaine–Gunn bound [72]) that is related to Pauli’s exclusion principle6. The strongest information comes from the dwarf spheroidal galaxies (dSphs), which are the most compact DM-dominated objects observed to date. In objects of this sort, the dynamics of the DM particles can be characterized by some coarse-grained primordial phase-space density $\mathcal{D}$ and the one-dimensional velocity $\sigma$; see [74] for a detailed discussion. Since the coarse-grained phase-space density either remains constant or diminishes, today we have [74] (see also [75])

$$\frac{{\rho }_s}{{\sigma }^3}\le 3^{3/2}{m}_s^4\mathcal{D}.$$

(20)

Therefore, writing ${\rho }_s=X_s{\rho }_{\mathrm{DM}}$, we obtain the lower mass bound

$$m_s\ge {\left(\frac{X_s{\rho }_{\mathrm{DM}}}{3^{3/2}\mathcal{D}{\sigma }^3}\right)}^{1/4}.$$

(21)

This result tells us that the phase-space bound is rescaled towards smaller values by $X_s^{1/4}\le 1$. A recent publication [76], which we use here, has set the phase-space bound $m_s\gtrsim 190$ eV at the $2\sigma$ level for $X_s=1$. It is worth mentioning that an axion-like particle with mass around the ${10}^{-22}$ eV scale (which is not the QCD axion of the $a\nu$MSM) can also be constrained with phase-space data [77].

Other important bounds on sterile neutrino DM come from the search of X-rays produced by the radiative decay ${\tilde{N}}_1\to \gamma {\nu }_{\alpha }$ [78,79]. Since the differential flux produced by the decay of sterile neutrinos depends on the product $X_s{sin}^2\left(2\theta \right)$, for each chosen value of $m_s$, the upper limit on ${sin}^2\left(2\theta \right)$ must weaken by decreasing $X_s$, rescaling exactly as $1/X_s$ (see also Ref. [80] for a related study). The most recent publications providing this type of X-ray bound used the data collected by the NuSTAR satellite, a space-based X-ray telescope, observing the Milky Way (see [42] for a review) and, more recently, the Andromeda galaxy (M31) [81]. We will take these bounds into account in Section 7.

Now, we describe in turn various production mechanisms for ${\Omega }_s$ as a function of $X_s$.

4.1. Non-Resonant Production

One important production mechanism of sterile neutrinos is the Dodelson–Widrow (DW) mechanism [38], which we will refer to as the non-resonant production. Precise calculations of ${\Omega }_s$ within this mechanism lead to [82,83]

$$m_s\simeq 3.28\times \mathrm{keV}{\left(\frac{{sin}^2\left(2\theta \right)}{{10}^{-8}}\right)}^{-0.615}{\left(\frac{{\Omega }_s}{0.26}\right)}^{0.5}\left\{0.547\times erfc\left[-0.969{\left(\frac{{\mathit{T}}_{\mathrm{QCD}}}{157\mathrm{MeV}}\right)}^{2.15}\right]\right\},$$

(22)

where erfc is the complementary error function. We have used here a normalization such that the argument of the curly bracket in (22) equals 1 for $T_{\mathrm{QCD}}\simeq 157$ MeV.

In addition to the phase-space and X-ray bounds already discussed, sterile neutrino DM is also subject to structure formation bounds. This is because the typical sterile neutrino momentum distribution exhibits a free-streaming length in the early universe, which modifies the formation of structures. These types of bounds are affected by considerable uncertainties related to, among other things, the simulation of non-linear structure formation as well as the difficulty to observe small-scale structures. Furthermore, this structure formation bound depends on the specific sterile neutrino production mechanism one considers. For the non-resonant production, a study for generic $X_s$ has, however, already been performed in [84]. To have an idea of the orders of magnitude, Ref. [84] found a bound on $m_s$ that is about 10 keV for $X_s=1$ and 1 keV for $X_s=0.1$ and so typically stronger than the phase-space bound.

4.2. Resonant Production

The second mechanism to produce sterile neutrino DM is a resonantly enhanced version of the DW mechanism, which relies on a non-vanishing lepton asymmetry L [39,85] and is based on the Mikheyev–Smirnov–Wolfenstein effect [86,87] (see Refs. [88,89,90] for more recent and precise calculations). In practice, the effective mixing in the plasma is enhanced by L, such that the abundance of active neutrinos allows the creation of sterile neutrinos more efficiently. This asymmetry can be generated dynamically by the heavier sterile neutrinos $N_2$ and $N_3$ if their masses are around the GeV scale [91,92,93]. Moreover, $N_2$ and $N_3$ provide a mechanism to generate baryon asymmetry through a different version of leptogenesis [94,95].

The literature to date has focused on the case in which all DM is due to ${\tilde{N}}_1$ (i.e., $X_s=1$), but in our model, the axion also contributes to DM so we need to find more general formulæ that hold for arbitrary $X_s$. In the non-resonant DW mechanism, the quantities ${\Omega }_s$, $\theta$ and $m_s$ are related by Equation (22), which has the form $f({\Omega }_s,\theta ,m_s)=0$, such that for each fixed value of ${\Omega }_s$, the DW mechanism is represented by a line in the $(\theta ,m_s)$ plane. In the resonant production, this function acquires an extra dependence on L, i.e., $f({\Omega }_s,\theta ,m_s,L)=0$, and the allowed region in the $(\theta ,m_s)$ plane is promoted to a band, which is limited by the DW line $f({\Omega }_s,\theta ,m_s,0)=0$. There exists another bound on this band, $f({\Omega }_s,\theta ,m_s,L_{\max })=0$, where $L_{\max }$ is the maximal value of L allowed by observations: the values $L>L_{\max }$ are ruled out because they would excessively change the abundances of light elements produced during Big Bang nucleosynthesis (BBN) [96]. To obtain this bound explicitly for each value of $X_s$, let us observe that the (dimensionless) yield $Y_s\equiv n_s/s$, where $n_s$ is the sterile neutrino density and s is the entropy density, is related to ${\Omega }_s$ through

$${\Omega }_s=\frac{m_s{n}_s}{{\rho }_{\mathrm{cr}}}=\frac{m_s{Y}_s}{{\rho }_{\mathrm{cr}}/s}$$

(23)

so

$$m_s=\frac{X_s{\Omega }_{\mathrm{DM}}{\rho }_{\mathrm{cr}}/s}{Y_s}.$$

(24)

Note that if we approximate $Y_s$ as a function of7 $\theta$ and L, we only reproduce the linear dependence of ${\Omega }_s$ on $m_s$ found in [39]. Using known results of the literature (see Ref. [42] for a review), one obtains, within this approximation

$$m_s\gtrsim \frac{X_s{\Omega }_{\mathrm{DM}}{\rho }_{\mathrm{cr}}/s}{Y_s(\theta ,L_{\max })}.$$

(25)

This formula tells us that the above-mentioned BBN bound in the resonant production band in the $(\theta ,m_s)$ plane is rescaled towards smaller values of $m_s$ by $X_s\le 1$.

5. Sterile Neutrino Dark Matter in a CPT-Symmetric Universe

It was recently pointed out that another mechanism to produce sterile neutrino DM is present if one constructs a CPT-symmetric universe [43,44] in the absence of inflation: the universe before the Big Bang is the CPT reflection of the universe after the Big Bang, so that the time evolution of the universe does not spontaneously violate CPT. In this scenario, a sterile neutrino cosmic abundance is produced according to late-time comoving observers like us just because the vacuum is time dependent. Therefore, unlike the production mechanisms of Section 4, a mixing of the sterile neutrino ${\tilde{N}}_1$ responsible for DM and the active neutrinos are not necessary. One can, therefore, set this mixing to zero, requiring the theory to be invariant under a $Z_2$ symmetry acting on ${\tilde{N}}_1$, which also makes ${\tilde{N}}_1$ exactly stable. As a result, the sterile neutrinos produced through this mechanism can easily avoid the X-ray bounds discussed in Section 4.

In our model, inflation can occur and can be triggered by the Higgs, therefore, we do not perform a general study of this possibility8. However, it is interesting to see how the calculations of [43,44] change in the presence of another DM component, which in our case is due to the axion.

As shown in [44], assuming that this production mechanism occurs in the radiation-dominated era, the yield of the sterile neutrino can be expressed in terms of its mass $m_s$:

$$Y_s=\frac{3I}{2{\pi }^2}{\left(\frac{15}{g_{*}}\right)}^{1/4}{\left(\frac{m_s}{\widehat{\mu }}\right)}^{3/2},$$

(26)

where

$$I\equiv \frac{1}{2{\pi }^2}{\int }_0^{\infty }dx{x}^2\left[1-\sqrt{1-e^{-x^2}}\right]\simeq 0.01276$$

(27)

and $\widehat{\mu }\simeq 5.966\times {10}^{18}$ GeV. In this case, the predicted sterile neutrino contribution to the DM energy density is

$$X_s{\rho }_{\mathrm{DM}}=m_s{n}_s=m_s{Y}_{s}s.$$

(28)

Using the known value of ${\rho }_{\mathrm{DM}}/s$, we find the sterile neutrino mass that is required to account for a fraction $X_s$ of the DM abundance:

$$m_s\simeq 4.8\times {10}^{8}GeV{X}_s^{2/5}{\left(\frac{g_{*}}{g_{*}^{\mathrm{SM}}}\right)}^{1/10},$$

(29)

where $g_{*}^{\mathrm{SM}}$ is the effective number of relativistic degrees of freedom in the SM ($g_{*}^{\mathrm{SM}}\simeq 106.75$ for $T\gg 100$ GeV). We note a dependence on $X_s$ and a (milder) dependence on $g_{*}$. We also observe that the phase-space lower bound discussed in Section 4 is always satisfied down to negligibly small values of $X_s$.

Note that in all the sterile neutrino production mechanisms that we have discussed in this section and Section 4, $m_s$ is generically well below (at least six orders of magnitude) the ${10}^{14}$ GeV scale. Then, from Equation (4), using the observational bounds on $m_i$ and v, it follows that the impact on the RGEs (see Appendix A) of the Yukawa couplings $Y_{ij}$ is generically negligible compared to the other contributions in Appendix A. This is a good thing because the $Y_{ij}$ contribute negatively to the $\beta$-function ${\beta }_{{\lambda }_H}^{\left(1\right)}$ of the Higgs quartic coupling. Moreover, when ${\tilde{N}}_1$ gives a sizeable contribution to DM, the value of $f_a$ required to reproduce the observed DM abundance decreases, as is clear from Section 3, and then so does $M_A$ (cf. Equation (5)). Therefore, the extra scalar starts stabilizing the EW vacuum from smaller energies [1,13,45,46]. It follows that requiring the sterile neutrino to contribute to DM naturally favors EW vacuum stability.

6. Primordial Black Holes as Dark Matter?

PBHs may be generated if the curvature power spectrum $P_{\mathcal{R}}$ has a peak of order ∼${10}^{-2}$ [58], about seven orders of magnitude larger than at ∼60 e-folds before the end of inflation. An enhancement of $P_{\mathcal{R}}$ generically occurs when the inflaton potential features a (quasi) inflection point9 [54,55,56,57,58]. This is the case in CHI [14,15,16,23,54,55], but in order to see if $P_{\mathcal{R}}$ reaches the required order of magnitude in the $a\nu$MSM, a study of this quantity together with other observables is required. We perform this study in this section.

As discussed, e.g., in [23], studying Higgs inflation in the unitary gauge, the potential of the canonically normalized Higgs field ${\varphi }^{\prime }$ is given by

$$U_H\equiv \frac{V_H}{{\Omega }_H^4}=\frac{{\lambda }_H\varphi {\left({\varphi }^{\prime }\right)}^4}{4{(1+{\xi }_H\varphi {\left({\varphi }^{\prime }\right)}^2/{\overline{M}}_{\mathrm{Pl}}^2)}^2},$$

(30)

where $V_H={\lambda }_H{\varphi }^4/4$, $\varphi$ is the Higgs field non-minimally coupled to gravity, which is related to ${\varphi }^{\prime }$ through

$$\frac{d{\varphi }^{\prime }}{d\varphi }={\Omega }_H^{-2}\sqrt{{\Omega }_H^2+\frac{3{\overline{M}}_{\mathrm{Pl}}^2}{2}{\left(\frac{d{\Omega }_H^2}{d\varphi }\right)}^2},$$

(31)

and ${\Omega }_H^2$ is defined by

$${\Omega }_H^2\equiv 1+\frac{2{\xi }_H{\left|H\right|}^2}{{\overline{M}}_{\mathrm{Pl}}^2}.$$

(32)

In a spatially flat Friedmann–Robertson–Walker geometry, the equations for the spatially homogeneous field ${\varphi }^{\prime }\left(t\right)$ and the cosmological scale factor $a\left(t\right)$ are

$${\ddot{\varphi }}^{\prime }+\frac{\sqrt{3{\dot{\varphi }}^{\prime 2}+6U_H}}{\sqrt{2}{\overline{M}}_{\mathrm{Pl}}}{\dot{\varphi }}^{\prime }+\frac{d{U}_H}{d{\varphi }^{\prime }}=0$$

(33)

and

$$H_I^2=\frac{{\dot{\varphi }}^{\prime 2}+2U_H}{6{\overline{M}}_{\mathrm{Pl}}^2},$$

(34)

where a dot represents the derivative with respect to cosmic time t and $H_I\equiv \dot{a}/a$. Inflation in general takes place when10

$$ϵ\equiv -\frac{{\dot{H}}_I}{H_I^2}<1.$$

(35)

Moreover, when

$$\delta \equiv -\frac{{\ddot{\varphi }}^{\prime }}{H_I{\dot{\varphi }}^{\prime }}$$

(36)

is small, one can neglect the inertial term in the inflaton in Equation (33) and reduce the problem to a single first order differential equation, leading to the useful slow-roll approximation where the parameters

$${ϵ}_H\equiv \frac{{\overline{M}}_{\mathrm{Pl}}^2}{2}{\left(\frac{1}{U_H}\frac{d{U}_H}{d{\varphi }^{\prime }}\right)}^2,{\eta }_H\equiv \frac{{\overline{M}}_{\mathrm{Pl}}^2}{U_H}\frac{d^2{U}_H}{d{\varphi }^{\prime 2}}$$

(37)

are small. These slow-roll functions can be constructed through the more general “horizon flow functions” of Ref. [100].

The number of e-folds is defined by

$$N_e\equiv {\int }_{t_b}^{t_e}dt{H}_I\left(t\right),$$

(38)

where $t_e$ is the time at the end of inflation and $t_b$ is the time when the various inflationary observables such as $P_R$, the corresponding spectral index $n_s$ and the tensor-to-scalar ratio r are determined through observations. In the slow-roll approximation $N_e$ is expressed as a function of the field ${\varphi }_b^{\prime }$ (at $t_b$) rather than as a function of time,

$$N_e={\int }_{{\varphi }_e^{\prime }}^{{\varphi }_{\mathrm{b}}^{\prime }}\frac{U_H}{{\overline{M}}_{\mathrm{Pl}}^2}{\left(\frac{d{U}_H}{d{\varphi }^{\prime }}\right)}^{-1}d{\varphi }^{\prime },$$

(39)

where ${\varphi }_e^{\prime }$ is the field value at the end of inflation, and $P_R$ at ${\varphi }_b^{\prime }$ can be computed through

$$P_R=\frac{U_H/{ϵ}_H}{24{\pi }^2{\overline{M}}_{\mathrm{Pl}}^4}.$$

(40)

At the quantum level, these inflationary formulæ remain approximately valid except that one must consider ${\lambda }_H$ and ${\xi }_H$ as functions of ${\varphi }^{\prime }$. In defining these functions, there are well-known ambiguities [5,6,15,101,102]. Here, we adopt the quantization used in [13], which can be embedded in a UV completion of gravity [103,104,105,106,107] (see Refs. [108,109] for reviews). In this approach, the ${\varphi }^{\prime }$-dependence of ${\lambda }_H$ and ${\xi }_H$ is obtained by solving the RGEs given in Appendix A as explained in [13]. The typical shape of the effective inflationary potential (close to criticality) computed in this way is the one shown in11 Figure 2.

We find that the above-mentioned requirement to generate PBHs is never satisfied, so PBHs cannot contribute to DM in our model. The reason is the following. Although $P_{\mathcal{R}}$ does have a peak at a time after inflation as a consequence of the inflection point, its height is several orders of magnitude smaller than ${10}^{-2}$ when one requires a plausible number of e-folds. This situation is illustrated in Figure 3. In that figure, we approach criticality by varying ${\lambda }_{HA}$, but varying other parameters leads to similar situations. We find that the slow-roll approximation is still reasonably good to give at least the order of magnitude of $P_{\mathcal{R}}$ because both $ϵ$ and $\delta$ are well below 1 around the peak of the power spectrum as shown in that figure. Indeed, the numbers of e-folds for ${\lambda }_{HA}$ and ${\lambda }_{HA}-\delta {\lambda }_{HA}$ are $N_e\simeq 65$ and $N_e\simeq 71$, respectively, while the corresponding values computed with the slow-roll approximation are reasonably close ($N_e\simeq 63$ and $N_e\simeq 70$, respectively). Although the height of the peak of $P_{\mathcal{R}}$ does increase by approaching criticality, it does so at the price of increasing $N_e$ above the bound of [110]: for a pretty low peak of order ${10}^{-7}$, the number of e-folds is already starting to be significantly above $\sim 60$. The more we approach criticality, the larger $N_e$ becomes.

In Figure 3, we set the parameters in a way to reproduce the observed neutrino oscillations and have a stable EW vacuum, a viable inflation and baryogenesis through leptogenesis. When all these requirements are satisfied, we always find that PBHs cannot contribute to DM in our model. Essentially, the reason is that the shape of the potential, although apparently able to produce PBH DM, does not have the right quantitative features to do so.

7. Axion–Sterile Neutrino Dark Matter

Having established that the only sources of DM in the $a\nu$MSM are axions and sterile neutrinos, we now identify the allowed parameter space in a combined axion–sterile neutrino DM scenario, taking into account all the previously discussed bounds.

Note that in our model, $X_s$ can then be expressed as

$$X_s=1-\frac{{\Omega }_a^{\mathrm{mis}}+{\Omega }_a^{\mathrm{string}}}{{\Omega }_{\mathrm{DM}}},$$

(41)

which relates $X_s$ and $f_a$. We recall that ${\Omega }_a^{\mathrm{string}}$ also depends on ${\lambda }_A$ (see Equation (8)).

In Figure 4, we show the region corresponding to the resonant sterile neutrino production in the $({sin}^2\left(2\theta \right),m_s)$ plane varying $X_s$. The plot includes the non-resonant production mechanism (the upper line) as a limiting case with vanishing lepton asymmetry (see Section 4.1 and Section 4.2). For each value of $X_s$, we also show the corresponding $f_a$ in two cases. The first case corresponds to a negligible ${\Omega }_a^{\mathrm{string}}$. In the second case, we give the value of the axion decay constant, which we call $f_a^{\mathrm{mis}+\mathrm{string}}$, taking into account both ${\Omega }_a^{\mathrm{mis}}$ and ${\Omega }_a^{\mathrm{string}}$ for ${\lambda }_A=0.1$. This is one of the values for which we cannot only account for the whole DM with axions and sterile neutrinos, but we can also reproduce the observed neutrino oscillation phenomenology, baryon asymmetry, have a stable EW vacuum and critical Higgs inflation (in agreement with Planck observations [111,112]) and solve the strong CP problem [13]. Note that moderate variations of ${\lambda }_A$ around this value produce very small changes in $f_a^{\mathrm{mis}+\mathrm{string}}$ because ${\Omega }_a^{\mathrm{string}}$ depends on ${\lambda }_A$ only logarithmically.

In Figure 5, the sterile neutrino production region of Figure 4 is compared with the X-ray and the phase-space bounds discussed in Section 4. In Figure 6, we also add the structure formation bounds discussed in Section 4.1. As shown in Figure 6, an allowed region for non-resonant sterile neutrino production only appears for $X_s\lesssim 0.3$.

Finally, regarding the CPT-symmetric universe discussed in Section 5, note that, interestingly, we can express $m_s$ as a function of $f_a$ with a (mild logarithmic) dependence on ${\lambda }_A$. This can be obtained by plugging Equation (41) into (29) and using the expressions for ${\Omega }_a^{\mathrm{mis}}$ and ${\Omega }_a^{\mathrm{string}}$ given in Section 3.

8. Conclusions

We have analyzed all DM candidates in the $a\nu$MSM, a simple extension of the SM, originally proposed in [1], which features three sterile neutrinos and the extra fields of the KSVZ QCD axion model. The $a\nu$MSM is well motivated because it not only accounts for DM, neutrino oscillations and baryon asymmetry, but it also solves the strong CP problem, stabilizes the EW vacuum and can implement CHI (in agreement with the most recent Planck observations).

We have ruled out PBHs as a possible source of DM in this model because $P_{\mathcal{R}}$ has a peak that is several orders of magnitude below the required height. Consequently, DM in this model is generically due to the axion and the lightest sterile neutrino. Imposing several constraints, this result allows us to relate the axion parameters such as $f_a$ and ${\lambda }_A$ to the neutrino parameters ($m_s$ and $\theta$).

Requiring the lightest sterile neutrino to contribute to DM in addition to the axion (the only candidate previously considered in the $a\nu$MSM) has several advantages. We have discussed how this requirement generically enlarges the parameter space with absolute EW stability and, as a result, that is where CHI occurs. This inflationary scenario does not suffer from a too low scale of perturbative unitarity breaking and fine tuning of initial conditions (before inflation). On the other hand, the sterile neutrino DM scenario benefits from the presence of an axion DM component because requiring the lightest sterile neutrino to account only for a fraction $X_s<1$ of the DM abundance relaxes all the existing constraints on this scenario. Therefore, one can say that axion and sterile neutrino DM mutually reinforce each other in the $a\nu$MSM.

We plan to keep testing the $a\nu$MSM with future astrophysical and, in particular, cosmological data, such as those regarding the cosmic microwave background and structure formation.