A CMB Search for the Neutrino Mass Mechanism and its Relation to the $H_0$ Tension

The majoron, a pseudo-Goldstone boson arising from the spontaneous breaking of global lepton number, is a generic feature of many models intended to explain the origin of the small neutrino masses. In this work, we investigate potential imprints in the Cosmic Microwave Background (CMB) arising from massive majorons, should they thermalize with neutrinos after Big Bang Nucleosynthesis via inverse neutrino decays. We show that Planck2018 measurements of the CMB are currently sensitive to neutrino-majoron couplings as small as $\lambda \sim 10^{-13}$, which if interpreted in the context of the type-I seesaw mechanism correspond to a lepton number symmetry breaking scale $v_L \sim \mathcal{O}(100) \, {\rm GeV}$. Additionally, we identify parameter space for which the majoron-neutrino interactions, collectively with an extra contribution to the effective number of relativistic species $N_{\rm eff}$, can ameliorate the outstanding $H_0$ tension.

The majoron, a pseudo-Goldstone boson arising from the spontaneous breaking of global lepton number, is a generic feature of many models intended to explain the origin of the small neutrino masses. In this work, we investigate potential imprints in the Cosmic Microwave Background (CMB) arising from massive majorons, should they thermalize with neutrinos after Big Bang Nucleosynthesis via inverse neutrino decays. We show that Planck2018 measurements of the CMB are currently sensitive to neutrino-majoron couplings as small as λ ∼ 10 −13 , which if interpreted in the context of the type-I seesaw mechanism correspond to a lepton number symmetry breaking scale vL ∼ O(100) GeV. Additionally, we identify parameter space for which the majoron-neutrino interactions, collectively with an extra contribution to the effective number of relativistic species N eff , can ameliorate the outstanding H0 tension.
Introduction: Despite unambiguous evidence that at least two of the known neutrinos have a non-zero mass, the Standard Model (SM) is still lacking of an explanation of their origin. Perhaps more concerning, however, is the question of why neutrino masses are so much smaller than those of charged leptons. While many models have been proposed over the years to explain both the origin and smallness of the neutrino masses (see e.g. [1][2][3][4][5][6]), perhaps the most compelling class of models are those which invoke the so-called seesaw mechanism [7][8][9][10][11]. In such scenarios, the SM is augmented by heavy right-handed neutrinos carrying a Majorana mass term m N , which naturally give rise to light neutrino masses m ν of the order ∼ y 2 N v 2 H /m N , where v H 246 GeV is the vacuum expectation value of the SM Higgs, and y N is the Dirac Yukawa coupling of the right-handed neutrinos. Generating the Majorana mass term necessary to implement the seesaw mechanism is often accomplished by introducing a new scalar that spontaneously breaks lepton number. Assuming that lepton number is a global symmetry, as in the SM, the spontaneous symmetry breaking (SSB) triggered by the scalar leads to the prediction of a pseudo-Goldstone boson, the so-called majoron [12] (see also [13][14][15]).
The majoron is notoriously difficult to probe since it interacts very weakly with all SM particles, particularly with charged fermions λ φe ∼ 10 −20 [12]. However, measurements of the Cosmic Microwave Background (CMB) have reached a level of precision where small modifications to the neutrino sector may be discernible [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34]. The effect of including majoron-neutrino interactions in the early Universe are twofold [17]: (i) they lead to a non-standard expansion history after Big Bang Nucleosynthesis (BBN) and prior to recombination (generically amounting to ∆N eff ∼ O(0.1)), and (ii) they act to suppress the neutrino anisotropic stress energy tensor, and hence reduce neutrino free-streaming [16]. The idea of identifying features in the CMB arising from the majoron, and thus providing an indirect probe of the neutrino mass mechanism, was proposed at the start of the century [17]. However, until now, no rigorous cosmological implementation of this idea has been performed 1 , nor has there been an analysis using real data.
Using Planck2018 data [35,36], we analyze a wellmotivated region of parameter space in which majorons thermalize with neutrinos after BBN via inverse neutrino decay. We show that neutrino-majoron couplings as small as 10 −13 can be robustly excluded with existing CMB data; future experiments, such as the Simons Observatory [37] and CMB-S4 [38], which are aiming to probe the effective number of relativistic species N eff at the sub-percent level, could have sensitivity to couplings as small as 10 −14 . If interpreted in the context of the type-I seesaw model, these couplings point toward a lepton number symmetry breaking scale of O(100) GeV and O(1) TeV, respectively. Thus, quite remarkably, the CMB is providing an indirect probe of the neutrino mass mechanism at collider energy scales (albeit unaccessible to colliders due to their small couplings), but using feeble interactions with neutrinos in the early Universe. While the ΛCDM model has been incredibly successful at describing both high-and low-redshift cosmological observations, a concerning tension has recently emerged between the value of the Hubble constant H 0 inferred using early Universe observations (with data either from the CMB [35], or by combining measurements from BBN with baryonic acoustic oscillations, i.e. BAOs [39][40][41]), and various local late Universe measurements performed 1 Refs. [20][21][22] explored the possibility that some component of radiation contained strong self-interactions; this was accomplished by artificially setting to zero the multiples ≥ 2 in the Boltzmann hierarchy for the interacting radiation. This approach, however, cannot be applied (or mapped) into the scenario of [17], since neutrino-majoron interactions rates are strongly time-dependent and not infinite in strength.
using observations of type-Ia supernovae (see e.g. [42][43][44][45][46][47]) and strong lensing [48][49][50][51] (see e.g. [52] for an overview of the various measurements). The most prolific of these discrepancies is between the value inferred by Planck, H 0 = 67.4 ± 0.5 km/s/Mpc [35], and that observed by SH 0 ES collaboration, which relies on cepheids to calibrate the distance to type-Ia SN, who find a value of H 0 = 74.0 ± 1.4 km/s/Mpc [43]. Depending both on the choice of distance calibration and how one chooses to combine datasets, the outstanding tension is determined to be at the level of ∼ 4−6 σ [52,53]. While it is of course possible that this tension is a consequence of unaccounted for systematics in either or both measurements, throughout this work we will take this discrepancy at face value and assume alternatively that this is an indication of new physics beyond the ΛCDM paradigm. Various groups have attempted to resolve this issue by including additional contributions to N eff [54][55][56][57][58], strong neutrino self-interactions [30,31], hidden neutrino interactions [59][60][61], exotic dark energy models [62][63][64][65][66][67][68][69][70][71][72][73][74][75], dark sector interactions [76][77][78][79], and modified theories of gravity [80][81][82]. Most of these solutions are either incapable of resolving the tension fully [83][84][85], are experimentally constrained [86], are highly fine-tuned, or lack theoretical motivation. Perhaps the most simple, and thus theoretically appealing, solution which can ameliorate the H 0 tension to the level of ∼ 3σ is simply to postulate the existence of non-interacting dark radiation producing a shift in the radiation energy density relative to the value predicted in the Standard Model of ∆N eff ∼ 0.25. A more appealing, albeit far more problematic, solution was introduced in [30], where it was shown that strongly interacting 2-to-2 neutrino scatterings together with a contribution to ∆N eff ∼ 1 was able to fully resolve the tension; unfortunately, this solution requires neutrino couplings that are not phenomenologically viable [86], a value of ∆N eff excluded by BBN [87], and is only successful at reducing the tension if CMB polarization data is neglected. Given that the majoron naturally contributes to ∆N eff at the level of ∼ 0.11 via late-time thermalization and decay, and damps neutrino free-streaming in a manner similar to that of the strongly interacting neutrino solution, it is natural to ask whether 2-to-1 neutrino-majoron interactions are capable of further reducing the H 0 tension, beyond what is simply accomplished with ΛCDM + ∆N eff . Indeed we show that including majoron-neutrino interactions broadens the posterior such that the H 0 tension can be further reduced, albeit only to the level of 2.5 σ, a level that is comparable with other viable solutions, such as early dark energy (see e.g. [70]).

