Sub-GeV dark matter as pseudo-Goldstone from the seesaw scale

Pseudo Nambu-Goldstone bosons (pNGBs) are naturally light spin-zero particles, which can be interesting dark matter (DM) candidates. We study the phenomenology of a pNGB theta associated with an approximate symmetry of the neutrino seesaw sector. A small coupling of theta to the Higgs boson is induced radiatively by the neutrino Yukawa couplings. By virtue of this Higgs portal interaction (i) the pNGB acquires a mass m_theta proportional to the electroweak scale, and (ii) the observed DM relic density can be generated by the freeze-in of theta-particles with mass m_theta ~ 3 MeV. Alternatively, the coupling of theta to heavy sterile neutrinos can account for the DM relic density, in the 1 keV-3 MeV window. The decays of theta into light fermions are suppressed by the seesaw scale, making such pNGB sufficiently stable to play the role of DM.


Introduction
The astrophysical and cosmological evidence for dark matter (DM) requires that a new particle is introduced, with constrained mass and interactions. These interactions must in particular maintain the DM candidate stable on cosmological time scales and generate the observed DM relic abundance. Many candidates can satisfy all constraints. Much fewer, though, provide a convincing explanation for the required values of the DM parameters. As examples, in most models (i) an ad-hoc symmetry is postulated to insure stability, and/or (ii) the DM mass is assumed to be of order (or below) the electroweak (EW) scale as demanded by the Weakly Interacting Massive Particle paradigm, and/or (iii) new interactions with strength similar to the EW ones, not motivated by other means, are just postulated. A scenario which is potentially predictive, and could provide motivations for some of these assumptions, consists in identifying the DM candidate with the pseudo-Nambu-Goldstone boson (pNGB) of a spontaneously broken global symmetry.
An interesting feature of such scenario is that the DM couplings to the Standard Model (SM) fields, which could render the DM candidate unstable on cosmological times, are suppressed by powers of the spontaneous symmetry breaking (SSB) scale. As a result, provided this scale is large, the decay width of the particle can be sufficiently small. A well-known example is the axion, the pNGB associated with the Peccei-Quinn U (1) symmetry, broken spontaneously at the scale f P Q . The axion is a viable DM candidate when f P Q lies in the interval ∼ (10 12 − 10 14 ) GeV. Another possibility is to consider SSB at the seesaw scale, which should be large to account for the smallness of the neutrino masses. One global symmetry related to the seesaw is U (1) B−L . In the case that it is spontaneously broken by an EW singlet, coupled to heavy sterile (right-handed) neutrinos, the corresponding NGB is known as the singlet Majoron [1]. Models with a pseudo-Majoron as DM have been studied [2]- [7]. More generally, we will consider any (family-dependent) global symmetry whose spontaneous breaking contributes to sterile neutrino masses. Symmetries associated with the seesaw sector have the advantage that the pNGB interactions are anticipated to exist anyway to account for neutrino masses.
Another remarkable property of the pNGB scenario, relevant for the DM phenomenology, is that the mass and the scalar potential interactions of a pNGB vanish in the limit of exact symmetry. Therefore, when the symmetry is explicitly broken by a unique (or dominant) term, the scalar interactions and the DM mass are necessarily related, and, if the relic density is essentially determined by such interactions, this implies a one-to-one correspondence between the DM relic density and mass. For instance we will consider a framework where the source of explicit breaking generates radiatively a 'Higgs portal' interaction of the form λθ 2 H † H/2, with θ the pNGB particle and H the SM Higgs doublet. This coupling generates a mass m 2 θ = λv 2 with v = 174 GeV and, at the same time, a relic density from the θ − H interactions. We will show how, given the known thermal distribution of SM particles in the early universe thermal bath, the Higgs portal can generate a DM relic density through the freeze-in mechanism (see e.g. Ref. [8]). The observed relic density is obtained for a unique possible value of the DM mass, which turns out to be m θ 2.8 MeV for a Higgs boson mass of 120 GeV. There is also a second possibility to generate the correct relic density, by the freeze-in of the interaction of θ with the sterile neutrinos. These production mechanisms differ from those previously considered for a pNGB candidate for DM.
A third interesting feature of the pNGB setup is that the DM mass is related to the scale of explicit symmetry breaking, which can be identified with a physical mass scale already present in the theory, analogously to the QCD scale for the axion. For instance, in the scenario we will consider below, the neutrino Dirac Yukawa couplings provide the source of explicit breaking of the global symmetry. Therefore, the DM mass is proportional to the EW scale, m 2 θ = λv 2 , and this provides a justification for the presence of a scalar DM particle at or below the EW scale. Moreover, the coupling λ can be computed as a function of the seesaw couplings. In turn, we will show that in this framework the DM mass can be protected from quadratically divergent corrections. The mechanism we invoke to remove these corrections is a variation of the one proposed long ago by Hill and Ross for very light pNGBs with both scalar and pseudo-scalar couplings [9]. We will consider a U (1) X symmetry such that its explicit breaking requires several Yukawa couplings at the same time, thus lowering the degree of divergence of the contributions to the pNGB effective potential. This mechanism has similarities with the collective breaking mechanism introduced to protect the EW scale in little Higgs models, where the SM Higgs is the pNGB (see Ref. [10] for a review).
All in all, we will show that the seesaw interactions (i) are associated with global symmetries broken spontaneously at the heavy neutrino mass scale, whose largeness guarantees the stability of the pNGB DM candidate, and (ii) have a source of explicit symmetry breaking built-in, and therefore induce the mass and scalar interactions of the DM. Moreover, they can lead to (iii) a one-to-one correspondence between the DM mass and the DM relic density, with (iv) a justification of the presence of DM at low scale, and with (v) a DM mass protected from large radiative corrections.
The paper is organized as follows. In section 2 we display the relevant effective interactions of our pNGB candidate for DM. In section 3 we show how a pNGB relic density of the order of the observed DM density can be obtained through freeze-in. In section 4 we present a class of seesaw models where the required pNGB interactions are generated naturally. In section 5 we compute the pNGB couplings to SM particles and analyze the corresponding constraints on the DM stability. We conclude in section 6.

Effective interactions of the pNGB
In this section we introduce the effective lagrangian of a pNGB θ, associated with a global symmetry of the neutrino sector, broken spontaneously at the seesaw scale f . We postpone to section 4 the detailed description of a class of models where such a light pNGB may emerge naturally. For the present purposes, it suffices to specify the effective interactions of θ at low energy (below f ).
In the class of models under consideration, θ couples to a sterile neutrino ν c as follows: where g is a Yukawa coupling and the 4-component Majorana neutrino is defined as usual by N ≡ (ν c ν c † ) T . For simplicity, we consider only one sterile neutrino, whose mass m N = where Λ is some cutoff scale at or above f . The light neutrino mass scale is given by the standard seesaw relation, m ν = y 2 v 2 /m N . It is useful to use Eq. (3) to express the coupling g as a function of the relevant energy scales in the theory: We will take the factor k ≡ log(Λ 2 /m 2 N )/(8π 2 ) to be of order one. Then, for a fixed value of m ν , m N and also f = √ 2m N /g are determined as a function of m θ and g, which are basically the only two free parameters. We require the scale f to be smaller than the Planck scale, M P = 1. 22 × 10 19 GeV. In addition, we assume for simplicity that f is larger than one TeV, in order for the effective theory to contain only the SM particles, the pNGB and the sterile neutrino, while the degrees of freedom involved in the SSB are decoupled.
These two requirements exclude the green shaded regions in the (m θ − g) plane shown in Fig. 1. The dotted lines in Fig. 1 correspond to constant values of the sterile neutrino mass, m N = m 2 θ /(g 2 m ν k) = 10 2 , 10 6 , 10 10 , 10 14 GeV, from top to bottom.

