Inverse Seesaw, dark matter and the Hubble tension

We consider the inverse Seesaw scenario for neutrino masses with the approximate Lepton number symmetry broken dynamically by a scalar with Lepton number two. We show that the Majoron associated to the spontaneous symmetry breaking can alleviate the Hubble tension through its contribution to $\Delta N_\text{eff}$ and late decays to neutrinos. Among the additional fermionic states required for realizing the inverse Seesaw mechanism, sterile neutrinos at the keV-MeV scale can account for all the dark matter component of the Universe if produced via freeze-in from the decays of heavier degrees of freedom.


Introduction
A significant and intriguing tension between local, late-time determinations of the Hubble rate and its preferred value when measured from early Universe probes persists. Analyses from type-Ia supernovae and strong lensing consistently favor values of H 0 significantly larger than those determined from the cosmic microwave background (CMB) and baryon acoustic oscillations data. A tension at the level of 4−6 σ [1,2], depending on the specific assumptions, exists between the Planck value from the CMB spectrum [3] and the one obtained by the SH 0 ES collaboration [4] from supernovae measurements. Among the many different solutions proposed [5], those that also address another open problem are particularly appealing. For instance, the authors of Refs. [6,7] have proposed that a light Majoron contributing to ∆N eff and decaying to neutrinos after Big Bang Nucleosynthesis (BBN) may alleviate the discrepant determinations of H 0 . This scenario would thus link the solution to the Hubble tension to the origin of neutrino masses and mixings.
Indeed, neutrino masses and mixings as required by the observation of the neutrino oscillation phenomenon, can be naturally accounted for through the Weinberg operator [8]. This operator violates Lepton number (L) by two units. Thus, if this breaking is dynamical, the Majoron, the Goldstone boson associated to the spontaneous symmetry breaking of Lepton number, would be intimately linked to the origin of neutrino masses and mixings.
Among the many possible realizations of the Weinberg operator, a particularly interesting one is the inverse Seesaw scenario [9][10][11], which assumes an approximate Lepton number symmetry. In this way, the smallness of neutrino masses is naturally explained by this symmetry argument and a very high scale for the new physics is not required, leading to potentially more interesting phenomenological consequences.
In this work we investigate a dynamical origin for the small Lepton number breaking of the inverse Seesaw scenario. Several constructions based on dynamical Lepton number breaking have been explored in the past [12][13][14][15][16][17][18][19][20]. In this work, we consider a Seesaw-like mechanism in the scalar sector so that a small vacuum expectation value is naturally induced for the scalar with Lepton number two [14]. Thus, neutrino masses will be proportional to this small parameter. The dynamical symmetry breaking will also imply the presence of a Majoron with the potential of alleviating the Hubble tension.
Furthermore, inverse Seesaw realizations may also lead to sterile neutrinos at the few keV scale which are good dark matter (DM) candidates [21][22][23]. While production via mixing [24] is ruled out 1 by the stringent constraints from searches of X-ray lines [28], the correct DM relic abundance could be achieved from the decay of heavier states instead [21][22][23][29][30][31][32][33]. After investigating this possibility, we find that the couplings of the scalars to neutrinos can lead to both the correct DM relic abundance through freeze-in via decays and the necessary Majoron population so as to alleviate the Hubble tension. This paper is organized as follows, in Section 2 we introduce the particle content and Lagrangian of the inverse Seesaw model with dynamical breaking of L considered. In Section 3 we discuss the dark matter production mechanism as well as its phenomenology and constraints in the parameter space. In Section 4 we analyze the conditions under which also the Hubble tension can be alleviated as proposed in Refs. [6,7]. Finally, in Section 5 we discuss and summarize the allowed regions of the parameter space while in Section 6 we conclude.

