Axion-Sterile-Neutrino Dark Matter

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.


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 Plank scale was presented in [1], where the SM was extended to include three right-handed neutrinos with a generic flavour structure and the extra fields of the simplest invisible QCD axion model, the KSVZ one [2]. The model of [1], which we refer to as the aν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 [3][4][5][6] and solve the strong CP problem through the Peccei-Quinn (PQ) symmetry 1 [8] at the same time.
In Ref. [9] it was found that Higgs inflation can be realized in its critical version [10][11][12] within the aν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 [13]. CHI is particularly interesting for two reasons. One is that it can occur with a moderate, O(10), non-minimal coupling ξ H between the Higgs and the Ricci scalar. Consequently, the scale of breaking of perturbative unitarity, which was noticed in [14], is pushed just below the Planck scale where anyhow new physics is required to UV complete gravity. Furthermore, in Ref. [15] it was shown that CHI, unlike standard Higgs inflation [16], 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 [17].
So far DM in this model has been accounted for exclusively through the axion. However, the aν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ν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 Sec. 2 we briefly review the aνMSM. The gauge group, SU(3) c × SU(2) L × U(1) 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 SU(3) c only. The gravitational sector includes non-minimal couplings of all scalars to gravity, which allow inflation to take place. Sec. 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).
Sec. 3 focuses on the axion contribution to DM. As explained there, we include the contribution from both the misalignment mechanism [18] and the decay of topological defects [19], which have been computed for the KSVZ model in [20]. 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νMSM is determined here explicitly.
Secs. 4 and 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 [21] and resonant [22] production mechanisms, which occur thanks to the mixing between such sterile neutrino and the active neutrinos of the SM (see Refs. [23][24][25] for reviews). The third one takes place in a recently proposed CPT-symmetric universe [26,27], 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ν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 β-function of the Higgs quartic coupling, as explained at the end of Sec. 5. Furthermore, the presence of a sizeable sterile-neutrino contribution to DM, as we will discuss explicitly, allows to reduce the mass of the extra scalar for fixed values of its couplings and so to stabilize the EW vacuum more efficiently [1,9,28,29]. All the sterile-neutrino production mechanisms, therefore, favor EW vacuum stability. This is another motivation for realising a fraction of DM through sterile neutrinos in the aν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νMSM, compared to the case 2 X s = 1 (where such region is quite narrow [35,36]): 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νMSM could be due to primordial black holes (PBHs): CHI features a (quasi) inflection point in the inflaton potential, which has been proposed in [37][38][39][40][41] as a potential trigger for PBH DM production (see Ref. [42] for a review). However, in Sec. 6 we show that, although this feature is qualitatively present, the aνMSM is not quantitatively able to account for any fraction of DM in the form of PBHs.
Therefore, the aν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 Sec. 7 taking into account the previously discussed bounds.
Finally, in Sec. 8 we offer our conclusions.