Relic density from freeze-in of the pNGB
We now study the role of the interactions in Eqs. (1) and (2) for the production of a θ relic density, which can play the role of DM. In this section, we implicitly assume that θ is stable on cosmological time scales. The region of parameters where this stability is achieved will be analyzed in detail in section 5. There are several mechanisms to produce the DM relic density. We begin by briefly recalling these various possibilities and confront them with the parameter space of our scenario. A well-known way to produce a cosmological density of a DM particle species is the usual freeze-out: it requires an annihilation cross section large enough to overcome the expansion of the universe and thermalize the DM particle. The interactions of θ come from the terms in Eq. (1), which we refer to as the 'sterile neutrino portal', and the terms in Eq. (2), that is, the 'Higgs portal'. These interactions may or may not lead to the thermalization of θ in the early universe, depending on whether the interaction rate Γ gets larger than the Hubble rate H. If it does, later on θ necessarily decouples from the plasma and, if stable, acts as DM with a certain energy density ρ θ , which depends in general on its mass and annihilation cross section.
Before determining ρ θ numerically, it is instructive to identify the thermalization region in first approximation, by evaluating the θ interaction rate Γ and requiring Γ > H. We start with the Higgs portal, considering for simplicity the h → θθ decay only, because in our scenario decays turn out to dominate with respect to scattering processes (see Fig. 2 below). The decay rate is This should be compared to the Hubble parameter in the radiation epoch, which is given by H = 1.66 g ρ * T 2 /M P , with g ρ * = 106.75 the degrees of freedom corresponding to the SM content (we shall assume such value in the rest of the paper). Taking T m h to be the relevant temperature, we find that, for a coupling the decays and inverse-decays are able to produce a thermal population of θ-particles. In other words, given the m 2 θ = λv 2 relation, one needs m θ > 44 MeV, which translates into the thin vertical line in Fig. 1. 1 1 Notice that Eq. (6) is the condition for DM to thermalize with the SM sector at T m h (where the rates reach the maximum, see Fig. 2 below). This condition is different from the one usually considered for the freeze-out scenario, where one requires that the scattering processes are in thermal equilibrium at T m θ . The latter condition leads to a value of λ about two orders of magnitude larger than the one in Eq. (6).
Similarly for the interaction with the sterile neutrino N , the process N N → θθ has a rate Γ = σv rel n N g 4 256π with n N the sterile neutrino number density, and we used T m N . This leads to thermal production of θ-particles provided where we made use of Eq. (4). This gives the thin diagonal line in Fig. 1. We do not discuss the effects of the couplings of θ to light fermions, because in the models under consideration they turn out to be far too small to thermalize θ (see section 5).
Outside the thermalization region, where Γ is smaller than H, one may still have a final Ω θ matching Ω DM (as usual we define Ω i ≡ ρ i /ρ crit , i.e. the ratio between the energy density ρ i and the critical density of the universe). The idea is that in the early universe a population of θ-particles is produced through pair production scattering processes, but with a density smaller than the thermal one. For appropriate values of m θ one may still reach a relic density equal to Ω DM . This process, called freeze-in, actually happens in our scenario as we shall describe below in detail (for generic properties of the freeze-in mechanism see Ref. [8]).
In the rest of this section we first introduce the ingredients to make a numerical calculation of the θ relic density, and we then apply them to the freeze-out and freeze-in processes. Afterwards we also discuss some non-thermal production mechanisms and end up with a summary of our results.

Boltzmann equation and reaction densities
From the Higgs portal, θ-particles are created from the following scattering processes: 2 σ(W W → θθ) = 1 9 σ(ZZ → θθ) = 1 9 where Γ h is the total Higgs decay width and n c is the number of colours of the fermion f . The pNGB coupling to the sterile neutrino, g, versus the pNGB mass, m θ , in GeV. We fixed m ν = 0.05 eV and k = 1. The upper (lower) green shaded region is excluded by requiring the SSB scale f to be larger than one TeV (smaller than the Planck scale). The four dotted lines correspond to constant values of the sterile neutrino mass: from top to bottom, m N = 10 2 , 10 6 , 10 10 , 10 14 GeV. The thin diagonal (vertical) line indicates the lower value of g (m θ ) that leads to a thermalization of θ. The thick line corresponds to a θ relic density equal to the observed DM relic density, as follows from the numerical solution of the Boltzmann equation (for m h = 120 GeV). Since θ is produced by freeze-in, its relic density grows with its couplings g and λ = m 2 θ /v 2 , therefore the region below and to the left (above and to the right) of the thick line corresponds to Ω θ < Ω DM (> Ω DM ).
Similarly, θ-particles can be created through the scattering with a sterile neutrino, where we neglected m θ in front of m N and √ s and we defined β N = (1 − 4m 2 N /s) 1/2 . In the following it will be convenient to consider separately the effect of the Higgs decay, Eq. (5), which means that, in order to avoid double counting, the on-shell part of the scattering processes must be subtracted (whenever there is one). In the narrow width approximation this means to perform in Eqs. (10)-(12) the substitution for i = W, Z, f . In the following, all scattering quantities we will consider are meant to be the subtracted ones, using Eq. (14). To calculate the relic density that one obtains from these processes, either through freeze-in or freeze-out, one has to integrate the Boltzmann equation, with z ≡ m h /T conventionally taken as the evolution parameter. Here H(z) is the Hubble parameter, s(z) the entropy density, and Y θ ≡ n θ /s with n θ the number density. The reaction density γ D h contains the effect of the Higgs boson decay, while γ annih is the sum of the (subtracted) reaction densities of the scattering processes in Eqs. (9)- (13). The reaction densities are given by γ(a b ↔ 1 2) = N θ dp a dp b f eq a f eq b dp 1 dp 2 (2π) 4 involving the Bessel functions K 1,2 . We have defined dp ≡ d 3 p/((2π) 3 2E). Here N θ = 2 is the number of θ particles produced per decay or annihilation, f eq i = (e E i /T ± 1) −1 e −E i /T is the Maxwell-Boltzmann energy distribution, |M| 2 is the amplitude squared summed over initial and final spins (with no averaging), s min = max[(m a + m b ) 2 , (m 1 + m 2 ) 2 ], and the reduced cross section is defined bŷ with σ the particle physics cross section of Eqs. (9)-(13), g a,b the number of degrees of freedom of the particles a, b and c ab a combinatorial factor equal to 2 (1) if a and b are identical (different).

Freeze-out
As it is well-known, if θ freezes out relativistically, that is at T m θ , its relic density is independent of the annihilation cross section, because it decouples when the thermal number density, Y eq θ ≡ n eq θ /s is still independent of the temperature. The formula for the relative density Ω θ is simply with g s * the number of degrees of freedom contributing to entropy at θ decoupling. This should match the value Ω DM h 2 0 = 0.11 ± 0.01 .
In these equations, h 0 is the reduced Hubble constant h 0 = H 0 /(100 Mpc km s −1 ) 0.70. As a result, the observed relic density can be obtained only for one DM mass value: m θ 0.15 keV, where we took g s * = 106.75, valid for a decoupling temperature above 100 GeV. Such a light θ can be only thermalized by the interaction with sterile neutrinos, which decouple at a temperature just below m N , leading to relativistic θ's. This result corresponds to the left vertical branch of the thick curve in Fig. 1, which is obtained by integrating numerically the Boltzmann equation.
If instead θ is heavier and thermalizes, to have a small enough relic density requires that it decouples non-relativistically. This can occur through the Higgs portal interaction. In this case the relic density essentially only depends on the annihilation cross section which must take the value σv rel 10 −26 cm 3 s −1 . To get a sufficiently small relic density, one needs to go to much higher values of m θ , where the coupling λ = m 2 θ /v 2 becomes large enough. The observed relic density is obtained for a unique value of m θ , because of the relation between m θ and λ, which is a consequence of the pNGB nature of θ. The freeze-out through the Higgs portal was considered in Ref. [12], assuming this restrictive relation, and it was found that the DM relic density is obtained with m θ = 50 − 70 GeV (for m h = 120 − 180 GeV).
However, as we will see in section 5, the freeze-out values, m θ = 0.15 keV as well as m θ = 50 − 70 GeV, require values of parameters which are excluded by the instability of θ on cosmological time scales. Therefore in our scenario both relativistic and non-relativistic freeze-out is not a viable way to produce DM. Instead, the freeze-in process turns out to be efficient in an allowed region of the parameter space, as we now discuss.