The model
The simplest extension of the Standard Model (SM) particle content to account for the neutrino masses and mixing evidence is the addition of fermion singlets, i.e. right-handed neutrinos. Furthermore, if a large Majorana mass is also included for the right-handed neutrinos, the smallness of neutrino masses is naturally explained through the hierarchy between this mass and the electroweak scale via the celebrated canonical type-I Seesaw [34][35][36][37]. Conversely, a Majorana mass significantly above the electroweak scale destabilizes the Higgs mass, worsening the electroweak hierarchy problem [38,39] and greatly hinders the testability of the mechanism.
It is therefore appealing to consider low-scale alternatives to the canonical Seesaw mechanism. This option was investigated in Refs. [40][41][42][43] showing that the smallness of neutrino masses can also be naturally explained by an approximate Lepton number symmetry. Two types of fermion singlets may be included according to their L assignment. The first option is Dirac pairs with L = 1 for which, in the limit of exact L, only the right-handed component may have a Yukawa coupling to the active SM neutrinos. The second option is Majorana sterile neutrinos with L = 0 which, for exact L symmetry, do not couple to any other fermion. At this level, three neutrinos remain massless. When the L symmetry is slightly broken, small neutrino masses can be induced, the Dirac neutrinos may split into pseudo-Dirac pairs, and additional suppressed couplings are allowed.
Following this principle, we extend the SM particle content with Dirac pairs with their corresponding right and left-handed components N R and N L . At least two of these pairs are required to reproduce the correct neutrino masses and mixings (in this case with the lightest neutrino massless). We also consider one Majorana sterile neutrino n L . The L violation that will induce the standard neutrino masses will be dynamical and originated through two scalars φ 1 and φ 2 . All the new states are singlets of the SM gauge group and their Lepton number charge assignment L is given in Table 1. According to this assignment the Lagrangian of the model can be parametrized as: where L L are the SM lepton doublets and H the Higgs doublet.

Scalar potential
In the inverse Seesaw mechanism, neutrino masses are protected by an approximate L symmetry. The smallness of the L-violating terms thus naturally explains the lightness of neutrino masses. Since we want to explore a dynamical breaking of this symmetry, we also consider a Seesaw-like mechanism in the scalar sector to avoid hierarchy problems and account for the smallness of the L-violating vacuum expectation value (vev) in a technically natural way. To this end, we mimic the type-II Seesaw [44][45][46][47] and assume that the vev of φ 2 will be induced by that of φ 1 as in Ref. [14]. In particular, the scalar potential is given by The physical pseudoscalars are given by where J is the Goldstone boson associated to the breaking of L, that is, the Majoron, and therefore massless from the scalar potential. Since L is expected to be broken from gravity effects [48], we will assume that a Majoron mass of the order of the eV scale is induced by them. The mixing angle β between a 1 − a 2 is:

Neutrino masses
When all the scalars develop their vevs, v H , v 1 and v 2 respectively, the neutrino mass matrix takes the inverse Seesaw form: arranging the states as ν L n L N L N c R . From now on we will work in the basis where M is diagonal. The approximate expressions for the flavour states in terms of the mass eigenstates are: where ν i , ν 4 , N + and N − are the mass eigenstates with masses with i = 1, 2, 3, θ ≡ m D M −1 characterizing the mixing between the active flavours ν L and the heavy states N ± , and U the unitary matrix diagonalising the light neutrino mass matrix after the block diagonalisation. We have assumed that M m D µ. In particular, we will assume that M is somewhat above the electroweak scale and that it controls the scale of the pseudo-Dirac pairs N ± . The splitting of the pseudo-Dirac pairs is only through the Majorana masses µ LL and µ RR . We will also assume that µ ss is at the keV scale and is the main contribution to m ν 4 , the dark matter candidate mass. For a summary of the approximate ranges of all the model parameters to correctly reproduce neutrino masses and mixings, the DM relic abundace and to improve on the Hubble tension see Table 2 where we sumarize our findings of the following sections. According to these values, all the µ parameters have been treated as a perturbation in the expressions above and the results are to leading order in perturbation theory. Furthermore, we have also approximated the results to leading order in θ to simplify the expressions. Notice that U , the rotation diagonalising the light neutrino mass matrix, m ν i , corresponds to the PMNS mixing matrix at leading order.