The aνMSM and generic observational bounds
We now give the details of the aνMSM that are needed for our purposes (see Refs. [1,9] 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] (two Weyl fermions q 1 , q 2 neutral under SU(2) L × U(1) Y and a complex scalar A) . Correspondingly, the SM Lagrangian, L SM , is extended by adding three terms, which we define in turn. L N represents the N -dependent piece: 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. L axion is the KSVZ piece: where ∆V (H, A) is the A-dependent piece of the classical potential 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, whereM Pl is the reduced Planck mass, R is the Ricci scalar, ξ H and ξ A are the non-minimal couplings of H and A to gravity and Λ is the cosmological constant. In our model the inflaton is identified with the Higgs; it is possible to do so with ξ H ∼ O(10), as discussed in Ref. [9], 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 = m D1 , m D2 , m D3 . The active-neutrino masses m i (i = 1, 2, 3) are obtained by diagonalizing the matrix We then express Y in terms of the M i and m i as done in Refs. [1,9]. On the other hand, the PQ symmetry breaking induced by A = f a / √ 2 leads to the quark mass M q = yf a / √ 2 and the scalar squared mass (2.5) Since f a 10 8 GeV (see Ref. [43] for a review), the O (v 2 /f 2 a ) term is very small and will be neglected.
Let us now discuss the other generic observational bounds that are relevant for our purposes 3 (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. [44,45] presented some of the most recent determinations of ∆m 2 21 and ∆m 2 3l (where ∆m 2 ij ≡ m 2 i − m 2 j and ∆m 2 3l ≡ ∆m 2 31 for normal ordering and ∆m 2 3l ≡ −∆m 2 32 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 [44,45] 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 172.5 GeV [46]. We take the values computed in [13], which expresses these quantities in terms of M t , the Higgs mass M h 125.1 GeV [47], the strong fine-structure constant renormalized at the Z mass, α s (M Z ) 0.1184 [48] and M W 80.379 GeV [47] (see the quoted literature for the uncertainties on these quantities).

Axion dark matter
As we will discuss, axion DM is produced in the aνMSM by two mechanisms: the misalignment one [18] and the decay of topological defects [19,20]. In order to determine these contributions to DM the topological susceptibility χ (given in terms of m a by χ = m 2 a f 2 a ) is needed. The axion mass m a and thus χ are complicated functions of the temperature T . Here we use the precise calculations of χ provided by [49,50].
As discussed in [20], the energy density ρ mis a due to axions produced by the misalignment mechanism contributes a fraction Ω mis a = ρ mis a /ρ cr given by 4 Requiring that the axion energy density does not exceed the total DM energy density ρ DM (and using Higgs inflation features a high reheating temperature, T RH 10 13 GeV, thanks to the sizable couplings between the Higgs and other SM particles [52,53]. Thus T RH f a and the PQ symmetry is restored after inflation in the aνMSM.
Therefore, here axion DM is also produced through decays of topological defects, which leads to a contribution ρ string a ≡ Ω string a ρ cr to the energy density, which in our model is given by [20] Ω string and the numerical simulations of [54] give a = 4 ± 0.7 and ζ = 1 ± 0.5. Therefore, the time t co can be computed by using the precise calculations of χ. However, since Ω string a 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 χ, where T QCD is the temperature of the QCD confining phase transition (T QCD 157 MeV), n 8. 16 and χ 0 0.0216 fm −4 (75.6 MeV) 4 [50]. Here our treatment starts to diverge from that of [20] because the quartic coupling λ A does not need to be tiny in our model (unlike in [20]). First, note that in the radiation dominated era the Friedmann equation can be written in the form where g * is the effective number of relativistic species. Using this result and Eqs. (3.4) and (3.5) one finds so where we used well-known determinations of g * in the SM (see e.g. [50]) 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 χ fits reasonably well the full lattice results already from temperatures of order of few hundreds of MeV (see Fig. 2 of [50]). Now, using again the Friedmann equation in (3.6) we have (3.9) We can equivalently use Eq. (3.8) or Eq. (3.9) to estimate the argument of the logarithm in (3.3) because (3.4) and (3.5) tell us (3.10) Therefore, Eq. (3.8) or Eq. (3.9) allows us to estimate Ω string a for each value of f a and λ A . So fixing λ A we obtain a maximal value of f a , which we call here f max a , from the requirement that the axion energy density does not exceed the total DM energy density ρ DM , namely where Ω DM ≡ ρ DM /ρ cr . If one considers for example λ A = 0.1 the argument of the logarithm in (3.3) becomes (3.12) and so which is significantly lower than the pure misalignment bound in (3.2). In Fig. 1 we show how f max a depends on λ A .

Sterile-neutrino dark matter
Another important source of DM in the aνMSM is the sterile neutrinoÑ 1 with the smallest mass m s . Generically, this is not exactly N 1 due to an active-sterile neutrino mixing. Indeed, in an inflationary 5 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 θ α1 , where α represents the flavour of the active neutrino ν α (α = e, µ, τ ). The θ α1 are the elements of the matrix Θ ≡ m D M −1 . It is convenient to introduce a total mixing parameter θ defined by [25] 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Ñ 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 ρ s of the lightest sterile neutrino can generically account for a non-negligible fraction X s of the total DM abundance, i.e. Ω s = X s Ω DM , where Ω s = ρ s /ρ cr and 0 ≤ X s ≤ 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 (a.k.a the Tremaine-Gunn bound [55]) that is related to Pauli's exclusion principle 6 . The strongest information comes from the dwarf spheroidal galaxies (dSphs), which are the most compact DMdominated objects observed so far. In objects of this sort the dynamics of the DM particles can be characterized by some coarse-grained primordial phase-space density D and the one-dimensional velocity σ: see [57] for a detailed discussion. Since the coarse-grained phase-space density either remains constant or diminishes, today we have [57] (see also [58]) Therefore, writing ρ s = X s ρ DM we obtain the lower mass bound This result tells us that the phase-space bound is rescaled towards smaller values by X 1/4 s ≤ 1. A recent publication [59], which we use here, has set the phase-space bound m s 190 eV at 2σ 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νMSM) can also be constrained with phase space data [60].
Other important bounds on sterile neutrino DM come from the search of X-rays produced by the radiative decayÑ 1 → γν α [61,62]. Since the differential flux produced by the decay of sterile neutrinos depends on the product X s sin 2 (2θ), for each chosen value of m s , the upper limit on sin 2 (2θ) must weaken by decreasing X s , rescaling exactly as 1/X s (see also Ref. [63] for a related study). The most recent publications providing this type of X-ray bounds used the data collected by the NuSTAR satellite, a space-based X-ray telescope, observing the Milky Way (see [25] for a review) and, more recently, the Andromeda galaxy (M31) [64]. We will take these bounds into account in Sec. 7. Now we describe in turn various production mechanisms for Ω s as a function of X s .

Non-resonant production
One important production mechanism of sterile neutrino is the Dodelson-Widrow (DW) mechanism [21], which we will refer to as the non-resonant production. Precise calculations of Ω s within this mechanism lead to [65,66] m s 3.28 × keV sin 2 (2θ) where erfc is the complementary error function. We have used here a normalization such that the argument of the curly bracket in (4.4) equals 1 for T QCD 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 exhibit a free-streaming length in the early universe, which modifies the formation of structures. This type 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 [68]. To have an idea of the orders of magnitude, Ref. [68] 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.

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 [22,67] and is based on the Mikheyev-Smirnov-Wolfenstein effect [69] (see Ref. [70][71][72] for more recent and precise calculations). In practice the effective mixing in the plasma gets enhanced by L, such that the abundance of active neutrinos allows to create 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 [73][74][75]. Moreover, N 2 and N 3 provide a mechanism to generate baryon asymmetry through a different version of leptogenesis [76,77].
The literature so far focused on the case in which all DM is due toÑ 1 (i.e. X s = 1), but in our model the axion also contributes to DM so we need to find more general formulae that hold for arbitrary X s . In the non-resonant DW mechanism the quantities Ω s , θ and m s are related by Eq. (4.4), which has the form f (Ω s , θ, m s ) = 0, such that for each fixed value of Ω s the DW mechanism is represented by a line in the (θ, m s ) plane. In the resonant production this function acquires an extra dependence on L, i.e. f (Ω s , θ, m s , L) = 0, and the allowed region in the (θ, m s ) plane is promoted to a band, which is limited by the DW line f (Ω s , θ, m s , 0) = 0. There exists another bound on this band, f (Ω s , θ, 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) [78]. To obtain this bound explicitly for each value of X s let us observe that the (dimensionless) yield Y s ≡ n s /s, where n s is the sterile neutrino density and s is the entropy density, is related to Ω s through Note that if we approximate Y s as a function of 7 θ and L only we reproduce the linear dependence of Ω s on m s found in [22]. Using known results of the literature (see Ref. [25] for a review) one obtains, within this approximation This formula tells us that the above-mentioned BBN bound in the resonant production band in the (θ, m s ) plane is rescaled towards smaller values of m s by X s ≤ 1.

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 [26,27] 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 Sec. 4, a mixing of the sterile neutrinoÑ 1 responsible for DM and the active neutrinos is not necessary. One can, therefore, set this mixing to zero requiring the theory to be invariant under a Z 2 symmetry acting onÑ 1 , which also makesÑ 1 exactly stable. As a result, the sterile neutrinos produced through this mechanism can easily avoid the X-ray bounds discussed in Sec. 4.
In our model inflation can occur and can be triggered by the Higgs, therefore, we do not perform a general study of this possibility 8 . However, it is interesting to see how the calculations of [26,27] change in the presence of another DM component, which in our case is due to the axion.
As shown in [27], 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 : andμ 5.966×10 18 GeV. In this case the predicted sterile-neutrino contribution to the DM energy density is Using the known value of ρ DM /s we find the sterile-neutrino mass that is required to account for a fraction X s of the DM abundance: where g SM * is the effective number of relativistic degrees of freedom in the SM (g SM * 106.75 for T 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 Sec. 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 Sec. 4 m s is generically well below (at least six orders of magnitude) the 10 14 GeV scale. Then from Eq. (2.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 β-function β (1) λ H of the Higgs quartic coupling. Moreover, whenÑ 1 gives a sizeable contribution to DM the value of f a required to reproduce the observed DM abundance decreases, as clear from Sec. 3, then so does M A (cf. Eq. (2.5)). Therefore, the extra scalar starts stabilizing the EW vacuum from smaller energies [1,9,28,29]. It follows that requiring the sterile neutrino to contribute to DM naturally favors EW vacuum stability. enhancement of P R generically occurs when the inflaton potential features a (quasi) inflection point 9 [37][38][39][40][41]. This is the case in CHI [10-12, 15, 37, 38], but in order to see if P R reaches the required order of magnitude in the aνMSM a study of this quantity together with other observables is required. We perform such study in this section.
As discussed e.g. in [15], studying Higgs inflation in the unitary gauge, the potential of the canonically normalized Higgs field φ is given by where V H = λ H φ 4 /4, φ is the Higgs field non-minimally coupled to gravity, which is related to φ through and Ω 2 H is defined by In a spatially flat Friedmann-Robertson-Walker geometry the equations for the spatially homogeneous field φ (t) and the cosmological scale factor a(t) arë and where a dot represents the derivative with respect to cosmic time t and H I ≡ȧ/a. Inflation in general takes place when 10 Moreover, when is small one can neglect the inertial term in the inflaton equation (6.4) and reduce the problem to a single first order differential equation, leading to the useful slow-roll approximation where the parameters are small. These slow-roll functions can be constructed through the more general "horizon flow functions" of Ref. [82]. 9 See also Refs. [80,81] for earlier works on PBH production in inflationary models. 10 As usual, the expansion of the universe is nearly exponential for < 1 and becomes exactly exponential as → 0. The number of e-folds is defined by 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 φ b (at t b ) rather than as a function of time, where φ e is the field value at the end of inflation, and P R at φ b can be computed through Pl .

(6.11)
At quantum level these inflationary formulae remain approximately valid except that one must consider λ H and ξ H as functions of φ . In defining these functions there are well-known ambiguities [4,5,11,83,84]. Here we adopt the quantization used in [9], which can be embedded in a UV completion of gravity [85][86][87][88][89] (see Refs. [90,91] for reviews). In this approach the φ -dependence of λ H and ξ H is obtained by solving the RGEs given in Appendix A as explained in [9]. The typical shape of the effective inflationary potential (close to criticality) computed in this way is the one shown in 11   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 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 Fig. 3. In that figure we approach criticality by varying λ 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 R because both and δ are well-below 1 around the peak of the power spectrum as shown in that figure. Indeed, the number of e-folds for λ HA and λ HA −δλ HA are N e 65 and N e 71, respectively, while the corresponding values computed with the slow-roll approximation is reasonably close (N e 63 and N e 70, respectively). Although the height of the peak of P R does increase by approaching criticality, it does so at the price of increasing N e above the bound of [92]: already for a pretty low peak of order 10 −7 the number of e-folds is starting to be significantly above ∼ 60. The more we approach criticality the larger N e becomes.
In Fig. 3 we set the parameters in a way to reproduce the observed neutrino oscillations, 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, it does not have the right quantitative features to do so. Note that in our model X s can then be expressed as which relates X s and f a . We recall that Ω string a also depends on λ A (see Eq. (3.3)). In Fig. 4 we show the region corresponding to the resonant sterile-neutrino production in the (sin 2 (2θ), 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 Secs. 4.1 and 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 Ω string a . In the second case we give the value of the axion decay constant, which we call f mis+string a there, taking into account both Ω mis a and Ω string a for λ A = 0.1. This is one of the values for which we can not only account for the whole DM with axions and sterile neutrinos, but we can also reproduce the observed neutrino oscillations phenomenology, baryon asymmetry, have a stable EW vacuum, critical Higgs inflation (in agreement with Planck observations [93]) and solve the strong CP problem [9]. Note that moderate variations of λ A around this value produce very small changes in f mis+string a because Ω string a depends on λ A only logarithmically. In Fig. 5 the sterile-neutrino production region of Fig. 4 is compared with the X-ray and the phase-space bounds discussed in Sec. 4. In Fig. 6 we also add the structure-formation bounds discussed in Sec. 4.1. As shown in Fig. 6 an allowed region for non-resonant sterile-neutrino production only appears for X s 0.3.
Finally, regarding the CPT-symmetric universe discussed in Sec. 5, note that, interestingly, we can express m s as a function of f a with a (mild logarithmic) dependence on λ A . This can be obtained by plugging Eq.

Conclusions
We have analysed all DM candidates in the aν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ν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 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 λ A to the neutrino parameters (m s and θ).
Requiring the lightest sterile neutrino to contribute to DM in addition to the axion (the only candidate previously considered in the aνMSM) has several advantages. We have discussed how this requirement generically enlarges the parameter space with absolute EW stability and, as a result, that 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νMSM.
We plan to keep testing the aνMSM with future astrophysical and, in particular, cosmological data, such as those regarding the cosmic microwave background and structure formation. Acknowledgments I thank G. Ballesteros and A. Urbano for useful discussions and J. Rubio for useful mail communications on primordial black holes in related models.

A Renormalization-group equations
For a generic coupling g defined in the MS renormalization scheme we write the RGEs as where d/dτ ≡μ 2 d/dμ 2 andμ is the MS renormalization energy scale. The β-functions β g can also be expanded in loops:  where g 3 , g 2 and g 1 = 5/3g Y are the gauge couplings of SM gauge group SU(3) c , SU(2) L and U(1) Y , respectively, y t is the top Yukawa coupling and λ H is the Higgs quartic coupling appearing in the term λ H (|H| 2 − v 2 ) 2 of the classical potential.
Since the SM couplings evolve in the full range from the EW to the Planck scale it is appropriate to use for them the 2-loop RGEs 12 , which, including the new physics contribution, read: Here we have corrected a missprint of the RGEs provided in Ref [9]: there are no g 2 3 y 2 and y 4 contributions to the RGE of λ H .
The matching at the mass thresholds due to the new scalar A and fermions N i , q 1 and q 2 is performed as explained in Ref. [9].