Freeze-in
In the freeze-in process θ is produced by the annihilation or decay of a heavier particle X, with rates small enough not to thermalize, leading to a less-than-thermal DM population, which reaches a plateau at T ∼ m X because at smaller T the number of X-particles is Boltzmann suppressed. If the X-particles are the SM particles, as it is the case through the Higgs portal, their number densities are known as a function of the temperature (simply given by their thermal distribution down to a temperature well below their mass). As a consequence, the created relic density depends only on the magnitude of the portal. 3 Let us first consider a value of the coupling g sufficiently small to neglect θ production at the seesaw scale (we will see that this is the case for g 10 −3 − 10 −4 , depending on m θ ). In this case only the Higgs portal interactions can account for the observed relic density. As said above, this mechanism is extremely predictive since, given the 'conformal' relation m 2 θ = λv 2 , the rates and hence the relic density depend only on the parameter m θ . Integrating numerically Eq. (15), it turns out that the observed relic density in Eq. (20) is obtained for where we took m h = 120 GeV for definiteness. This result corresponds to the right vertical branch of the thick curve in Fig. 1.  Figure 2: The decay thermalization rate, γ Dh /(n eq θ H) (black), compared to the various scattering thermalization rates, γ a annih /(n eq θ H) for a = W, Z, h, t, b (in red, blue, green, orange and purple respectively), as a function of z = m h /T and for m h = 120 GeV, m θ = 2.8 MeV. For freeze-in these rates remain always well below one, which is the thermalization value.
An important difference between freeze-in and freeze-out is that, while the latter is generally dominated by the off-shell annihilation, the former is naturally dominated by the on-shell annihilation, i.e. by the decay of the mediator involved in the annihilation. For a freeze-out, unless the mediator has a mass ∼ 2 m DM , the decay has a small effect because, when the off-shell annihilation goes out-of equilibrium, it has already decoupled (being more Boltzmann suppressed). For a freeze-in instead, the DM production occurs mostly at T of order the mediator mass, or slightly below, where the Boltzmann suppression is still mild. As a result the decay, which involves less (small) couplings, is naturally dominant, γ D h > γ annih . This is shown in Fig. 2 where the various reaction densities are plotted for m h = 120 GeV and m θ = 2.8 MeV. In the approximation where only the decay is included, the Boltzmann equation can be integrated analytically, giving [8] where g h = 1. This analytic result turns out to differ from the numerical result by less than 1%.
To understand Eq. (22) let us first discuss why the freeze-in is infrared dominated, with most of the θ production occurring at T m h (just before n eq h gets too much Boltzmann suppressed). The θ number density produced at temperatute T is essentially given by the number of decays per unit time occurring at this temperature, n eq h (m h /T )Γ(h → θθ), times the number of θ-particles produced per decay, N θ , times the Hubble time, 1/H. Therefore, taking into account thermal averaging properly, one has Y θ γ Dh /(Hs), which is nothing but the rate of thermalization γ Dh /(n eq θ H) shown in Fig. 2, times the relativistic value of Y eq θ . Thus one obtains Y θ ∝ M P m h Γ(h → θθ)/T 3 , which increases as T decreases. Hence one finds Y θ = c[N θ n eq h Γ(h → θθ)/(Hs)] T =m h , with c a numerical factor which turns out to be equal to (3π/2)/K 2 (1) 2.9. It is larger than one because, as Fig. 2 shows, the maximum value of the decay reaction density occurs at T ∼ m h /3.5 rather than at T ∼ m h . 4 With this in mind, and taking into account the dependence Γ(h → θθ) ∝ λ 2 /m h ∝ m 4 θ /m h , the DM relic density from the Higgs decay scales as m 5 θ /m 3 h , so that the value of m θ needed varies very little with the exact value of the relic density, and the larger the Higgs mass is, the larger m θ has to be. For example for m h = 140, 180, 300 GeV one gets the observed relic density for m θ = 3.0, 3.6, 4.9 MeV instead of 2.8 MeV in Eq. (21).
It is interesting to discuss what is the contribution of each scattering channel separately, even though they have a small effect with respect to the decay. As can be seen from Fig. 2, among the scattering channels, the W one gives a contribution slightly larger than the Z and h ones, and much larger than the b and t ones. Numerically, for m h = 120 GeV, the W, Z, h, b, t scatterings contribute to the total number of θ-particles produced in the proportions 1 : 0.4 : 0.4 : 4 × 10 −3 : 10 −1 respectively, while the W scattering alone gives a contribution 500 times smaller than the decay channel.
Similarly to the decay, at a temperature T m a , a scattering process aa → θθ gives a number density Y a θ γ a annih /(Hs), which is infrared dominated too. Again, this is nothing but the thermalization rate γ a annih /(n eq θ H), displayed in Fig. 2, up to the multiplicative constant Y eq θ . One gets Y a θ = c a [γ a annih /(Hs)] T =ma with c a a numerical factor which, unlike for the decay, does not take a unique value. It depends on m θ as well as on the s dependence of the scattering cross section (which fixes the position of the peak of the reaction densities in Fig. 2). For the W , Z and h channels c a 2.2, whereas for the fermion channels c t 1.7 and c b 6. One understands consequently why the W , Z and h channels, which give contributions Y θ ∝ M P /m W,Z,h , dominate over the b channel, which gives Y θ ∝ m 3 b M P /m 4 h . As for the top channel, it is suppressed too because for s 4m 2 t it is more threshold suppressed than the other channels, see Eqs. (9)- (12), and for higher s it gets quickly suppressed by the m 2 t /s 2 asymptotic behaviour of the top cross section (as compared to the 1/s asymptotic behaviour of the W , Z and h cross sections).
Let us now move to the region m θ < 3 MeV, where the Higgs portal alone would lead to a too small relic density. In this case the N annihilation process can do the job, as long as the reheating temperature is large enough to produce a thermal population of sterile neutrinos, which we assume. 5 The annihilation rate depends on two free parameters only, g and m θ , with the sterile neutrino mass m N determined by Eq. (4). From the freeze-in of 4 Actually if one takes the value of γ Dh /(Hs) at its maximum value, i.e. at T ∼ m h /3.5, one obtains a number density equal to 70% of the exact result. Therefore, the production of a particle through freeze-in of a decay is essentially given by the maximum value of the corresponding thermalization rate. Note also that if one takes instead N θ n eq h Γ(h → θθ)/(Hs) at its maximum value, which occurs for T ∼ m h /3.0, it turns out that one gets an even better approximation, as this overestimates the exact result by only 3%. 5 It is well-known that, if a sterile neutrino has Yukawa couplings inducing a light neutrino mass of the order of the atmospheric or solar neutrino mass scales, the decays/inverse decays of the sterile neutrino are in thermal equilibrium at T ∼ m N . There are various cosmological issues which can put constraints on thermal sterile neutrinos, see e.g. [14]. However they become relevant for sterile neutrino masses below the ∼ GeV scale, which are not pertinent in our case, see Figs. 3,4 below. the reaction N N → θθ one obtains as usual Y θ [γ annih /(Hs)] T =m N . One then finds that the DM relic density scales as g 6 m ν /m θ , as confirmed by the numerical integration of the Boltzmann equation.
Our numerical result in the (m θ − g) plane is shown in Fig. 1 for m ν = 0.05 eV. The diagonal branch of the black thick line corresponds to a freeze-in through the sterile neutrino portal, with the correct value of the DM relic density. The required value of the θ−N coupling is given approximatively by g 2 · 10 −3 (m θ /MeV) 1/6 (eV/m ν ) 1/6 . Above (below) the line the relic density is too large (small). Since the relic density is proportional to the sixth power of g, the predicted value of g depends very weakly on the cutoff Λ, which we fixed in all the numerical analysis, k = log(Λ 2 /m 2 N )/(8π 2 ) = 1.