Majoron Interactions:
We parametrize the majoronneutrino interaction as: where ν corresponds to a light neutrino mass eigenstate. The coupling λ, taken here to be universal, is typically intimately related to the mass of the active neutrinos m ν and the scale at which lepton number is spontaneously broken, v L . For example, in the type-I seesaw mechanism, λ can be expressed as where U is the mixing between sterile and active neutrinos, and the last line follows from a condition in the type-I seesaw that U 2 ∼ m ν /m N [6]. Interestingly, for values of v L ∼ v H and neutrino masses consistent with current constraints, the value of λ within this model can e.g. naturally be of the order of 10 −12 , which happens to be around the region where inverse neutrino decays (νν → φ) can thermalize light majorons after BBN, but prior to recombination. In what follows we will treat λ as a free parameter to remain as model-independent as possible, and when appropriate, relate v L to λ by considering the atmospheric mass splitting m ν ∼ |∆m 2 atm | 0.05 eV [88]. Namely, v L 1 TeV (10 −13 /λ).
The Majoron Mass: Quantum gravity is expected to break all global symmetries [92,93], and hence the majoron should acquire a small but non-zero mass. Naively, one might expect the majoron mass to arise from dimension-five (dim-5) Planck scale suppressed operators [94,95]. Should these dim-5 operators involve only the Higgs and the scalar responsible for the SSB of lepton number, the majoron mass is expected to be where β is the coupling constant of a given operator at the Planck scale -which, for concreteness, we have assumed to be the same for all relevant dim-5 operators [95]. Of course, the actual details of the breaking of global symmetries by gravity depend upon the unknown quantum nature of the gravitational theory at the Planck scale [96]; thus we treat m φ as a free parameter in this work, centered approximately around the keV scale, but allowed to vary from 0.1 eV to 1 MeV. Fig. 1 contains a depiction of the majoron parameter space relevant for this work. In addition to highlighting parameter space currently excluded by SN1987A [90,91], BBN (see Supplementary Material, and e.g. [97]), and KamLAND-Zen [89], we identify masses and couplings for which the majoron is consistent with arising from dim-5 Planck scale suppressed operators 2 . We defer discussion of the remainder of this plot to later sections.
Model Extensions: Looking forward, it may be interesting to consider the possibility that one of the active FIG. 1. Majoron parameter space. The left and right vertical axes correspond to the majoron-neutrino coupling and the scale at which lepton number is spontaneous broken in the type-I seesaw model respectively. Current constraints from KamLAND-Zen [89], BBN (see text), and SN1987A [90,91] are shown in grey. The pink region demarcates parameter space for which the majoron fully thermalizes after neutrino decoupling, leading to ∆N eff = 0.11. The green band highlights the region of parameter space in which the majoron mass could arise from dim-5 Planck suppressed operators (3). Shown in blue is the parameter space excluded in this work using Planck2018 data at 95% CL. The parameter space below the black dotted line is excluded if there was a small but primordial population of thermal majorons. The region labeled 'H0' is the preferred 1σ contour for resolving the Hubble tension.
neutrinos is exactly massless, as this would decouple the lightest neutrino form the majoron, changing the cosmological evolution of the system. One could also conceive of the possibility of a multi-majoron system resulting from the SSB of a more complex flavor symmetry group in the neutrino sector [17]. In such a scenario, one could produce a more complicated thermalization history which produces step-like features in the evolution of the energy density, and damps the perturbations in a nontrivial manner. While these models are beyond the scope of the current work, they provide a clear extension of the ideas and prospects studied here.
Early Universe Cosmology: The collision terms governing the evolution of the neutrino and majoron phase space distributions are determined by the decay rate of the majoron into two neutrinos φ →νν, given by where in the last step we have considered m ν m φ . In order to model the time-dependent evolution of the number density and energy density of the system, we follow [98] (see also [99]) in assuming that all relevant species are characterized by a temperature T i and chemical potential µ i , and solve for their time evolution ac-counting for all relevant interactions 3 (see Supplementary Material for details). If the majoron is sufficiently heavy and interactions sufficiently strong, the majorons may begin to thermalize prior to or during BBN, leading to an enhanced expansion history of the Universe that would alter the formation of the light elements. For small couplings and masses (λ 10 −5 and λ 10 −10 MeV/m φ ), majorons thermalize with neutrinos after BBN, and when the majorons become non-relativistic at T ν ∼ m φ /3, they decay out of equilibrium to neutrinos leading to a small enhancement in N eff , which asymptotes to ∆N eff = 0.11. We identify in Fig. 1 a shaded pink region for which full thermalization is achieved after BBN. For yet smaller couplings, partial thermalization can be achieved; the dashed pink line in Fig. 1 identifies majorons that never thermalize, but augment N eff to a level that may be observable with CMB-S4 experiments [100].
We model the phase space perturbations by considering the coupled neutrino-majoron fluid, and approximate the entire system as being massless 4 . Despite the fact that the temperature of the Universe eventually becomes similar to the majoron mass, the majoron contribution to the energy density of the neutrino-majoron system is never larger than 10%. We have explicitly verified that the equation of state ω = (p φ + p ν )/(ρ φ + ρ ν ) and the speed of sound c 2 s = δ(p φ + p ν )/δ(ρ φ + ρ ν ) deviate by less than 3% with respect to that of an ultra-relativistic fluid, i.e. ω = c 2 s = 1/3 (see Supplementary Material). Additionally, we adopt the relaxation time approximation for the collision term [105], which has been shown to accurately reproduce the full solution in similar scenarios [26,27]. The above simplifications allow us to express the density contrast δ, the fluid velocity θ, the shear σ, and the higher anisotropic moments in the synchronous gauge as [105,106]: Here, h and η account for the metric perturbations, k is a given Fourier mode, F νφ represents the th multipole, a the scale factor, and Γ is the interaction rate accounting for inverse neutrino decays and majoron decays, given by where K 1 is the modified Bessel function of the first kind. For convenience one can approximate e µν Tν 1, and T γ /T ν 1.4 -we have verified that this introduces a negligible error in the final result. In Eqns. (5) all derivatives are understood to be with respect to conformal time.
Analysis: In order to efficiently scan the parameter space of interest, we define an effective interaction Γ eff in terms of the majoron mass and coupling as This effective interaction is defined such that for Γ eff 1 majorons thermalize in the early Universe. We perform runs with two distinct sets of priors: the first is used to place constraints on majoron models producing strong mann hierarchy is expected to be entirely negligible given current constraints on mν < 0.12 eV [35], see also [101][102][103][104].
modifications to the neutrino perturbations, and the second is used to identify parameter space for which the H 0 tension can be ameliorated. For both sets of runs, we adopt log-flat priors in λ or Γ eff and m φ spanning and respectively. In addition to these two parameters, we also allow for the possibility of extra relativistic and noninteracting degrees of freedom. We allow ∆N eff to vary linearly between −2 ≤ ∆N eff ≤ 4, and treat this additional radiation as free streaming. This additional contribution to N eff should not be considered ad hoc, but rather a natural expectation of majoron models. For example, should the reheating temperature be above the mass of right handed neutrinos, a thermal population of majorons produced in the early Universe may come to dominate the energy density of the Universe, producing nearly arbitrarily large contributions to ∆N eff . Such an effect becomes increasingly important for feeble interactions, such that an effective lower bound can be placed on the the neutrino-majoron interaction -needless to say, however, this bound is inherently dependent on pre-BBN cosmology. We include in Fig. 1 a line, labeled ∆N * eff , that identifies parameter space for which the contribution to ∆N eff from a primordial population of majorons would be excluded by Planck and measurements of large scale structure. We include a more comprehensive discussion of this effect in the Supplementary Material.
Results and Conclusions: After implementing the above modifications to both the energy density and neutrino-majoron perturbations in CLASS [107,108], we perform an MCMC with Montepython [109,110] using the Planck-2018 TTTEEE+lowlTT+lowE+lensing likelihood [36], including data on BAOs from the 6DF galaxy survey [111], the MGS galaxy sample of SDSS [112], and from the CMASS and LOWZ galaxy samples of BOSS DR12 [113], both including and excluding a Gaussian contribution to the likelihood on H 0 from SH 0 ES [43], taken to have a mean value and standard deviation of 74.0 and 1.4 km/s/Mpc. All MCMCs have been run until the largest Gelman-Rubin coefficient was R − 1 < 0.03 or better. In Table I we outline all relevant cosmological parameters for the analyses of Planck 2018+BAO+SH 0 ES data.
In Fig. 1 we show the 95% exclusion contours derived in this work, and the 1σ contour for parameter space preferred from the fit including the SH 0 ES likelihood. We derive the 95% CL exclusion contours using only Planck data in order to remain conservative, and note that including e.g. BAO data leads to a minor strengthening  Table I for best-fit values and 1σ uncertainties. The red solid line roughly corresponds to H0 = 68.0 ± 1.9 km/s/Mpc and hence is in 2.5σ tension with the SH0ES measurement.
of this contour. Interestingly, the results obtained here illustrate that Planck has already begun to significantly probe well-motivated regions of parameter space in which the majoron mass arises from dim-5 Planck scale suppressed operators. If interpreted in terms of the type-I seesaw model, current CMB observations are now probing lepton symmetry breaking scales O(100) GeV, with future CMB experiments potentially reaching the level of ∼ 10 TeV. Before continuing, we would like to emphasize that the constraints derived in this work are both stringent and robust over wide regions of parameter space. For example, a majoron of m φ = 1 eV and λ = 10 −11 is excluded by more than 5σ.
In Fig. 2 we show the posterior distributions for ΛCDM, ΛCDM with a floating value of ∆N eff , and the majoron + ∆N eff , each including and excluding the SH 0 ES likelihood. The SH 0 ES posterior is shown for comparison. Including the majoron broadens the posterior and induces a minimal shift of the central value to large H 0 , an effect which is more visible when the SH 0 ES likelihood is included. While the difference induced by including the majoron is not enormous, the H 0 tension can be reduced from 4.4 σ to 2.5 σ when neutrinomajoron interactions, and an additional contribution to dark radiation, are included.
By performing a MCMC including the SH 0 ES likelihood, we find that a scenario with ∆N eff = 0.52 ± 0.19, 0.1 eV < m φ < 1 eV, and coupling strengths λ ∼ (10 −14 − 10 −13 ) (eV/m φ ) -as highlighted in red in Fig. 1 -would render a posterior for H 0 of 71.9 ± 1.2 km/s/Mpc and an overall improvement of ∆χ 2 −12.2 with respect to ΛCDM. We remind the reader here that, because of the residual 2.5σ tension, it may not be entirely meaningful to combine the partially discrepant datasets, and thus care should be given in the interpretation of this region. Notice that the improvement in the χ 2 does not exclusively arise from the shift in H 0 ; this can be seen from the fact that the contribution of the CMB likelihood in the Majoron+∆N eff is less than that of ΛCDM. Interestingly, this region of parameter space corresponds to lepton number symmetry breaking scales in the type-I seesaw near the electroweak scale. Furthermore, it is worth emphasizing that unlike the strongly interacting neutrino solution proposed in [30] (defined by a 2-to-2 neutrino contact interaction), the solution proposed here is robust to the inclusion of polarization data, is phenomenologically viable, and is theoretically motivated.
An important comment on the consistency of this type of solution is necessary. If the contribution to N eff is of primordial origin, then successful BBN excludes values of ∆N eff 0.4 at T ∼ MeV [87,114,115]. In addition, including a floating value of ∆N eff in the CMB analysis can induce a shift in the preferred value of Ω b h 2 , which is also constrained by BBN. In the Supplementary Material, we address the extent to which the parameter space in the ∆N eff − Ω b h 2 plane preferred by the CMB fit is compatible with expectations of BBN.
Evidence for the existence of the majoron, arising from the spontaneous breaking of global lepton number, would provide a strong clue to the origin of the neutrino masses. In this work we have looked at the extent to which CMB measurements have probed the existence of such a particle through its impact on the expansion history of the Universe and its interactions with neutrinos. We show that there exists a broad range of well-motivated parameter space that is now excluded using Planck2018 measurements of the CMB power spectrum. Furthermore, we identify a region in which the majoron interactions help ameliorate the outstanding H 0 tension to a level that is beyond what is simply accomplished by including ∆N eff . If confirmed, the H 0 tension could be providing the first insight into the origin of the small neutrino masses.