Dark matter
The Majorana fermion singlet n L , with its L = 0 charge assignment, can only mix with the other neutrinos via L-violating, and therefore suppressed, parameters. Furthermore, its allowed interactions with N L and N R are via φ 1 through the Y Ls and Y Rs parameters respectively. These two parameters, together with η, are all protected by an additional symmetry. Indeed, setting the three of them to zero a new U (1) transformation for φ 1 , independent from L, becomes a symmetry of the Lagrangian. We will therefore consistently assume small values for these three parameters. As previously discussed, a small value of η guarantees that the induced vev v 2 will be suppressed and thus naturally explain the smallness of neutrino masses. Small values for Y Ls and Y Rs in turn imply that interactions and decays of n L are very suppressed, making it an ideal DM candidate via freeze-in production. In this way, the same symmetry behind the smallness of neutrino masses also guarantees DM stability in a natural way. In the following we will discuss the production mechanism as well as the main constraints from it and other observations on the parameter space of the model.

Dark matter production
The DM candidate in our model is the mass eigenstate ν 4 which is approximately aligned with the fermion singlet n L with only suppressed mixings with the rest of the interaction eigenstates given the approximate L symmetry, as shown in Eq. (2.11). While the active flavour eigenstates ν L do contain an admixture of the DM candidate ν 4 as given by Eq. (2.11), processes that produce ν L such as decays of the Z and W or of the heavy neutrinos N ± via their Yukawa interactions with the Higgs, will not contribute to the production of ν 4 beyond the standard Dodelson-Widrow mechanism [24]. Indeed, the active flavour eigenstates ν L are already in thermal equilibrium in the early Universe and additional contributions such as these will not modify their abundance. In other words, the thermal masses of the active neutrinos ν L are very relevant in the early Universe and dominate over the keV-scale mass of ν 4 , suppressing the mixing [49]. That is, the interaction eigenstates are approximately the effective mass eigenstates [50]. Therefore, in this regime, it is more convenient to work in a mixed basis with N ± together withν L andn L : the "incomplete" flavour states ν L and n L at energies below m N ± . In this intermediate basis the original interaction eigenstates read: Sincen L does not share the relevant contributions to the thermal masses withν L , it is through processes in whichn L is produced where contributions to the final DM abundance of ν 4 beyond the Dodelson-Widrow mechanism can be achieved. The main interactions ofn L are with the heavy pseudo-Dirac pairs N ± and the new scalar particles via the couplings Y Ls and Y Rs . Thus, given the smallness of Y Ls and Y Rs , the main production channel for DM is through freeze-in [51] decays of the heavy neutrinos to a DM state and S, with S = ϕ 1(2) , A, J any of the physical scalar degrees of freedom. We find that the total dark matter production rate is where Θ is the Heaviside step function and c 12 ≡ cos α 12 and c β ≡ cos β. We have made the approximation m N + ∼ m N − ≡ m N . If the heavy neutrinos thermalize with the SM plasma, which happens for all the values of θ we will consider, the relic density is given by where M P = 1.22 · 10 19 GeV is the Planck mass, s 0 and ρ 0 c are the present entropy density and critical energy density respectively, h is the present Hubble constant expressed in units of 100 km s −1 Mpc −1 and g (m N ) is the number of radiation degrees of freedom during the N i decays, which we approximate as g (m N ) = 106.75 for m N 100 GeV. This expression can be simplified to study the analytical scaling of the relic density in two opposing limits of the following ratio r ≡ µ Rs µ Ls .
Nevertheless the full expression of the relic density has been taken into account in the numerical results. Assuming small r 1 and that the decay width is dominated by the J and ϕ 1 final states, the relic abundance scales as Ω ν 4 h 2 (r 1) 0.13 m ν 4 10 keV The y-axis is the mixing angle between the dark matter mass eigenstate and the SM neutrino flavour eigenstates and the x-axis the dark matter mass. Black dashed lines represents the correct relic density for a fixed ratio µ Rs /µ Ls . The three blue shaded regions dubbed "X-ray" represent respectively, from left to right, constraints from XMM-Newton [52], NuSTAR [28] and INTEGRAL [53] on spectral photon lines generated by decaying dark matter. The grey region corresponds to Lyman-α constraints from the matter power spectrum on light free-streaming dark matter particles, estimated from [54]. The red region represents the parameter space for which the dark matter lifetime is shorter than the age of the Universe, estimated using Eq. (3.9). In this figure we fixed m N = 150 GeV, where we have neglected the scalar masses with respect to m N and approximated c β ∼ c 12 ∼ 1 as well as written µ Ls in terms of U α4 ∼ θµ Ls /m ν 4 , the mixing between the DM candidate ν 4 and the active neutrinos ν Lα with α = e, µ, τ (see Eq. (2.11)) for which strong constraints exist from X-ray searches.
In the opposite limit r 1 and assuming dominant decays to ϕ 2 and A, again neglecting their masses with respect to m N , the relic density scales approximately as In Fig. 1 we show the allowed regions of the parameter space allowed by the X-ray searches and Lyman-α forest constraints on neutrino DM as well as the values of µ Rs /µ Ls for which the correct relic abundance would be obtained for different values of θ. The parameter θ represents the mixing between the active neutrinos and the heaviest mass eigenstates and can be bounded to be θ ≤ 10 −2 from flavour and electroweak precision tests of the unitarity of the PMNS matrix [55]. As can be seen, for values of θ close to the current upper bound in Fig. 1, a sizable hierarchy of about four orders of magnitude between µ Rs and µ Ls would be needed in order to obtain the correct relic abundance through µ Rs . Indeed, the stringent X-ray constraints require sufficiently suppressed active neutrino-DM mixing, which is mainly 3 dominated by µ Ls . Conversely, this hierarchy is avoided for smaller values of θ. This choice may be considered more natural, since there is no reason for a significant hierarchy among these parameters from the charge assignments of the fields. However, a lower bound on θ can be extrated from the requirement of perturbative unitarity for the Yukawa coupling Y LL by the relation implying that for v 2 1 GeV, θ has to be larger than O(10 −5 ) to ensure perturbativity for Y LL . Moreover, small values of θ would reduce the testability of this region of the parameter space, at least through unitarity constraints of the PMNS matrix or direct searches of the heavy neutrinos at colliders. In this regime, the dominant phenomenology of the model would rather correspond to DM searches via X-rays as well as through cosmology from its impact on the H 0 tension and contributions to ∆N eff .