Non-thermal production
In general, one may think of non-thermal, more model dependent, θ production mechanisms. A class of mechanisms is related to the phase transition at the scale of SSB f . At this epoch (part of) the energy stored in the scalar potential false vacuum may be transferred to the pNGBs. There could also be cosmic strings produced at the phase transition, which would decay into pNGBs, but their exact density depends on model details. Here we assume that these contributions to Ω θ are negligible. This is the case, in particular, if the universe underwent an inflationary phase at temperatures smaller than f (note that the reheating temperature can still be larger than m N = gf / √ 2, since the coupling g is much smaller than one in our scenario).
A potential source of non-thermal production of pNGB dark matter could be provided by the oscillations of the field θ around the minimum of its effective potential, in case the value of θ at high temperature is displaced from such minimum, a mechanism well-studied in the case of the axion. Let us summarize the axion case. The axion mass is very suppressed above Λ QCD , therefore at high temperatures the axion behaves as an exact NGB. In this case, at the end of the associated Peccei-Quinn phase transition at scale f P Q , the value of the axion field a can lie in any of the equivalent vacua described by 0 ≤ a/f P Q < 2π. Later, when m a (T ) becomes important with respect to the Hubble parameter H(T ), the field a begins to oscillate around zero, producing a coherent state of particles at rest with an associated Ω oscill a (for a review see Ref. [15]). One may naively think that this picture applies also to our scenario, with a negligible m θ (T ) at high temperature (where the EW symmetry is restored), and the oscillations beginning only at T ∼ TeV, when m θ is generated. If this were the case, taking the values of f and m θ in the range we considered for the freeze-in production, one would conclude that Ω oscill θ overcloses the universe, unless the initial value of θ is tuned to be much smaller than f , or inflation takes place at temperatures below the TeV scale.
However, the temperature dependence of m θ (T ) is very different from the one of the axion, whose mass is generated non-perturbatevely by an anomaly at Λ QCD . In our case, m θ is generated, instead, by an explicit breaking of U (1) X and therefore it does not vanish at high temperatures. On the contrary, we expect m θ (T ) to receive large thermal corrections. Even if the EW symmetry is restored for T TeV, there are contributions to m θ that are not proportional to the Higgs vev v, such as δm 2 θ ∼ g 2 y 2 Λ 2 /(16π 2 ) 2 (see the discussion at the end of section 4). These lead to m 2 θ (T ) constant · T 2 at high temperature. A quantitative estimate of a possible non-thermal production of θ would then require a non-trivial study of the thermal evolution of the field θ during and after the phase transition at scale f . Still, we just notice that m θ (T ) is typically larger than H(T ) already at T ∼ f . Thus, we argue that the field θ does not acquire a random initial value of order f , but rather it sits in the minimum of the potential already at high temperature, or in other words coherent oscillations are strongly suppressed. On the basis of this argument, we assume that Ω oscill θ is negligible, and that the thermal freeze-in dominates the DM relic density.

Summary
Our results on the θ relic density are summarized in Fig. 1, where we show the line Ω θ = Ω DM in the (m θ − g) plane. The behaviour of the line with the correct relic density can be easily understood. For m θ 0.15 keV, g must be larger than (or equal to) the value needed to thermalize θ. For progressively larger m θ , the value of g should correspond to less and less thermalization. For 0.15 keV < m θ < 3 MeV, one can see that the required coupling g is progressively smaller than the thermalization value, indicated by the thin diagonal line. When one approaches m θ 3 MeV, g becomes rapidly smaller as the Higgs portal begins to contribute to the θ number density. Once the Higgs portal produces by itself the correct amount of θ-particles, g must be small enough to make the sterile neutrino freeze-in negligible.
We conclude that the observed DM relic density can be obtained thermally for where the exact upper value depends on the Higgs mass, as discussed in section 3.3. Note that this upper bound holds as an absolute prediction as soon as the Higgs portal dominates the DM production. This is necessarily the case, in particular, if the reheating temperature lies below the sterile neutrino mass scale. More generally, the prediction m θ 3 MeV holds for any pNGB whose Higgs portal interaction λ gives the dominant contribution to its mass and relic density, independently of the associated global symmetry and of the source of explicit symmetry breaking which induces λ. In particular the SSB scale may be different from the seesaw scale. Of course the prediction for m θ would change if the number of DM particles was not negligible already at temperatures higher than the EW scale. In this case one has to produce less of them through freeze-in, which means that m θ has to be smaller. In other words, m θ 3 MeV constitutes an absolute upper bound for the freeze-in mechanism, and it holds as soon as the freeze-in production dominates over the initial population.