Miguel Escudero and Samuel J. Witte
The Supplementary Material section contains additional information justifying various comments and statements asserted in the text, and outlining various computational details relevant for the reproducibility of this work. We also expand briefly on various phenomenological aspects. We begin by providing details on the computation of the background evolution and CMB phenomenology. We then discuss the derivation of the BBN constraint shown in Fig. 1, and finally discuss the implications of a primordial majoron population produced in the early Universe, which can be relevant should the reheating temperature of the Universe be sufficiently high.
Background Evolution: We follow [98] (see also [99]) and assume throughout that the distribution function for all relevant species can be characterized by their temperature T i and chemical potential µ i . The time evolution equations for each of such quantities reads [98]: where n, ρ, and p correspond to the number, energy and pressure density of the given species, H is the Hubble parameter, and δ t ρ and δ t n are the energy and number density exchange rates. Here the chemical potentials are the same for neutrinos and antineutrinos since they are produced at the same rates. Since we are exclusively interested in 1 ↔ 2 processes, within the Maxwell-Boltzmann approximation, we can express the energy and number density exchange rates as [98]: Conservation of energy and number of particles in the φ →νν process implies The above system of equations are solved 5 starting from a sufficiently large temperature such that the majoron population is negligible in the plasma, and with initial conditions obtained from neutrino decoupling within the SM [98]: This system is evolved until the maximum time between T γ = m φ /20 and t = 20 × τ φ to ensure that the majoron population disappeared from the Universe. We have ensured that the continuity equation dρ tot /dt = −H(ρ tot + p tot ) is fulfilled at each integration time step with a relative accuracy of 10 −5 or better. Should the majoron thermalize with the neutrinos while relativistic, occurring for τ φ 1/H(T = m φ /3), one can solve for the resulting temperature and chemical potential of the joint neutrino-majoron system. Imposing conservation of energy and number density: one finds that the equilibrium temperature and chemical potential are given by These values in turn imply that the maximum energy density contained in the majoron species is ρ φ (T eq , µ eq ) 0.09 × ρ ν (T eq , µ eq ) < 0.045 × ρ tot . (S10) In the case in which the majoron thermalizes with the neutrinos (i.e. τ φ < 1/H(T = m φ /3) equivalently to Γ eff > 1) and then decays, one finds the following asymptotic values for the temperature and chemical potential: Eq. (S11) can then be used to compute ∆N eff at the time of recombination and the energy density stored in neutrinos today (Ω ν h 2 ), the values of which are given by: In Fig. S1 we compare the interaction rate of a 1 keV majoron to the Hubble expansion rate for various values of λ. When Γ/H 1, the majoron equilibrates with the neutrinos. For the 1 keV candidate shown, this occurs for couplings λ 4 × 10 −12 , corresponding to Γ eff 1. The right panel of Fig. S1 illustrates the evolution of the energy density in the majoron system for the same 1 keV candidate. One can see both from the evolution of ∆N eff and from the evolution of the energy density that equilibrium is indeed attained for λ 4 × 10 −12 , as was expected from the simple comparison of the interaction rate. We have verified using the full solutions that the approximations adopted above are valid to high precision.
At the moment, it is not practical from a computational perspective to implement the evolution of the background for every sampled point in parameter space. To avoid this issue we derive fitting formulas to map the evolution of ρ νφ , valid for arbitrary values of m φ and Γ eff . These equations have been implemented into CLASS for a rapid evaluation of the background evolution. For the sake of reproducibility, we provide the fitting formula for the energy density of the neutrino-majoron system, expressed in terms of the majoron mass and Γ eff : where Γ eff is as defined in Eq. (7) and with: , (S14) (S17) CMB Phenomenology: For the numerical analysis presented in the main body of the text, we treat the neutrinos and majorons as a joint massless system. In order for this adopted treatment to be valid, the fractional shift in the energy density and equation of state from an ultra-relativistic system should be small. We illustrate in Fig. S2 that indeed this approximation holds to extremely high degree, thus validating the joint treatment of these two species within a single massless fluid. In Figs. S3 and S4, we illustrate the impact of the majoron on the TT and EE power spectra for m φ = 1 eV and various values of Γ eff (Fig. S3), and Γ eff = 10 4 with various values of m φ (Fig. S4). For sufficiently light majorons, presence of interactions enhances both the TT and EE spectra, and induces periodic oscillations in the C 's. For large masses, the impact of the perturbations vanish and the remaining signature is simply that induced by the presence of an additional contribution to ∆N eff . For completeness, we also show in Fig. S5 the relative change in the linear matter power spectrum induced at small scales for the same candidates.
As discussed in the primary text, including a floating value of ∆N eff shifts the preferred value of Ω b which is probed by the CMB and BBN. Obtaining a coherent cosmological picture requires ensuring compatibility of these two distinct probes with local measurements of H 0 . N eff can be modified after BBN and prior to recombination, as e.g. is done in the case of the majoron. Larger values of N eff , as preferred to resolve the H 0 tension, naturally shift Ω b to larger values, however the degeneracy of these parameters in the CMB and BBN is not exact for the case of ΛCDM + ∆N eff . This is shown explicitly in Fig. S6, and we note that the case of the majoron + ∆N eff is quite similar to the case of theΛCDM + ∆N eff . It is interesting that as the value of H 0 shifts toward the locally measured value (as occurs when one includes the SH 0 ES dataset in the likelihood), the preferred central value derived from the CMB analysis produces an increasing tension with the values inferred from BBN. While this tension is mild, it is important to bare in mind that the central value of H 0 in the Planck+BAO+SH 0 ES analysis is still reasonably below the central value preferred by the SH 0 ES data itself.
Finally, we show the two-dimensional posterior corner plot in Fig. S7. In the Γ eff vs m φ two-dimensional posterior a double peak structure in the majoron mass can be seen. This is a result of the fact that majoron-neutrino perturbations particularly affect the CMB spectra when m φ ∼ 1 eV (as can be appreciated from Figure S4 and from the posterior). This leads to stringent constraints on Γ eff for m φ ∼ 2 eV as highlighted in Figure 1, and to a double peak posterior on m φ with maximums at m φ ∼ 0.3 eV and m φ ∼ 30 eV.
Big Bang Nucleosynthesis Constraints: We set constraints on the neutrino-majoron coupling by requiring successful BBN by using the constraint on the effective number of relativistic degrees of freedom during BBN. The latest analysis finds [87]: This means that the one-sided 95% CL upper limit is N BBN eff < 3.33. Since we do not explicitly solve for the light element abundances within our modified cosmology, we conservatively adopt an upper limit of N BBN eff FIG. S6. 1σ and 2σ contours for Ω b h 2 and ∆N eff using the measured primordial element abudances [87] (black dashed), and compared with the preferred regions in a ΛCDM+∆N eff scenario obtained using the Planck+BAO and Planck+BAO+H0 likelihood analysis. induce a shift in N eff at the level of ∆N eff = N eff − N SM eff = 4/7 0.57 (where N SM eff = 3.045 [98,116,117]), which is clearly excluded by the measured primordial element abundances.
The main effect of a positive contribution to ∆N eff at the time of BBN is to induce a higher expansion rate of the Universe during the formation of the primordial elements with respect to the SM. Thus, the bound from BBN on ∆N eff can be interpreted as a time constraint on the generation of the primordial element abundances. In particular, in a Universe with ∆N eff = 0.45, the time at which Deuterium forms (corresponding to T D γ 0.07 MeV [118][119][120]) is t = 256.69 s 6 . Consequently, we derive the BBN constraint shown in Fig. 1 by requiring Any appreciable change in the expansion history of the Universe for T γ T D γ has been shown to render a minor impact on any relevant primordial nuclei abundance [115]. We find that imposing the constraint in Eq. (S19) leads to the following bound on the majoron-neutrino coupling strength: where the first term in the denominator results from the majoron production by inverse neutrino decays and the second term (i.e. λ < 4 × 10 −5 ) arises from the production of majorons via neutrino-neutrino annihilations (see [34]).
Primordial Majoron Abundance: Majoron interactions with the Standard Model arise through the active-sterile neutrino mixing [12], implying that majorons have non-negligible interactions with heavy sterile neutrinos. If the Universe was reheated to sufficiently high temperatures, it is conceivable that a primordial thermal population of majorons can be produced as a result of these interactions [95]. Here we comment under which conditions this occurs and the cosmological implications of such a primordial majoron population. Within the singlet majoron scenario, sterile neutrinos decay into an active neutrino and a majoron at a rate [3,121]: Which implies that sterile neutrinos will have sizable decays to majorons provided that v L > v H , or equivalently λ 10 −13 . The production rate of sterile neutrinos from the SM plasma is Γ ∼ 4 × 10 −3 y 2 N T [122][123][124], where y N is the sterile neutrino Dirac Yukawa coupling -which, within the type-I seesaw is y N ∼ 4 × 10 −8 m N /GeV m ν /0.05 eV. By comparing the Hubble parameter H ∼ 1.66 √ g T 2 /M pl with the sterile neutrino production rate, it is easy to show that sterile neutrinos (with couplings capable of generating the observed neutrino masses) are brought into thermal equilibrium at temperatures and would disappear from the plasma soon after they become non-relativistic, at T ∼ m N /3. Clearly, if such sterile neutrinos decay into majorons, they will produce a primordial thermal population of these particles. This statement is, however, dependent upon the unknown thermal history of the Universe -for example, this can be trivially avoided if the reheating temperature T RH < m N , as it would prevent sterile neutrinos from ever being thermalized in the early Universe.
Once sterile neutrinos decay/annihilate away from the thermal plasma, the majoron bath decouples from the SM model plasma. The majoron temperature after electron-positron annihilation is simply given by entropy conservation and reads: which corresponds to ∆N eff = 0.027 at the time of BBN, provided m φ < 1 MeV.
In Fig. S8 we compare the relative contribution to N eff as a function of Γ eff , assuming (i) a small pre-existing thermal population of majorons present at early times (that decoupled at T 100 GeV) and (ii) majorons are only produced via inverse decays of neutrinos. Projected sensitivity from the Simons Observatory [37] and CMB-S4 [100] are shown for comparison. For small majoron interactions, a pre-existing thermal population comes to dominate the energy density and produces a large shift in ∆N eff that can be easily constrained by observations of the CMB.
In Fig. 1, we include a black dotted line that denotes the region of parameter space for which a primordial majoron population would produce ∆N eff ≥ 1 at recombination, and would thus be excluded by Planck. To be concrete, we determine the bound at: Finally, notice that presence of a primordial population of majorons would lead to enhanced damping to the neutrino anisotropic stress and an additional contribution to ∆N eff (such that if Γ eff > 1, ∆N eff would asymptote to 0.16 rather than 0.11) -in this scenario, we expect that the constraints derived in this work to strengthen.