Dark matter decay ν 4 → ν i + γ
Sterile-neutrino like dark matter can decay into a neutrino and a photon producing a monochromatic spectral line. The dark matter mixing with active neutrinos, as given by Eq. (2.11), is constrained by the International Gamma-Ray Astrophysics Laboratory (INTEGRAL) [53] by looking for DM decaying in the Milky Way halo, as well as from NuSTAR [28] and XMM-Newton [52]. These constraints correspond to the blue regions in Fig. 1.

Dark Matter lifetime
Notice that, apart from the usual decay channels to three light neutrinos or a neutrino and an X-ray photon, DM may also decay to a Majoron and a light neutrino. Thus, the associated lifetime of the DM needs to be larger than the age of the Universe. The decay rate is given by (3.9) The stability condition τ ν 4 ≡ Γ (ν 4 → J + ν i ) −1 > τ universe excludes the parameter space corresponding to the red region depicted in Fig. 1. Nonetheless, notice that the constraints on the mixing from X-ray searches are always stronger than those from the DM lifetime. In our scenario the DM density is generated from the decay of non-relativistic heavy neutrinos thermalized with the SM plasma. The effect of the resulting non-thermal phase space distribution on the matter power spectrum has been studied in various works [23,29,[63][64][65][66].

Constraints from the power spectrum and Lyman-α
The Lyman-α constraints on our DM candidate can be expressed, following the procedure of [54], as where g s (T ) is the temperature-dependent effective number of entropy degrees of freedom. This constraint is represented by the grey band in Fig. 1 for the reference value m Ly−α WDM = 3 keV but can be straightforwardly translated to a different value using Eq. (3.11).

The Hubble tension
The solution proposed in Refs. [6,7] to alleviate the present Hubble tension contains two key ingredients. The first is a contribution to ∆N BBN eff ∼ 0.4 that should already be present during Big Bang Nucleosynthesis. The second ingredient is an interaction rate between the Majoron and neutrinos that would exceed the Hubble rate between the BBN and CMB epochs. Thus, Majorons will thermalize with neutrinos for temperatures close to the Majoron mass T ∼ m J by decay and inverse decay processesν i ν i ↔ J. After becoming non-relativistic, Majorons would subsequently decay into neutrinos, resulting in a slight increase of ∆N eff . In addition to this extra late radiation component, Majoron-neutrino interactions cause a damping of the neutrino free streaming by suppressing their anisotropic stress and therefore affect the determination of the Hubble constant from the CMB.

Majoron contribution to ∆N eff
Since the scalar ϕ 1 mixes with the SM Higgs and also couples to the Majoron, the latter will be produced both from interactions with SM fermionic states ψ mediated by virtual ϕ 1 as well as from ϕ 1 decays when it is present in the bath through the following couplings, respectively: Such couplings allow to maintain Majorons thermalized with the SM plasma until the freezeout temperature T FO below which the Majoron population decouples from the thermal bath and behaves as background radiation, potentially leading to a contribution to ∆N eff that can alleviate the Hubble tension. Both scatterings and (inverse) decays can allow the light scalars to thermalize and are investigated in the following.
Thermalization via scattering In order to estimate the freeze-out temperature in this case, one can compare the expansion rate of the Universe to the typical momentum-exchange rate induced by the coupling from Eq.  Figure 2: Parameter space allowing to alleviate the Hubble tension (white). The blue region corresponds to values than cannot alleviate substantially the Hubble tension, while the green area represents BBN constraints on ∆N eff estimated in Ref. [7]. The orange region is excluded by constraints from ATLAS [67] on Higgs invisible decay as detailed in Sec. 4.3. The dashed purple line is determined from the (inverse) decay thermalization critera of Eq. (4.7) and the dash-dotted purple line corresponds to m ϕ 1 < T scat FO with T scat FO obtained from Eq. (4.4). The purple region represents the parameter space beyond the range of validity of the analysis described in Sec. 4.1 to determine the contribution to ∆N eff and for which a more elaborated estimate should be performed.
where we follow the notations and conventions from the appendices of Ref. [18]. In the relativistic limit, this quantity can be estimated as As argued in Ref. [68], such an estimate of the freeze-out temperature is rather reliable given the large temperature dependence of Eq. Thermalization via (inverse) decay A population of ϕ 1 could also be produced by inverse decay of SM fermions ψ. Given the fact that Γ ϕ 1 →JJ > Γ ϕ 1 →ψψ , if the ϕ 1 production rate induced by inverse decay is sizable enough to ensure ϕ 1 thermalization with the SM plasma, the Majorons J should thermalize as well. Once the scalars ϕ 1 become non-relativistic and their abundance is exponentially suppressed, that we estimate to be around z ≡ m ϕ 1 /T 5, thermal equilibrium with the SM plasma is lost and the population of J also freezes-out. As the coupling between the scalar ϕ 1 and ψ is proportional to m ψ , this process would be relevant for our parameter range mostly when the ϕ 1 decay channel to muons is open, i.e. for m ϕ 1 > 2m µ . The decay width of ϕ 1 to a pair of SM fermions is given by where c ψ is a colour factor. As detailed in [18], the Boltzmann equation relevant for ϕ 1 production by (inverse) decays can be expressed in term of the yield Y ϕ 1 ≡ n ϕ 1 /s is the thermal-equilibrium expected value of the yield and K 1,2 (z) are modified Bessel functions of the second kind. The values of the coupling ϕ 1 − ψ allowing to reach thermal equilibrium at z = 5 are given by the condition This condition has been checked numerically and yields a rather conservative constraint on the parameter λ 1H . For larger couplings than the benchmark point of Eq. (4.7), a J population would thermalize with SM fermions and freeze-out at T FO m ϕ 1 /5.
Contribution to ∆N eff In the parameter space for which the condition of Eq. (4.7) is satisfied, the (inverse) decay processes are more efficient than scatterings to maintain thermal equilibrium. The resulting contribution to ∆N eff is [18] ∆N eff 0.29 with g s (m µ ) 17.6. The ∆N eff range found in [6,7] that alleviates the Hubble tension is between 0.2 and 0.5, with a preferred value of 0.37. In Fig. 2 we depict the values of λ 1H and m ϕ 1 that would lead to such a contribution bounded by the blue line corresponding to ∆N eff = 0.2 and the green line corresponding to ∆N eff = 0.5. The white region represents the parameter space that allows to alleviate the Hubble tension. We emphasize that the boundaries of this white region of parameter space might be subject to small corrections given the order of magnitude estimate presented in this section. Nevertheless, a region with 0.2 < ∆N eff < 0.5 will be present in that area, since we verified that at least one of the two processes analyzed would allow Majorons to thermalize down to the temperature required. In particular, in Fig. 2, in the region between 2m µ < m ϕ 1 < 700 MeV the (inverse) decays of ϕ 1 allow them to thermalize with the SM bath and Majorons as long as λ 1H is above the dashed purple line. Hence, the vertical boundaries from the green and blue areas correspond to when we estimate that ϕ 1 becomes Boltzmann-suppressed and decouples so that also the Majorons freeze-out with 0.2 < ∆N eff < 0.5. Conversely, in the white triangle below the dashed purple line, scatterings with SM fermions are able to keep the Majorons in equilibrium instead with a final contribution to ∆N eff in the same range. In some areas of the parameter space both processes may be relevant simultaneously but such a detailed analysis is beyond the scope of this work and should not lead to sizable deviations from Fig. 2.

Majoron interactions with neutrinos
The authors of Ref. [7] define an effective width normalized such that for Γ eff 1 Majorons do thermalize with the active neutrinos as required to alleviate the Hubble tension: where λ ν is the dimensionless Majoron-neutrino coupling: In our setup this parameter corresponds to where we have neglected the contribution from µ Ls to m ν i . In this approximation The best fit for Γ eff found in Ref. [7] depends slightly on the number of active neutrinos interacting with the Majoron. Indeed, notice from Eq. (4.11) that λ ν is proportional to the neutrino mass. Thus, if the lightest neutrino is very light or massless, for instance if only two N R -N L pairs are considered, its coupling to the Majoron would be negligible. In particular the best fit changes from Γ eff = 67.6 to Γ eff = 59.9 when 2 or 3 neutrinos are considered to interact with the Majoron respectively. The dependence of Γ eff with ∆N BBN eff was found to be stronger. Indeed, the best fit Γ eff = 67.6 corresponding to ∆N BBN eff = 0.37 jumped to Γ eff = 678 for ∆N BBN eff = 0.48, although this larger contribution to N BBN eff significantly worsened the fit. In all cases the best fit for the Majoron mass was m J ∼ 0.3 eV. The preferred values of ∆N BBN eff = 0.37 and Γ eff ∼ 60 can easily be achieved as shown in Fig. 2 and Eq. (4.12).

Constraints from Higgs invisible decay
A coupling between the Higgs and the light scalars J, ϕ 1 is generated via mixing from the kinetic terms of φ 1 and the λ 1H term of the scalar potential as Via these couplings, the Higgs can decay invisibly to a Majoron or ϕ 1 pair with a decay rate where we replaced the mixing angle α 1H by its analytical approximation in the limit of small mixing. The invisible branching ratio of the Higgs is constrained to be B(h → inv) < 0.11 (0.19) from ATLAS [67] (CMS [69]) which translates into

Summary of the available parameter space
Taking into account all the constraints discussed in the previous sections, we sketch in Tab. 2 the ranges for the parameters of the model in which all conditions may be satisfied so that the correct neutrino masses and mixings and dark matter relic density can be recovered together with an improvement of the Hubble tension. Neutrino masses are controlled by the product θ 2 µ LL . The parameter θ represents the mixing between the active neutrinos and the heavy pseudo-Dirac pairs and is bounded to be θ ≤ 10 −2 from tests of the PMNS unitarity via precision electroweak and flavour observables [55]. Conversely, if θ is too small, the heavy pseudo-Dirac pairs, which populate the DM abundance via their decays, would not thermalize. This fixes the range for this parameter between roughly 10 −4 and 10 −2 . We have shown in Fig. 1 that the correct relic abundance can be obtained for θ = 10 −2 , 10 −3 , but it can also be recovered for smaller values of θ. The parameter µ LL should then correspond to m ν i /θ 2 , with values in the keV to MeV range. The value of µ LL in turn comes from the breaking of L by two units of the vev of φ 2 , v 2 . Thus, assuming order one Yukawas, v 2 and µ LL will have a similar range as reflected in Tab. 2. Finally, v 2 is induced by the vev of φ 1 , v 1 , through the η cubic coupling so that v 2 ηv 2 1 /m 2 ϕ 2 . The most natural choice for these parameters is to assume that v 1 and m ϕ 2 are close to the electroweak scale, so as to avoid hierarchy problems. Thus, the suppression in v 2 stems from the relative smallness of η, since this parameter is protected by the additional U (1) symmetry which is gained when this parameter together with µ Ls and µ Rs are set to zero.
Regarding the generation and properties of DM, the most stringent constraint is on the parameter µ Ls . Indeed, the mixing of DM with the active neutrinos U α4 ∼ θµ Ls /m ν 4 induces its decay to X-rays, for which stringent limits exist as shown in Fig. 1. In particular, for m ν 4 ∼ 10 keV, µ Ls is constrained to be between the eV and keV scales, depending on the value of θ. On the other hand, the DM relic abundance is controlled by µ Ls and µ Rs and a value around 1 keV is required. Thus, as shown in Fig. 1, either µ Rs is significantly larger than µ Ls , or θ ≤ 10 −4 . In addition, the DM abundance is induced by the decay of the heavy neutrinos to DM and some scalar degree of freedom. Therefore, the mass of the heavy neutrinos m N should not be much higher than the TeV scale to avoid further suppressing the DM relic abundance. For large θ, these neutrinos could be searched for at colliders.
Finally, in order to alleviate the Hubble tension two main ingredients are necessary. The first is a sufficient contribution to ∆N eff from the freeze-out of the Majorons. The main parameters controlling this are the mass of ϕ 1 , m ϕ 1 and its coupling to the Higgs λ 1H . Indeed, the Majoron is mainly aligned with the angular component of φ 1 and it can be kept in thermal equilibrium most efficiently via the mixing of ϕ 1 with the Higgs. In order to reach ∆N eff ∼ 0.4, the Majoron must decouple roughly with the muons, which can happen for m ϕ 1 < 1 GeV, as shown in Fig. 2. Regarding the coupling, the lack of evidence for an invisible Higgs decay at LHC requires λ 1H 0.01. On the other hand, λ 1H ≥ 10 −4 is necessary to keep the Majorons in thermal equilibrium until sufficiently late times. The second ingredient required to alleviate the Hubble tension is a coupling between the Majoron and the active neutrinos that allows  them to thermalize after BBN and a Majoron mass around the eV scale so that it will decay to neutrinos and contribute to ∆N eff at CMB. This decay width depends on the ratio of the neutrino masses over v 1 , as well as on the mass of the Majoron itself and the correct value is obtained for v 1 ∼ 100 GeV for m J ∼ 1 eV.

Conclusions
Extending the Standard Model particle content with right-handed neutrinos is arguably the simplest extension able to account for the evidence of neutrino masses and mixings. In order to also provide a natural explanation to the smallness of neutrino masses, two options emerge.
In the canonical, high-scale type-I Seesaw the ratio between the electroweak scale and the large Majorana mass provides naturally the required suppression. Conversely, in low-scale realizations, such as the linear or inverse Seesaw, the Lepton number L symmetry that protects neutrino masses is instead exploited. If this symmetry is approximate and only broken by small parameters, these will also naturally suppress the generation of neutrino masses. These low-scale realizations have the twofold advantage of enhancing the relevant phenomenological impact of the model and hence its testability, as well as avoiding a contribution to the Higgs hierarchy problem. We have explored the possibility that the small breaking of the L symmetry in the inverse Seesaw is dynamical. Its smallness emerges from a Seesaw-like structure in the scalar sector in which the vev of the L = 2 scalar responsible for neutrino masses is only indirectly induced by a vev around the electroweak scale, as in the type-II Seesaw. The parameter linking the two is small in a technically natural way since it is protected by an additional symmetry.
This spontaneous breaking of L in turn leads to the existence of a Majoron. We have explored the parameter space and conclude that this Majoron may contribute to the number of relativistic degrees of freedom in the early Universe as well as couple to the active neutrinos with the required values as to significantly alleviate the Hubble tension. This possibility is mainly constrained by the invisible decay of the Higgs, since the Majoron production critically depends on the mixing between the new scalar that breaks the L symmetry and the Higgs. Nevertheless, an order of magnitude smaller mixings than presently allowed by LHC constraints would still allow for a solution to the Hubble tension.
Among the new neutrinos introduced in low-scale Seesaws, two options exist due to the approximate L symmetry. The first are pseudo-Dirac pairs in which the left-handed component has a sizable mixing with the active neutrinos. The second are Majorana sterile neutrinos with couplings suppressed by the L-breaking parameters. For keV-scale masses, these Majorana sterile neutrinos may be sufficiently stable to be good dark matter candidates. While pro-duction via mixing through the Dodelson-Widrow mechanism is excluded by X-ray searches, we have shown that the correct relic abundance may be obtained for appropriate values of the model parameters via the freeze-in decays of the heavier pseudo-Dirac pairs to the new scalars and dark matter. These same couplings also control the mixing of the active neutrinos with dark matter as well as with the heavy pseudo-Dirac pairs. The main constraints on these mixings come from searches of the dark matter decays to X-rays and from unitarity tests of the PMNS mixing matrix from precision electroweak and flavour observables respectively. While the combination of these two probes rules out significant parts of the allowed parameter space, the correct relic abundance can still be obtained from the parameter that controls the mixing of dark matter with the right-handed component of the pseudo-Dirac pair. This mixing is more difficult to constrain, since the SM active neutrinos mainly mix to the left-handed component. Two possibilities are then viable. If the mixing between the active neutrinos and the heavy pseudo-Dirac pairs is sizable, close to their PMNS unitarity constraints, then the parameter that controls the dark matter-right-handed neutrino mixing needs to be significant. This implies a hierarchy of four or five orders of magnitude with respect to the mixing with the left-handed component, which may be considered fine-tuned. Conversely, if the two couplings are similar, the mixing of the heavy neutrinos with the active states needs to be very suppressed, reducing the testability of the model through PMNS unitarity deviations and, eventually, direct production at colliders.
Finally, the heavy pseudo-Dirac pairs could possibly explain the baryon asymmetry of the Universe through the ARS baryogenesis via leptogenesis mechanism [50,[70][71][72][73][74][75][76][77][78][79][80]. While we have assumed that the heavy pseudo-Dirac neutrinos thermalize, the ARS leptogenesis mechanism requires that at least some of them do not reach thermal equilibrium. This could be an option, since we only require one of them to thermalize in order to populate the DM abundance via its decays. Moreover, the correct DM density might also be obtained without thermalization of the heavy pseudo-Dirac pairs. This possibility together with the impact of the additional interactions of the heavy pseudo-Dirac pairs in the context of leptogenesis would be an interesting extension of the present study.
To summarize, we have shown that the SM extension considered with a dynamical breaking of the L symmetry characterizing the inverse Seesaw, is able to account simultaneously for the observed neutrino masses and mixings in a natural way as well as to provide a dark matter candidate with the correct relic abundance and alleviate the present Hubble tension between CMB and supernovae observations. The main constraints on the allowed parameter space come from unitarity tests of the PMNS mixing matrix through precision electroweak and flavour observables, searches for invisible Higgs decays at the LHC and X-ray searches for this decay mode of the sterile neutrino dark matter candidate. by the IFT Centro de Excelencia Severo Ochoa Grant SEV-2016-0597. This project has received support from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska -Curie grant agreement No 860881-HIDDeN.