Approximate symmetries of the seesaw sector
In this section we consider the SM augmented with sterile neutrinos ν c i , with a global symmetry broken spontaneously at the seesaw scale, and we discuss in some detail the generation of a mass for the associated pNGB.
Let us consider the most general Yukawa interactions to be added to the SM in the presence of gauge singlet fermions, where l α are the lepton doublets and H is the Higgs doublet, whose neutral component acquires a vev v = 174 GeV. Here and in the following, the mass parameters m αj and M ij are intended in a generalized sense as dynamical scalar fields that may or may not acquire a non-zero vev. Thus, the above lagrangian has a global U (1) L symmetry with the following lepton number assignments: where we dropped flavour indices. When M acquires a vev, U (1) L is spontaneously broken and a massless NGB appears in the spectrum of the theory, the singlet Majoron [1]. One can write more explicitly where Φ is a complex scalar with L = 2, the vev ρ = f breaks spontaneously lepton number, and θ is the Majoron. There are various possible sources of U (1) L explicit breaking coming from other sectors of the lagrangian. There may be soft terms in the scalar potential involving Φ which break lepton number; their mass scale, whose size is arbitrary, determines the induced Majoron mass (see e.g. [7]). Another possible source of U (1) L explicit breaking are quantum gravity effects at the Planck scale M P , which can break in general all non-gauge symmetries. Assuming these effects are suppressed by powers of M P , they can be used to generate a Majoron mass at the keV scale [2,3], which has been extensively studied as DM candidate [3,4,5,6].
In this paper we will consider a different source of explicit breaking of global symmetries, provided by the set of the Yukawa couplings [9]. This allows to relate the size of the pNGB parameters to the fermion mass scales already present in the theory. We will focus on symmetries of the lagrangian in Eq. (24) other than U (1) L , with each lepton carrying in general a different charge. In this case the symmetries are respected only by some matrix elements m αj and M ij , and they are explicitly broken if some other matrix elements are nonzero. The mass and couplings of the pNGBs will be completely determined by the seesaw parameters and by the choice of the cutoff, since they can receive cut-off dependent quantum corrections.
The pNGB mass is, in general, quadratically sensitive to the cut-off. In order to understand the origin of quadratic divergences and the mechanism to remove them, consider first the explicit breaking of a U (1) symmetry in a theory with only one sterile neutrino ν c : where M a,b are real mass parameters. Here M a is generated by the spontaneous breaking of the U (1) symmetry associated with the NGB θ, while M b breaks this symmetry explicitly. It is instructive to compute the fermion loops generating the pNGB mass term, m 2 θ θ 2 /2, using the right-hand side of Eq. (27): if M b is zero, the two relevant one-loop diagrams cancel each other, as expected for an exact NGB. However, in the presence of the explicit breaking, a non-zero quadratically divergent contribution is left, which is given by The effective theory below the scale f may still contain a light scalar θ, since m θ can be parametrically small, but its NGB nature is obscured by the quadratic dependence on the details of the ultraviolet completion.
In the presence of more than one fermion family one can define several family-dependent U (1) symmetries. In general, they are explicitly broken by some mass matrix entries. However, certain U (1)'s will be broken only when several entries are non-zero at the same time: this is the key to reduce the degree of divergence of the radiative contribution to m θ , where θ's are the associated pNGBs [9]. This fact can be understood in the language of the effective potential V ef f , considering the full (active + sterile) neutrino mass matrix M: the term in V ef f quadratic in the mass matrix, ∼ Tr(MM † )Λ 2 , is invariant under certain U (1) symmetries, that is, it does not depend on the associated pNGBs. Thus, the potential of these pNGBs contains at most terms quartic in the mass matrix, ∼ Tr(MM † MM † ) log Λ 2 , which are only logarithmically sensitive to the cutoff (for an application of this idea to eV scale sterile neutrinos, with the pNGB playing the role of dark energy, see [16]). Of course non-abelian symmetries are also possible, with several NGBs and, potentially, qualitatively different phenomena, but this extension will not be needed for our purposes and will not be considered in this paper.
In order to identify the combination of matrix entries that induces a pNGB mass, we rewrite Eq. (24) with the replacement l α (H/v) → ν α , in the minimal case of two sterile neutrinos, which are sufficient for realistic light neutrino masses: Here a sum over n f = 3 active flavours is understood (α = e, µ, τ ). Let us identify the matrix entries whose phases can be removed. There are 2n f + 3 mass terms, which in general have a phase. One can absorb n f + 2 of these phases by redefining the n f + 2 neutrino fields present in Eq. (29). As a result, there remain n f + 1 complex matrix entries in Eq. (29) (either n f phases in m and one in M , or equivalently n f − 1 phases in m and two in M ). When one sets these n f + 1 entries to zero, there are no physical phases left. Now, suppose these n f + 1 complex entries vanish and the non-zero entries are generated by a set of scalar fields Φ a ≡ ρ a × exp(iθ a /f a )/ √ 2, acquiring a vev ρ a = f a . Then, by an appropriate redefinition of the lepton fields, one can remove all the phase fields θ a from the Yukawa lagrangian of Eq. (29). Thus, these fields will have only derivative couplings to the leptons, resulting from the redefined lepton kinetic terms. This means that the lagrangian has exact U (1) symmetries, broken spontaneously by Φ 's, and θ a are exact NGBs.
When one of the zero matrix entries is switched on, however, it is no longer possible to remove all the phases. This means that one U (1) is explicitly broken, with the associated pNGB θ acquiring non-derivative couplings. Thus, a non-zero m θ is generated by neutrino one-loop diagrams. To estimate m θ , consider the set of the (n f + 3) non-zero matrix entries. One can check that θ could be rotated away from Eq. (29) if any out of four of these entries were switched off. This means that a U (1) symmetry is recovered when any of them is put to zero, and therefore m θ must be proportional to the product of the four entries. As a consequence, quadratically divergent contributions to m θ turn out to be absent.
The θ mass can be controlled by the entries of M only, by those of m only, or by both. Let us specify three models with a U (1) X symmetry, which are representative of these three generic possibilities: (i) The symmetry is broken explicitly in the singlet neutrino sector. Two independent entries of the Majorana mass matrix M are allowed by U (1) X , while the third is forbidden. For example take X(ν c 1 ) = 0, X(ν c 2 ) = 1 and one scalar field Φ with X(Φ) = −2. Then M 11 is allowed and M 22 is generated by the vev of Φ, while M 12 = 0. The NGB coupling to neutrinos reads M 22 exp(iθ/f )ν c 2 ν c 2 /2. (The Dirac neutrino sector is not relevant here: the m entries may or may not respect U (1) X , in any case the leading symmetry breaking effects are controlled by M .) When M 12 is different from zero, neutrino loops give a non-zero contribution to m θ as long as they contain four appropriate mass insertions: where we included a loop suppression factor and µ is the renormalization scale, which can be taken of the order of the sterile neutrino masses. The cancellation of quadratic divergences as well as the structure of the non-vanishing contribution to m θ are dictated by our symmetry argument, and they can be both explicitly checked by computing the loops.
(ii) The symmetry is broken explicitly within the Dirac neutrino sector only. The entries on a given row of the Dirac mass matrix m are both allowed, while in the other rows only one entry is allowed. For example take X(ν e ) = 1, X(ν µ,τ ) = 3, X(ν c 1 ) = −1, X(ν c 2 ) = 1 and a scalar field φ with X(φ) = −2. Then m e1 is allowed, m e2 , m µ1 and m τ 1 are generated by the vev of φ, and m µ2 = m τ 2 = 0. (For simplicity we suppose that the scalar fields φ and Φ that can contribute to m and to M belong to different sets. The singlet neutrino sector is assumed to respect U (1) X : in the present example M 12 is allowed, thus giving an equal mass to the two sterile neutrinos.) This case reproduces closely the scenario discussed by Hill and Ross [9], which dealt with the quark Dirac mass matrix. Once the entries explicitly breaking U (1) X are switched on in m, one obtains (iii) The symmetry is broken explicitly only by the interplay of the Dirac and Majorana neutrino mass matrices. One possibility is that U (1) X allows for all the entries of M , e.g. taking the charges X(ν c 1 ) = −1, X(ν c 2 ) = 1 and X(Φ) = 2; this symmetry can allow only one between m α1 and m α2 , depending on the charge of ν α . Another possibility is that U (1) X allows for all the entries of m, but it forbids two entries of M , e.g. taking the charges X(ν α ) = X(ν c 1 ) = 0 and X(ν c 2 ) = −X(φ) = 0. As usual the mass of the NGB θ is generated when one zero entry is switched on. One finds Before moving to the other properties of the pNGB relevant for DM, some comments are in order to assess the soundness of the above estimate for the pNGB mass m θ . We have seen that, in general, in the presence of explicit breaking of the associated global symmetry, m θ can be sensitive quadratically to the cutoff of the theory. However, if the global symmetry is broken only by the contemporary presence of several couplings, then m θ shall be proportional to the product of all these couplings. In this case the Feynman diagrams contributing to the pNGB mass will have a lower degree of divergence and thus quadratic divergences vanish. This possibility was often employed in model-building in the past, and recently was extensively used in the context of little Higgs models, under the name of "collective breaking" (see Ref. [10] for reviews). The motivation is to stabilize the electroweak scale against large quantum corrections. In these models the Higgs is a NGB of a symmetry broken spontaneously at the scale f ew ∼ TeV, whose mass is generated at one-loop level by two couplings y 1,2 , whose contemporary presence breaks explicitly the symmetry. As a consequence, m 2 h ∼ y 2 1 y 2 2 f 2 ew log(Λ 2 /µ 2 )/(16π 2 ), that is, the sensitivity to the cut-off Λ is only logarithmic. However, in general, the quadratic divergence reappears at the two-loop level, with a correction to the Higgs mass δm 2 h ∼ y 2 1 y 2 2 Λ 2 /(16π 2 ) 2 , which indicates that Λ cannot be larger than ∼ 4πf ew , in order for the theory to remain natural.
It is instructive to compare this little Higgs scenario with our scenario, where the DM candidate θ is a pNGB associated with SSB at the scale f , of the order of the seesaw scale. In this case the pNGB mass generated at one-loop can be written schematically as m 2 θ ∼ g 2 y 2 v 2 log(Λ 2 /µ 2 )/(16π 2 ). The mass m θ is not proportional to f , but rather to the electroweak scale v. Besides, it can be much smaller than v because of the four powers of Yukawa couplings and the loop suppression. This is true even with a huge cutoff Λ, since the dependence on it is only logarithmic. However, at higher orders the quadratic divergence may reappear, e.g. through a two-loop diagram with a virtual Higgs exchanged across the neutrino loop, leading to m 2 θ ∼ g 2 y 2 Λ 2 /(16π 2 ) 2 . This result is not surprising, since the mechanism we adopted explains why the pNGB mass is related to the EW scale, but it does not address the stability of the EW scale against radiative corrections. In other words, as already remarked by Hill and Ross [9], the Higgs sensitivity to quadratic corrections also enters in m θ at higher order. In order for this two-loop correction to be negligible with respect to the one-loop estimate, one needs the cutoff of the Higgs boson loops, Λ H , to be smaller than ∼ 4πv. (By enlarging the global symmetry, it may be possible to remove also the two-loop quadratic divergence, see e.g. [16], but in general higher orders reintroduce the problem.) If one wants to stabilize a theory which includes a scale much larger than v (e.g. the scale f ), one must address the usual hierarchy problem, e.g. postulating supersymmetry broken at Λ susy ∼ TeV, or a strongly interacting sector that condenses at Λ c ∼ TeV generating dynamically the EW scale. In this paper we assume the stability of the EW scale is realized, and thus we can adopt our one-loop estimates for the pNGB mass and couplings, in order to study phenomenology.

A seesaw model leading to the Higgs portal
Let us consider in some detail a specific U (1) X symmetry of the neutrino sector, such that the pNGB θ acquires radiatively a coupling to the SM Higgs. Take the sterile neutrino charges X(ν c 1 ) = −1 and X(ν c 2 ) = 1, and a scalar field Φ with charge X(Φ) = 2, whose vev breaks the symmetry spontaneously. The interaction lagrangian involving the sterile neutrinos and the pNGB θ reads Here M 12 is a mass term allowed by U (1) X , while M 11 and M 22 are generated after SSB. The symmetry is broken explicitly by either m α1 or m α2 , depending on the U (1) X charge assigned to the lepton doublet l α . 6 Note that θ can be rotated away by rephasing the neutrino fields if m α1 = 0 or m α2 = 0. In this case the symmetry is restored and θ is a true NGB. The dependence on θ can be removed also when two independent entries of M are vanishing, therefore the symmetry breaking effects must vanish also in this limit. As we discussed in the previous section, this need of four different couplings to break the symmetry is the key to cancel quadratic divergences. 7 The interactions in Eq. (33) generate an effective potential for the pNGB θ. We assumed that M ij and m αi are real and positive (see discussion below) and we performed the one-loop computation of the effective potential, which can be written as 6 The choice of X-charges is made to realized a model of type (iii), in accordance with the classification of the previous section. This choice is not unique. An equivalent possibility is a U (1) X symmetry with charges X(ν c 1 ) = −1, X(ν c 2 ) = 1, X(l α ) = 0, broken spontaneously by a scalar field of charge X(φ) = 1. Then in Eq. (33) one should replace m α1 → m α1 e iθ/f and m α2 → m α2 e −iθ/f , while the entries M 11 and M 22 are independent from θ and represent the source of explicit symmetry breaking. It is easy to check that the two realizations are equivalent, since the lagrangian is the same up to a rephasing of the neutrino fields. 7 Note that the field θ can be rotated away from Eq. (33) by phase redefinitions, as long as the n f entries that break explicitly U (1) X are set to zero. This is slightly different from the general case discussed below Eq. (29), where n f + 1 entries must vanish, for the trivial reason that here the same phase field appears (with opposite sign) in two independent entries. The only physically relevant fact still holds: the product of four independent entries is needed to generate a mass for θ.
We find that the quadratically divergent contributions cancel explicitly, as expected, and the logarithmically divergent ones give with a renormalization scale µ ∼ M ij and up to tiny corrections of higher order in m αi /M kl . This type of models is particularly predictive, because the pNGB mass is generated by the same loops that generate λ, that is, m θ is obtained from Eq. (34) by replacing the Higgs with its vev, m 2 θ = λv 2 . By taking a common value m N = gf / √ 2 (m = yv) for each entry of the Majorana (Dirac) neutrino mass matrix, Eq. (35) reduces to Eq. (3), up to a factor 2n f accounting for the sum over flavour indexes. As we saw in section 3, the value of the coupling λ has a crucial role for the (partial) thermalization of θ and thus for the determination of its relic density.
A comment on the CP symmetry is in order. In general, the mass terms M ij and m αi may be complex, that is, they may carry phases that cannot be removed by a redefinition of the fields. These phases correspond to an explicit breaking of CP. In this case θ would have both scalar and pseudo-scalar couplings to the fermions, with characteristic phenomenological signatures [9]. Moreover, the effective potential would contain terms odd in θ, such as µθH † H (with µ v), that induce a vev for θ and a small mixing with the Higgs. This may endanger the stability of θ, since it would couple linearly to SM particles. Still, these couplings can be very small and in addition the θ decays may be kinematically forbidden if m θ is sufficiently small. In this paper we do not investigate this more complicated CPviolating possibility, and we rather assume that CP is a good symmetry of the seesaw sector. In this case M ij and m αi in Eq. (33) are all real and θ preserves its pseudo-scalar nature, since Eq. (33) is invariant under a CP transformation with θ → −θ. Usually the DM stability requires an additional (discrete) symmetry that forbids its decay into SM particles. In the present scenario the CP symmetry itself guarantees that θ couples linearly only to the heavy sterile neutrinos, and we will see in section 5 that this makes θ sufficiently long-lived.

Constraints on the pNGB lifetime
We first derive the couplings of θ to the SM fermions and gauge bosons, and compare its decay width into light SM particles with the lifetime of the universe. Subsequently, assuming that θ-particles account for the whole DM density, we discuss the more stringent astrophysical and cosmological constraints which exist on the decay of θ to neutrinos, electrons and photons.

θ couplings to the SM particles
The pNGB lifetime is determined by its couplings to light fermions. These come from the θ-N interaction, through the ν-N mixing induced by the Dirac neutrino masses. In this way θ can decay to light neutrinos (at tree level) and to charged fermions (at one-loop level). In turn, these couplings to SM fermions could induce, through triangle loop diagrams, couplings to SM gauge bosons. These decays should be sufficiently slow to make the DM lifetime longer than the age of the universe, τ 0 5·10 17 s. This provides interesting constraints on the seesaw parameters, as we now discuss.
Neutrinos. At tree level, θ couples only to light neutrinos as follows: Here X ij is the power of e iθ/f associated with the sterile neutrino mass matrix entry M ij . In the singlet Majoron model X ij = 1 for all i, j, therefore one finds µ ν = m ν ≡ −mM −1 m T , which is the usual seesaw formula. On the contrary, in our models based on a familydependent U (1) X symmetry, θ couples differently to each entry of M , as for example in the model of Eq. (33). This lagrangian leads to a total decay width (into both neutrinos and antineutrinos) where we neglected (m νi /m θ ) 2 corrections in the phase space, m νi being the light neutrino mass eigenvalues. In the Majoron case one obtains g 2 θνν = i m 2 νi /f 2 . This width is smaller than 1/τ 0 for g θνν 3 × 10 −19 (MeV/m θ ) 1/2 . Such tiny coupling is natural for θ, because the pNGB couplings are suppressed by the SSB scale f : one needs f 3 × 10 9 GeV(m ν /eV)(m θ /MeV) 1/2 , where we used the one-family approximation, g θνν m ν /f . To translate this bound in the (m θ − g) plane we insert the relation f = √ 2m N /g in Eq. (4), which gives The condition 1/Γ(θ → νν) > τ 0 excludes the region above the blue dashed line in Figs. 3,4. Charged fermions. The pNGB θ couples also to charged fermions, through EW oneloop diagrams. This effect arises because of the mixing between the sterile and the weaklyinteracting neutrinos, in particular it is also operative in the singlet Majoron model [1,17]. The coupling to quarks is induced by a one-loop θ − Z mixing diagram, with neutrinos in the loop. The coupling to charged leptons is generated by an analog Z-exchange diagram plus a triangle diagram with W -exchange. The resulting coupling can be written in a compact form in the one-family approximation, as follows: where the sign is + (−) for up quarks and charged leptons (down quarks), and G F is the Fermi coupling constant. This effective coupling can be suppressed only taking a small g = √ 2m N /f , since all the other parameters in g θff can be determined experimentally. The decay width is given by We remark that all models where the singlet Majoron is given a mass larger than one MeV and plays the role of DM candidate are constrained (or already excluded) by such decays into charged fermions. In our scenario, the requirement 1/Γ(θ → ff ) > τ 0 leads to an upper bound g 2(MeV/m θ ) 1/2 (eV/m ν ) (MeV/m f ), where we have assumed m θ 2m f for simplicity. This bound excludes the region above the red dashed line in Fig. 4. As a consequence, taking into account that g and m θ are related by Eq. (3), and that the Dirac mass m = yv cannot be larger than ∼ TeV, we conclude that m θ close to the EW scale would imply a lifetime shorter than τ 0 . In fact, even stronger constraints on the decay width into charged fermions come from astrophysical and cosmological observations (see section 5.2), which will lead to a stronger constraint, m θ < 1 GeV. In particular, in this class of models θ cannot be the ∼ 50 GeV DM candidate produced by the freeze-out of the λ interaction with the Higgs, discussed in section 3.2. One is left with the possibility of a sub-GeV DM candidate, because in this case the decays into charged fermion pairs are sufficiently slow (or forbidden kinematically), and the correct relic density can be generated by the freeze-in mechanism.
In the realistic three-family case, the coupling g θff is generically of the same order, but with a complicated dependence on flavour parameters. In particular one may argue that some cancellation can take place, to reduce the θ decay width. In addition, we have seen in section 4 that the pNGB mass generation depends crucially on the interplay between the flavour structures of the matrices m and M . Therefore, if one were to study the whole parameter space, one should know the explicit dependence of g θff on the mass matrix entries.
For illustration, we display the result for the model in Eq. (33), considering for simplicity only one lepton doublet. Writing the Dirac mass matrix as m ≡ (m 1 m 2 ) and keeping terms up to order (m i /M jk ) 2 , the effective coupling of θ to fermions through the mixing with the Z gauge boson is given by where T 3f is the third isospin component of the left-handed part of the fermion f , M 1,2 are the eigenvalues of the matrix M , which is diagonalized by a rotation of angle δ, defined by tan 2δ = 2M 12 /(M 11 − M 22 ), and we denoted c ≡ cos δ and s ≡ sin δ. Finally, the loop function K = K(M 1 , M 2 ) is given by Concerning the θ decay into charged leptons, one needs to add the contribution of the Wexchange diagram. The main phenomenological constraint comes from θ → e + e − . Therefore, in our simplified calculation we identify the active neutrino with the electron neutrino ν e . Then, the additional contribution to the effective coupling of θ to electrons is given by that carries a relative factor −2 with respect to the Z-exchange contribution. This is why the sum of the two contributions in Eqs. (41) and (45) gives g θeē /m e = g θuū /m u = −g θdd /m d , consistently with the one-family result in Eq. (39). Also, the functions F 1,2 are of the order m 2 i /M j ∼ m ν , so the couplings are of the same order as in Eq. (39). Still, cancellations between the various terms are possible leading to a suppression of the θ decay width for special flavour structures. This possibility may deserve a future investigation, since an appropriate family symmetry could in principle raise the pNGB lifetime, so that θ could become a viable DM candidate even for m θ > 1 GeV, a region where freeze-out could lead to the observed relic density. In this paper we do not invoke a family symmetry for the suppression of g θff to happen, but rather we rely on the one-family estimate given in Eq. (39).
Gauge bosons. The last effective coupling of the pNGB that may have important phenomenological consequences is the one to photons. Let us discuss first the limit of exact global symmetry, and then comment on the effect of explicit breaking. In general, the couplings of a NGB to gauge bosons are controlled by the gauge anomalies of the associated global symmetry. In the case θ is the Majoron, the global symmetry can be identified with B − L, because θ can be rotated away from the Yukawa interactions by rephasing all the SM fermions with a B − L transformation. Since B − L is anomaly-free with respect to the SM gauge symmetries, this rephasing does not generate anomalous couplings of the type θFF , rather the only leftover interaction is derivative, (θ/f )∂ µ J µ B−L , which is generated by the redefinition of the fermion kinetic terms. One concludes that the Majoron has no anomalous couplings to gauge bosons. 8 This discussion is easily generalized to the U (1) X symmetry that we consider in this paper. First, note that the sterile neutrino charges are irrelevant for gauge anomalies. In order to allow for fermion masses, we can take X(l α ) = −X(e c α ) = x α for α = e, µ, τ , and zero X-charge for the quarks. Then it is easy to check that U (1) X is anomaly-free with respect to electromagnetism and colour, and thus θ does not couple to photons and gluons. The anomaly with EW interactions is proportional to α x α and can also be taken equal to zero for simplicity. In summary, in the limit of exact symmetry our NGB does not couple to gauge bosons.
It is far more difficult to compute the pNGB couplings to gauge bosons in the presence of explicit symmetry breaking sources. One expects that such couplings will arise at some level, since there is no symmetry arguments to prevent them at all orders. However, the 8 Equivalently, one may remove θ from the Yukawa interactions by an L transformation; in this case one is left with the couplings ( , where c W,Y = 0 account for the EW anomalies of the lepton number symmetry. There is no physical difference between the two pictures. In particular one can check that (i) θ does not couple to weak gauge bosons at one-loop even in the L case, because the triangle diagrams generated by the derivative coupling to J µ L cancel c W,Y exactly; (ii) θ does not decay to two on-shell quarks even in the B − L case, because the quarks are vector-like under B − L. computation of the lowest order non-vanishing contribution is non-trivial. We have shown explicitly above that, in the present framework, θ couples at tree-level to neutrinos only, while, at one-loop, couplings to the Higgs boson as well as to charged fermions are induced. At two-loop there is a number of diagrams that connect θ to two gauge bosons. These diagrams do not necessarily add up to zero: on the one hand, they involve the neutrino mass parameters that break explicitly U (1) X , thus they are not expected to respect the anomaly argument above; on the other hand, however, these parameters may not be sufficient to induce an operator θFF already at two-loop order.
The computation of these two-loop diagrams is beyond the purpose of the present paper, because we will show that a decay θ → γγ induced at this order would be irrelevant for phenomenology anyway: once the constraints from θ → ff are imposed, the surviving DM parameter space is not further reduced by constraints on photons. To see this, we estimate the size of a two-loop contribution by taking the effective one-loop couplings in Eq. (39), and computing the usual fermion triangle diagrams with two final state photons. The effective lagrangian can be written as and the corresponding decay width is where G(x) = | arcsin 2 x|/x 2 [18]. Note that even in the limit of exact NGB, with m θ m f and thus G(x) ≈ 1, the sum does not give zero when the sign in Eq. (39) is taken into account. Therefore, we stress again that this is just a conservative estimate, and the computation of all the two-loop contributions may lead to a further cancellation. Neglecting the mass dependence in G, the requirement Γ(θ → γγ) < 1/τ 0 gives the order of magnitude constraint g 4×10 3 (MeV/m θ ) 3/2 (eV/m ν ). This bound is weaker than the one obtained below Eq. (40) from Γ(θ → ff ), as long as m θ 1 GeV, which is the region relevant in the present scenario.

Cosmological and astrophysical bounds on θ couplings
In practice, the DM lifetime has to be larger than the lifetime of the universe, because late DM decays affect several cosmological and astrophysical observations. To derive the corresponding bounds on our scenario we make the assumption that θ is all the DM in the universe. For the sake of simplicity, here too we take the approximation of one lepton family, barring large cancellations between the flavour parameters. A sub-GeV particle could decay into neutrinos, electrons and photons. We discuss these decay channels in turn.
Let us first consider the decay θ → νν. As it was noted in Ref. [5] such a decay could affect the expansion history of the universe, because it represents energy transfer from a nonrelativistic θ to relativistic neutrinos. Using SNIa and CMB data, one obtains the bound [19] Γ(θ → νν) < 4.5 × 10 −20 s −1 (or τ > 700 Gyr) .
Using Eq. (37), the corresponding bound on the coupling of θ to neutrinos is This upper bound can be improved in the range 30 MeV < m θ < 200 MeV [20], using searches for the diffuse neutrino supernova background by Super-Kamiokande. For masses m θ > 200 MeV the best limit come from atmospheric neutrino observations. Since the observed spectrum coincides with theoretical estimates, one can set an upper bound on g θνν [20,21]. We use the relation in Eq. (38) to translate these constraints in the (m θ − g) plane.
The blue region shown in Figs. 3,4 is excluded by the combination of the observational bounds discussed above.
Let us now consider the θ → e + e − decay. We adapt to our case the analysis performed in Ref. [21]. For m θ 20 MeV, the dominant physical process which constrains the parameter space is annihilation at rest, contributing to the 511 keV line. The limit is approximately given by MeV m θ Γ(θ → e − e + ) < 5 × 10 −27 s −1 .
Using Eq. (40), this limit leads to For 20 MeV m θ 1 GeV the dominant process is internal bremsstrahlung, i.e., photons radiated from the final electron or positron. One gets a slighly more stringent bound than the one shown in Eq. (50) [21]. The predicted value of the coupling g θee in our model, as follows from Eq. (39), is The corresponding observational bound is plotted in Figs. 3,4 as the red shaded region. Note that the value m θ 3 MeV, predicted by the Higgs portal, could lead to an excess of 511 keV γ rays from the galactic center of the Milky Way. As discussed in Ref. [22], a ∼ 1.5 MeV electron or positron, produced by a θ → e + e − decay, is enough non-relativistic to subsequently annihilate at rest in the galactic center, leading to a γ line at 511 keV. In particular, with such a mass one can easily obtain a γ flux excess of the order of the one observed by the INTEGRAL γ observatory [22,23]. However, unless the DM galactic profile is much more cuspy than the usually considered profiles, a DM decay gives a flux that is not sufficiently peaked around the galactic center [22,24] to be able to reproduce the morphology of the signal observed by INTEGRAL.
We discuss now the decay θ → γγ. For small masses, m θ keV, the photons energy is absorbed by the baryonic gas in the early universe and the processes of recombination and reionization are affected. One can use the analysis done in Ref. [25] (see also [26]) to bound the coupling g θγγ . For larger masses, m θ keV, photons are no longer absorbed and propagate freely. Their contribution to the isotropic diffuse photon background allows also to bound g θγγ [26,27].
For the region we are most interested in, stronger bounds can be obtained from the gamma-ray line emission limits from the Milky-Way central region. Indeed, in the range 40 keV< m θ <16 MeV [27], one obtains Γ(θ → γγ) < 10 −28 s −1 m θ MeV .
Using Eq. (47), this corresponds to Taking the estimate for g θγγ given in Eq. (47), all these constraints translate in an upper bound on g, which is shown in Figs. 3,4 by the brown dot-dashed curve. We should mention that there are astrophysical constraints on the couplings g θνν [28], g θee [29], and g θγγ [29,30], based on stellar energy loss due to θ emission, provided m θ is low enough to be produced in stellar interiors. These limits are valid without the need to assume that θ is DM. However, these astrophysical limits are in general weaker than the ones mentioned above.
We also would like to note that the region we consider for pNGB DM masses includes the keV scale, corresponding to warm dark matter (WDM). Since at the epoch of structure formation WDM has free-streaming lengths below the Mpc scale, having WDM at least as a non-negligible DM component can alleviate some of the disagreement between the standard cold DM scenario and a variety of galactic observations at small scales [31]. The ideal observation which could place a limit on WDM is Lyman-α forest, i.e. the Lymanα absorption produced by intervening neutral hydrogen in the spectra of light emitted by distant quasars. Using Ly-α observations together with other cosmological data sets, different groups have put lower limits on m W DM , assuming WDM has a thermal distribution and that it is the whole of DM. The bounds in the literature [32] differ by factors of a few, ranging from 0.5 keV up to 4 keV. In our case, the DM species is not in thermal equilibrium and the bounds can slightly change compared to the WDM thermal relic case. In addition, if WDM is only a part of the whole of DM the bounds are relaxed. We should mention that there are some potential problems in the obtention of these bounds, and they should be regarded as controversial. On the one hand there could be large systematic errors, and on the other hand the Ly-α analysis has to be performed at scales which are already in the non-linear regime, where calculations are less reliable. If one disposes of the Ly-α data these bounds on WDM disappear altogether [33].

Conclusions
We proposed a new pseudo-scalar gauge-singlet DM candidate θ with mass in the keV -MeV range. Its couplings to the SM particles are feeble, because they are mediated by new physics at a large scale f , which we identify with the seesaw scale. The θ relic density can be produced by the scattering with the heavy particles (the sterile neutrinos), at temperatures of the order of f , or alternatively it can be generated at the EW scale through the Higgs portal, by the tiny θ-H coupling λ ∼ 10 −10 , which is induced by the seesaw interactions. Today, θ decays into light neutrinos and, if it is heavier than one MeV, into e + e − , with rates that can saturate the present upper bounds.
We argued that such a candidate is theoretically well-motivated. The heavy new physics sector is generically associated with several global symmetries. Some of the corresponding pNGBs may remain light, if the explicit symmetry breaking effects are sufficiently small. For concreteness, we demonstrated that a U (1) X family symmetry of the sterile neutrino sector can be broken collectively by a set of neutrino Yukawa couplings, so that the pNGB mass m θ is proportional to the EW scale and to the (small) product of the Yukawa couplings. In such scenario m θ is not quadratically sensitive to the cutoff (at leading order), therefore the presence of a light DM scalar below the EW scale is justified.
In order to calculate the θ relic density, we computed the rate of the θ interactions with a sterile neutrino N and, through the Higgs portal, with the SM particles, and we studied numerically the Boltzmann equation for the θ number density. In this framework there are only two independent parameters, the mass of the pNGB, m θ , and the coupling of the pNGB to the sterile neutrino, g, a feature which makes our scenario especially predictive. We find that the Higgs portal produces the desired DM relic density, through the freeze-in mechanism, for a unique value of the DM mass, m θ 3 MeV. This prediction relies only on the relation between the DM mass and its coupling to the Higgs, m 2 θ = λv 2 , which is well-justified in the case of our pNGB. As long as the reheating temperature is sufficiently high, the sterile neutrino portal can produce the required relic density by freeze-in for smaller values of m θ , from 3 MeV down to 0.15 keV. These results are summarized in Fig. 1.
The constraints from Γ(θ → νν), however, exclude the region m θ 1(100) keV, for a light neutrino mass m ν = 0.05(1) eV. In addition, the constraints from Γ(θ → e + e − ) put an upper bound on the coupling g in the region m θ 1 MeV: for m ν = 0.05(1) eV and g ∼ 10 −3 (10 −4 ), the expected electron-positron flux is close to the present sensitivity. These results are summarized in Figs. 3,4. In turn, the seesaw scale can be constrained, by requiring our pNGB to be a viable DM candidate in agreement with all the bounds above. In the case of freeze-in through the sterile neutrino portal, we find 10 5 GeV m N 10 10 GeV, and a SSB scale f ∼ 10 3 m N . In the case of the Higgs portal, one has instead 10 10 GeV m N 10 14 GeV, and correspondingly 10 13 GeV f M P .
Finally, we note that our analysis of the parameter space was mostly performed in the one-family approximation, that is, neglecting possible hierarchies among the sterile neutrino mass parameters and the neutrino Dirac Yukawa couplings. Some of our main results are independent from this approximation, such as the allowed range for the DM mass. On the contrary, the bounds on g and m N are obviously sensitive to the flavour structure. A rough idea of these dependence can be grasped by comparing the case m ν = 0.05 eV (Fig. 3) with the case m ν = 1 eV (Fig. 4). A more detailed exploration of the flavour parameter space of this scenario will be desirable, in particular, if one wants to compare with neutrino flavour models.