The Hubble Tension as a Hint of Leptogenesis and Neutrino Mass Generation

The majoron, a neutrinophilic pseudo-Goldstone boson conventionally arising in the context of neutrino mass models, can damp neutrino free-streaming and inject additional energy density into neutrinos prior to recombination. The combination of these effects for an eV-scale mass majoron has been shown to ameliorate the outstanding $H_0$ tension, however only if one introduces additional dark radiation at the level of $\Delta N_{\rm eff} \sim 0.5$. We show here that models of low-scale leptogenesis can naturally source this dark radiation by generating a primordial population of majorons from the decays of GeV-scale sterile neutrinos in the early Universe. Using a posterior predictive distribution conditioned on Planck2018+BAO data, we show that the value of $H_0$ observed by the SH$_0$ES collaboration is expected to occur at the level of $\sim 10\%$ in the primordial majoron cosmology (to be compared with $\sim 0.1\%$ in the case of $\Lambda$CDM). This insight provides an intriguing connection between the neutrino mass mechanism, the baryon asymmetry of the Universe, and the discrepant measurements of $H_0$.


I. INTRODUCTION
The Hubble Tension, circa early 2021. Despite its simplicity, the standard cosmological model (i.e. ΛCDM) has proven to be remarkably successful in describing the vast array of cosmological observations at hand. However in recent years, a growing discrepancy has emerged between the value of H 0 as predicted by ΛCDM [1][2][3][4][5][6], H 0 67.4 ± 0.5 km/s/Mpc and local observations that favor a significantly larger value, with a central values ranging between 70 H 0 74 km/s/Mpc, coming from e.g. type Ia supernovae [7][8][9][10][11][12][13][14][15], strong gravitational lensing [16][17][18][19][20], surface brightness fluctuations [21], and megamasers [22] (see e.g. [23,24] for recent reviews). This tension has now reached a significance quantified at the 4 − 6 σ level [24], and appears across an array of different datasets with seemingly independent systematics. It is thus necessary to now consider the very real possibility that this discrepancy is arising from a failure of the ΛCDM model to accurately describe the evolution of the Universe.
In this light, a large number of potential solutions have been proposed which typically fall into one of two categories: those which modify the Universe at late times (z 1) and those which modify the dynamics and evolution near recombination (10 3 z 10 5 ) -we refer the reader to [25][26][27][28][29][30][31][32] and to [33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49] for several recent proposals of each type, and to [50] for a recent comprehensive review of models. Cosmological observations from type Ia supernovae and baryonic acoustic oscillations (BAO) severely limit the viability of late-time solutions (see e.g. [51][52][53]). In order for early Universe solutions to be successful in raising the inferred value of H 0 , one must modify the expansion rate near recombination. While these models often struggle to maintain the high-quality fit to CMB data that is obtained in ΛCDM, this classification of solutions has proven thus far to be most successful in producing agreement between all cosmological datasets -see [53]. Finding meaningful ways to motivate the novel physics required, however, has been challenging -many of the proposed solutions require complicated and un-motivated models which are fine-tuned in order to ensure that the effects 'turn-on' at the correct epoch and produce sufficiently sizable shift in H 0 .
The Hubble tension and a light Majoron. Recently, in Ref. [54] (see also [55] for a short summary), the authors illustrated that a light majoron, naturally arising in neutrino mass models from the spontaneous breaking of a global lepton number symmetry [56][57][58][59], could partially counteract the effect of additional dark radiation (i.e. ∆N eff ), pushing the inferred value of H 0 to larger values while maintaining a good fit to the CMB data. The presence of the majoron has two effects [60]. First, the interactions damp the free streaming nature of neutrinos, which in turn damps the anisotropic stress; since the anisotropic stress sources the metric, the net effect is a time-dependent modification to the growth of potential wells [61]. Second, for large enough interactions, the majorons can thermalize with neutrinos between Big Bang Nucleosynthesis (BBN) and recombination. Upon becoming non-relativistic these majorons decay back into neutrinos, producing a net enhancement in ∆N eff [54,60]. Ensuring the effects of the majoron occur at the correct epoch requires a mass m φ near the eV scale; this scale, however, is not arbitrary, and can easily be motivated should global lepton number be explicitly broken by physics at the Planck scale [62,63], as might be expected in theories of quantum gravity [64][65][66][67][68][69]. In order to generate a sufficient level of damping one requires    At early times (high temperatures), a global U (1)L symmetry is spontaneously broken, generating sterile neutrino masses and giving rise to a pseudo-Goldstone boson: the majoron (φ). Sterile neutrinos start to be sizeably produced (but do not equilibrate) at T ∼ 10 6 GeV. Then, at T ∼ [10 6 − 10 4 ] GeV the CP violating oscillations of these sterile neutrinos generate a net primordial lepton asymmetry in the Standard Model. Soon after the electroweak phase transition (at T ∼ 130 GeV) sphalerons freeze-out and yield a final baryon asymmetry from the initial lepton asymmetry. After sphaleron freeze-out, sterile neutrinos and majorons thermalize with the plasma, and later decouple when sterile neutrinos decay. In particular, for ∼ GeV scale sterile neutrinos this occurs at temperatures below the QCD phase transition T 100 MeV. Finally, right before recombination, majorons with m φ ∼ 1 eV re-thermalize with active neutrinos (νν → φ) before decaying (φ →νν), generating a larger inferred cosmological value of H0.

j A F C W G T y z B R D F 7 K y I j r D A x N q C S D c F b f v k v a Z 3 V v H r t / K 5 e a V z l c R T h C I 6 h C h 5 c Q A N u o A k + E H i
neutrino-majoron couplings λ ∼ 10 −13 [54,55] 1 . This coupling, when interpreted in the context of the type-I seesaw favors a lepton symmetry breaking scale slightly above the electroweak scale (v L ∼ 1 TeV). Arguably, the only unmotivated aspect of this proposed solution is the apparent ad hoc contribution of ∆N eff , preferring values ∼ 0.5, which are in mild tension with BBN [73,74].
Primordial Majorons from Leptogenesis. In this work we attempt to source the additional dark radiation required to resolve the H 0 tension from a primordial population of majorons. We show explicitly that these particles can be produced from the decays of GeV-scale sterile neutrinos in the early Universe. Coincidentally, 1 The model discussed here has, on occasion, been confused with that of the strongly interaction neutrino solution proposed in [44,45]. In light of this, we take the opportunity here to highlight the many differences. First, the solution of [44,45] requires a neutrino self-interaction cross section 10 orders of magnitude larger than that present in the Standard Model. This, in turn, requires a new MeV-scale neutrinophilic boson with order one couplings. These values are not motivated in neutrino mass models, and are robustly excluded by experimental data unless the boson interacts only with τ neutrinos [70][71][72]. Next, the solution requires an additional contribution of ∆N eff ∼ 1, a value robustly excluded by BBN [73,74] -see also [75,76] for a recent assessment of the BBN bounds and [77,78] for models trying to evade these constraints. Finally, the observed shift in H 0 only occurs when polarization data is not included in the fit [44,[79][80][81], while the results for the majoron model discussed here are robust to the inclusion of this dataset. Thus, while the proposed models both involve neutrinophilic bosons, they are in fact remarkably different.
sterile neutrinos at the GeV scale are precisely those required for a successful implementation of low-scale leptogenesis via sterile neutrino oscillations, i.e. ARS leptogenesis [82] (see also [83][84][85]). We verify explicitly that symmetry breaking scales v L ∼ (0.01 − 1) TeV required to resolve the Hubble tension can be made fully consistent with conventional ARS leptogenesis, so long as the Higgs mixing is small enough so as to avoid thermalizing the scalar responsible for breaking lepton number, and that the lepton number phase transition occurs at T > 10 4 − 10 6 GeV. The scenario proposed here thus offers an intriguing connection between the H 0 tension, the neutrino mass mechanism, and the generation of the baryon asymmetry of the Universe. Fig. 1 shows a sketch of the thermal history, highlighting the main ingredients of our proposal.
This manuscript is organized as follows. We begin by introducing the well-known singlet majoron model in Section II. In Section III we first discuss the requirements in order to successfully produce the baryon asymmetry of the Universe via the ARS leptogenesis mechanism, and then compute the thermal evolution and subsequent decays of the sterile neutrinos responsible for sourcing the primordial majoron abundance. Section IV describes the cosmological evolution of the majoron-neutrino system, and presents the results of a MCMC performed using Planck2018 + BAO data. We present a summary and our conclusions in Section V. We finish in Section VI by discussing some interesting avenues for future work, and we refer the reader to the Appendices for various technical details.

II. THE SINGLET MAJORON MODEL
Throughout this manuscript we work with the wellknown singlet majoron model [56] (see also [59]) in which Majorana masses for the right-handed neutrinos N R are generated from the spontaneous breaking of a global U (1) L lepton number symmetry 2 . In this set-up, the small neutrino masses then arise from the type-I seesaw mechanism [88][89][90][91][92].
This model is realized by augmenting the Standard Model (SM) with n ≥ 2 right handed neutrinos with lepton number L = +1 and a scalar field Φ with L = +2, all singlets under the SM gauge group. This particle content and charge arrangements, together with the requirement of renormalizability of the interactions, leads to the following Lagrangian where λ Nij are the Φ-N Yukawa couplings, and h αi are the Higgs-Lepton-N Yukawa couplings. The scalar potential V Φ is given by Upon spontaneous symmetry breaking (SSB) of the U (1) L symmetry, the scalar Φ acquires a vacuum expectation value v L and will generate Majorana masses for the sterile neutrinos M N = λ N v L . Since U (1) L is a global symmetry, a pseudo-Goldstone boson appears on the spectrum: the majoron φ [56]. After SSB, it is convenient to parametrize Φ as where ρ is a CP even scalar, which in the limit λ ΦH → 0 will have a tree level mass given by m 2 ρ = 2λ Φ v 2 L . Only after the SSB of the electroweak symmetry will Dirac neutrino masses appear, [m D ] αi = h αi v H / √ 2 with v H = 246 GeV. Diagonalizing the neutrino mass matrix in the limit m D M N yields light active Majorana neutrinos with masses of the order: 2 Global symmetries are also expected to be explicitly broken by quantum gravity [64][65][66][67][68][69]. This has two important consequences. First, this breaking will generate a mass for the majoron; if the breaking is perturbative, one can show that dimension 5 Planck-suppressed operators naturally generate masses in the range 1eV m φ 100 keV [62,63], which coincidentally overlaps with the range of interest for the Hubble tension (see [54] for a discussion). Additionally, an explicit breaking of the symmetry guarantees that topological defects will naturally decay on short timescales [86,87], and thus consequently pose no threat of over-closing the Universe.
Assuming no strong cancellations occur, one typically expects the heavy sterile neutrinos to have very small mixings with light active states (see [93] for a general case), with a magnitude roughly given by and with Higgs' Yukawa couplings of the order One of the conventional appeals of the seesaw mechanism is that the fine-tuning of the Yukawa couplings required to generate the active neutrino masses can be ameliorated when M N > v H . Successful low-scale leptogenesis, however, requires GeV-scale sterile neutrino masses, which are capable of reducing, but not removing the aforementioned fine-tuning problem (see Eq. (6)).
While this naively appears to remove at least part of the original appeal, it is worth emphasizing that there exist mixing textures (i.e. non-trivial h αi ) 3 which can substantially enhance |θ| 2 relative to the value quoted in Eq. (5) (and further remove any need for tuning), however for the sake of simplicity and concreteness we choose to work within the prototypical seesaw limit.

A. Neutrino-Majoron Interactions
With Eq. (1) in hand, one can enumerate the novel interactions that arise between the neutrinos (both active and sterile), the majoron, and ρ. Working in the mass basis of N and ν, the relevant interactions in the seesaw limit are [95,96] where the couplings are given by 4 Here we have omitted generation indices, however this is valid when considering the three light neutrino mass eigenstates as λ ν is a diagonal matrix up to tiny O(|θ| 2 ) corrections [59]. It is worth highlighting that interactions between the majoron and charged fermions are both neutrino mass and loop-suppressed, and are thus expected to appear at the level of λ φee 10 −20 [56]. The majoron is thus, for all intents, a truly neutrinophilic boson. There are two interactions in particular that are worth mentioning at this point. In order to generate a sizable primordial population of majorons in the very early Universe (T > 10 MeV), we must require sterile neutrinos to efficiently decay into majorons, and in order to ensure majorons thermalize with neutrinos at late-times, we must ensure the inverse decay of active neutrinos to majorons exceeds the Hubble rate near recombination. As we shall see, both of these requirements are satisfied for values of M N ∼ GeV and v L ∼ v H .
Let us begin by assuming sterile neutrinos have thermalized at temperatures T M N , and identify the condition necessary to generate a sizable branching fraction to majorons. The decay rate of sterile neutrinos into a majoron and active neutrino is This should be compared with the decay rate into SM particles, which for sterile neutrino masses 0.1 GeV < M N < 10 GeV is roughly given by Γ(N → SM) ∼ 10 × Γ(N → 3ν) (see e.g. [97]), where Therefore, the ratio between the decay rates is implying a branching ratio Br(N → νφ) 1 for all relevant parameter space. We will show explicitly in the next section that in order for the majoron solution to remain viable (at least in the seesaw limit), we require v L to be 2 TeV and M N 3 GeV; the former is necessary in order to ensure neutrinomajoron interactions can damp the neutrino anisotropic stress, and the latter is required in order to generate a sufficiently large contribution to ∆N eff (if sterile neutrinos decay at earlier times the energy density of any majoron population produced will be diluted by subsequent entropy dumps). Eq. (13) then shows that sterile neutrino decays will be efficient in generating majorons across the entirety of the parameter space of interest. Now let us turn our attention to the late-time phenomenology, where we must require majorons to thermalize with active neutrinos near recombination. For the small couplings and masses of interest (λ ∼ 10 −13 , m φ ∼ 1 eV), this process proceeds via the inverse decays of active neutrinos (νν ↔ φ). The efficiency of this process is governed by the majoron decay rate intoνν: Notice that the interaction strength in this model is directly proportional to the light neutrino masses. Neutrino oscillation measurements imply [98,99]: 0.05 eV, depending of whether the neutrino mass ordering is normal or inverted, respectively. In [54], the authors considered the scenario in which the majoron interacted equally with all three neutrinos -which given the observed mass splittings is realized when m ν 0.15 eV. The two alternative limiting cases are realized when the lightest neutrino is approximately massless. In normal ordering this will correspond to m 1 ∼ 0, and then majorons interact almost exclusively with the most massive eigenstate m 3 (i.e. the number of interacting neutrinos is N int = 1). On the other hand, in inverted ordering, m 3 = 0 and then majorons interact with the two, nearly degenerate, neutrino eigenstates 1 and 2 (i.e. N int = 2). While it is far from obvious, we will show in what follows that cosmological observables are not strongly sensitive to the difference in N int , provided of course that N int ≥ 1.
For convenience in what follows, we will define here an effective width parameter where the normalization has been chosen such that for Γ eff 1 majorons thermalize with active neutrinos via inverse neutrino decays. This parameter is thus more intimately connected with the cosmological implications of majorons at late-times than the direct coupling λ ν itself.

III. ∆N eff AS A PRODUCT OF LEPTOGENESIS
The majoron solution to the H 0 tension, presented by the authors in [54], requires three novel ingredients: (1) a majoron with mass m φ ∼ eV (which could be generated by an explicit perturbative breaking of the U (1) L symmetry via dimension 5 Planck-scale suppressed operators [62,63]), (2) a coupling to active neutrinos at the level λ ∼ 10 −13 (corresponding to v L ∼ v H ), and (3) an additional contribution to ∆N eff ∼ 0.5. In this section we show that the novel contribution to ∆N eff , which had previously been introduced in an ad hoc manner, can be sourced directly from a primordial population of majorons. This occurs naturally if the masses of the sterile neutrinos are roughly M N ∼ GeV, which is coincidentally exactly the mass scale required for a successful implementation of ARS leptogenesis.
Here, we begin by outlining the main ingredients of the ARS leptogenesis framework, focusing in particular on whether any feature of the singlet majoron model could prevent or inhibit the generation of the baryon asymmetry of the Universe. We then argue that sterile neutrinos in the ARS leptogenesis framework inevitably lead to a thermal population of majorons after the electroweak phase transition, which subsequently decouples as sterile neutrinos decay, yielding a sizable majoron primordial abundance as relevant for CMB and BBN observations.

A. ARS Leptogenesis and the Majoron
The idea behind ARS leptogenesis is as follows [82] (see also [83,84], and [85] for a review): 1. One assumes that at sufficiently high temperatures (i.e. after reheating) there are no sterile neutrinos in the plasma. At temperatures T T EW ∼ 160 GeV an out-of-equilibrium population of sterile neutrinos is slowly produced via the small Dirac Yukawa couplings h αi .
2. Having been produced, these sterile neutrinos will undergo efficient CP-violating oscillations when t osc ∼ 1/H. For degenerate sterile neutrinos, this corresponds to temperatures where ∆M is the mass difference between a pair of sterile neutrinos.
3. These oscillations will generate lepton asymmetries in each of the sterile neutrinos individually, but in such a way that the total lepton number asymmetry is still zero. However, sphaleron processes will only convert the SM lepton asymmetries into a baryonic one. Then, when sphalerons freeze out at T ∼ 130 GeV, the baryon asymmetry present at that temperature is frozen and remains constant until today. Thus, as long as one sterile neutrino has not thermalized by T EW (such that there are non-vanishing SM leptonic asymmetries), a nonzero baryon asymmetry will have been generated.
Clearly, the combination of these steps meet the three Sakharov conditions and allow for successful baryogenesis. Rigorous calculations for the case of two sterile neutrinos have shown that the baryon asymmetry of the Universe can be successfully generated in the context of the seesaw limit for 0.1 GeV M N 10 GeV and ∆M/M N ∼ 10 −7 − 10 −5 , see e.g. [100][101][102][103]. For the case of three sterile neutrinos, similar calculations have shown that successful baryogenesis can be achieved without such strong mass degeneracy [104]. Importantly, ARS leptogenesis has not yet been rigorously investigated in the context of the singlet majoron model. To our knowledge, the only reference to have discussed this issue focused on identifying a minimal, but model-dependent, set of requirements for successful leptogenesis [105]. In order to ensure that the majoron model of interest here can indeed generate the observed baryon asymmetry of the Universe, we revisit the requirements identified in [105] using a relaxed set of assumptions.
There are three key requirements in order to maintain the efficiency of the ARS mechanism within the majoron model. Firstly, sterile neutrinos cannot thermalize with the Standard Model plasma at temperatures T > T EW , otherwise the lepton asymmetry (and thus also the baryon asymmetry) will vanish. Secondly, the sterile neutrinos must undergo CP-violating oscillations. Thirdly, such oscillations must be coherent at T lepto because it is then when the primordial lepton asymmetry is generated.
The thermalization of sterile neutrinos can occur via processes of the type φφ →N N , ρρ →N N and ρ → NN . If either of these scalar states have thermalized at high temperatures, then avoiding thermalization of the sterile neutrinos amounts to requiring small sterile neutrino couplings, or equivalently large vevs. In particular, one finds v L > 10 5 − 10 6 GeV [105]. Obviously this is in conflict with the requirement for majorons to re-thermalize with neutrinos near recombination, which requires v L 2 TeV. It is reasonable, however, to expect that these states would not have thermalized at early times since they are inherently a part of the sterile neutrino sector -which in ARS leptogenesis are assumed not to be produced during reheating. Should that be the case, we have verified explicitly that none of the processes mentioned above will generate a thermal sterile neutrino population at T > T EW , provided that |λ ΦH | < 10 −7 (see Appendix A).
In order for the CP violating oscillations of sterile neutrinos to be effective in the early Universe, the U (1) L symmetry should be broken at T c > T lepto ∼ (10 4 − 10 6 ) GeV. In the majoron model the sterile neutrino mass is a time dependent parameter controlled by the vacuum expectation value of the Φ field: , therefore one must also ensure that the U (1) L symmetry is spontaneously broken at T > T lepto (see Eq. (16)). By studying the 1-loop thermal corrections to the U (1) L potential, see Appendix A for the details, we have shown that the condition T c > T lepto ∼ (10 4 − 10 6 ) GeV can be translated into a bound on the Higgs-scalar mixing at the level of: Finally, within the ARS mechanism, the primordial lepton asymmetry is mainly generated when sterile neutrinos start to oscillate at T lepto ∼ 10 5 GeV. In this FIG. 2. Contribution from a primordial majoron population to ∆N eff at BBN (blue), assuming m φ 1 MeV, as a function of decoupling temperature T d (obtained using the entropy density in the SM from [108]). For majorons with Γ eff 1 and m φ 0.1 eV, the contribution to ∆N eff at recombination (red) can be greatly enhanced due to the fact that majorons do not decay immediately after becoming relativistic; this point is illustrated for Γ eff = 10 −1 and m φ = 10 (dashed). In green we show the approximate parameter space relevant for ARS leptogenesis, see Eq. (19). We show two constraints from BBN (dark and light blue horizontal regions; to be compared with blue solid line) which arise from two different determinations of Yp [109,110], see Appendix C.
stage, it is key that the coherence of such oscillations is maintained. In particular, processes of the type φN ↔ φN should not be faster than the oscillation rate Γ ∆M 2 /4E at T ∼ T lepto . Explicitly, see Appendix A, this requirement can be translated into a rather weak bound on the λ N coupling given by: which we can clearly appreciate implies a rather mild hierarchy between M N and v L . Collectively, these conditions imply that ARS leptogenesis in the context of the singlet majoron model will likely remain unaltered, so long as the Higgs-portal coupling is sufficiently small (|λ ΦH | 10 −7 ). While this requirement may appear at first sight to be tuned, we would like to point out that the smallness of this coupling is maintained by quantum corrections. Using SARAH [106,107], we have calculated the two-loop beta function for the running of λ ΦH , which has contributions in the The former of these is inherently small, and the latter is suppressed by the active-sterile neutrino mixings. Thus, the requirement |λ ΦH | 10 −7 is stable under radiative corrections.

B. Primordial Majoron Population from Sterile Neutrino Decays
Successful ARS leptogenesis requires sterile neutrino masses 0.1 GeV M N 10 GeV, which are generically expected to thermalize with the SM plasma at temperatures 1 GeV T 80 GeV [111]. For the sterile neutrino masses and vevs of interest, sterile neutrino annihilations will be efficient in thermalizing a majoron population during this epoch (note that Γ( Eventually, as the sterile neutrinos decay, the thermalized majoron population will decouple from the plasma and freeze out while relativistic. Comparing the decay rate to the Hubble expansion we can estimate the temperature at which majorons decouple (see Appendix B for details): which holds for sterile neutrinos with v L v H . In the event that majoron decoupling is instantaneous, one can compute the energy density stored in the majoron population at BBN by simply accounting for the entropy dilution in the SM plasma after decouplingthis is shown in Fig. 2 (blue line labeled 'Primordial') as a function of decoupling temperature T d . Here, we express the energy density in terms of and N SM eff = 3.044 [112][113][114][115][116]. It is worth noting that very light particles which decouple at extremely late times near T d ∼ 10 MeV may be in tension with constraints from Big Bang Nucleosynthesis (BBN). There is, however, some ambiguity as to where these constraints truly lie; the reason being that the leading local cosmology-independent determinations of the primordial helium abundance Y p rely on spectroscopic observations of HII in metal poor galaxies, which could suffer from systematics. The most recent estimate of Y P is [110]: Y P = 0.2453 ± 0.0034, but there are other independent analyses that report substantially larger values [109]: Y P = 0.2551 ± 0.0022. In order to account for the possibility of additional systematics in the determination of Y p , we present throughout constraints on ∆N eff from these two distinct analyses -we defer a more detailed discussion of the BBN systematics and the sensitivity to ∆N eff to Appendix C. These constraints are plotted in Fig. 2  Should sterile neutrinos decay at T T QCD ∼ 200 MeV, the effect of non-instantaneous majoron decoupling can lead to significant changes in the expected energy density of the primordial population. The reason here being simply that the SM plasma undergoes stronger and more rapid entropy dumps arising from the QCD phase transition. In order to more properly estimate the energy density in the parameter space of interest we have solved for the thermodynamic evolution of the N and φ number densities in the early Universe, and evolved this system to temperatures T ∼ 1 MeV for a wide array of sterile neutrino masses and interaction strengths. This was done using the methods of [112,117], in which one assumes each population can be approximately described by a thermal distribution with a time-dependent temperature and chemical potential. Following this approach, we show the contribution to ∆N eff as a function of M N and v L in Fig. 3. We also highlight the preferred parameter space for the H 0 tension (red region). The details of this analysis are contained in Appendix B.

IV. MAJORON COSMOLOGY
Having shown in the previous section that the same sterile neutrinos responsible for leptogenesis can also generate a sizable primordial population of majorons, we now turn our attention to the cosmological evolution of majorons after BBN. The story here is more straightforward: majorons with Γ eff 1 thermalize with neutrinos at temperatures T ∼ m φ viaνν ↔ φ processes, and subsequently decay into neutrinos after becoming nonrelativistic a short time later. There are two effects of this process: (1) a further increase in ∆N eff [54,60] and (2) a damping of the neutrino free streaming [60,61]. The former point has been illustrated in Fig. 2, where the contribution to ∆N eff at CMB (red) is compared for two different choices of parameters to the contribution to ∆N eff at BBN (blue). Assuming the majoron decays prior to recombination, models with Γ eff 1 lead to larger energy densities, as weakly coupled majorons remain non-relativistic for longer times before decaying. While the effect of damping neutrino free streaming is more subtle, these interactions can lead a noticeable imprint on the CMB multipoles. We outline our treatment of both of these effects in this section, and discuss the results of performing MCMCs to current cosmological data.

A. Background Evolution
In order to describe the background evolution of majorons and neutrinos in the early Universe we once again use the formalism developed in [112,117] (described above). This was shown in [112] to provide a highly accurate description of late-time thermalization of neutrinophilic bosons. We defer the details of these calculations to Appendix D.
Given an initial primordial majoron population characterized by ∆N BBN eff , we can solve for the evolution for the neutrino-majoron system by accounting for majoron decays and inverse neutrino decays φ ↔νν. Notice that this system is controlled by three parameters: the number of interacting neutrinos N int , the effective decay width Γ eff given in Eq. (15), and the majoron mass. We highlight in Fig. 4 the impact of late-time neutrinomajoron thermalization on the expansion rate of the Universe (via the contribution to ∆N eff ), for various choices of these parameters. From the upper panel of Fig. 4 we can see that i) the presence of neutrino-majoron interactions enhances ∆N CMB eff with respect to ∆N BBN eff , and that ii) for Γ eff < 1 the growth of ∆N eff is substantially larger (as mentioned above, this happens because for such values of Γ eff majorons decay out of equilibrium at T ν < m φ /3). Finally, from the lower panel of Fig. 4 we can see the impact of accounting for the three possible values of N int . The results are fairly similar, but perhaps contrary to what could be expected, ∆N eff is larger for smaller N int . This happens simply as a result of equilibrium thermodynamics: smaller values of N int lead to a more degenerate final state ν population (because in practice more neutrinos have gone into that sector), and therefore have a larger energy density.
An enhancement of the expansion rate prior and close to recombination has long been appreciated as a key ingredient in models attempting to resolve the H 0 tension. The reason being that the characteristic angular size of fluctuations in the CMB has been measured with very high precision, ∼ 0.03% [1]. For a flat FLRW Universe, this angular scale can be written as where z represents the redshift of last scattering, c s (z) is the photon-baryon sound speed, and H(z) is the expansion rate of the Universe. Given that θ s is measured to such high precision, and that c s (z) and z are constrained via alternative observables, the most clear modification that allows for larger values of H 0 is to enhance H(z) prior to recombination [53]. Using the approximate relationship between ∆N eff at recombination and H 0 derived in [118]: one can see that fully resolving the H 0 tension requires values of ∆N CMB eff ∼ 1. The problem is that values of ∆N CMB eff 0.3 are disfavored by Planck data, regardless of whether this radiation is dark and free-streaming [1] or strongly interacting [42]. Our set up, however, is completely different from either of these cases because the majoron-neutrino interactions are only efficient for a limited period of time, roughly between 3 m φ T ν m φ /10. For the mass range of interest, this implies majoron-neutrino interactions will only alter CMB multipoles 1000. This effect allows for an additional increase in ∆N eff relative to that of ΛCDM without spoiling the fit to the Planck observations.

B. Planck Analysis
To study in detail the effect of majoron-neutrino perturbations on the CMB we have modified the cosmological Boltzmann code CLASS [119,120]. In order to analyze the subsequent evolution of majorons as relevant for CMB observations we shall make a number of approximations, which we enumerate in what follows for the sake of clarity. First, we shall neglect neutrino masses in the evolution of the background energy density and at the level of the cosmological perturbations. We take this approximation because including neutrino masses at the perturbation level is rather complicated, see [121]. Nevertheless, this approximation is well justified provided that neutrino masses are at the level m ν 0.2 eV. This being said, large neutrino masses may be capable of reducing the value of σ 8 , as was the case for the strongly interacting neutrino solution (see [44]). Second, for the purpose of describing the neutrino-majoron perturbations we assume that they form a single coupled fluid. This is a good approximation because all species are effectively relativistic except when the majorons decay. Their contribution to the equation of state of the system is always below 13% for the models presented, and typically 5% for the parameter space of interest 5 . Finally, we approximate the collision term at the perturbation level by the relaxation time approximation [122]. This approach has been shown to be accurate for scenarios in which the interacting particles subtend large angles after collisions in the plasma [123]. In the parameter space we consider, majorons will only interact once they are mildly nonrelativistic, by which time the typical angular separation of neutrinos in the cosmic frame is large, θ > 10 • , which justifies our approach. Although it is beyond the scope of this paper to account for the exact collision term, we refer to [121] for a recent study dealing with the exact collision term in the context of invisible neutrino decays. With the exception of φBestFit, all models use ΛCDM best-fit parameters, adjusted to ensure Ω b , zmr, and H0 are constant. Notice that the best fit majoron cosmology has χ 2 φ − χ 2 ΛCDM = −4.5, implying a better fit than ΛCDM to Planck+BAO data.
In Fig. 5 we show how the various effects described here effect the TT power spectrum. Specifically, we isolate each effect, and show the residuals relative to that of ΛCDM for a majoron with m φ = 0.35 eV and Γ eff = 10 2 . These lines are produced by fixing Ω b , z mr , and H 0 to their ΛCDM values. We compare the background only contribution to a model in which ∆N eff is fixed to 0.4 to illustrate that the late-time thermalization only effects low multipoles. We also show the bestfit majoron cosmology for N int = 2 and T d = 50 MeV in red (∆N BBN eff = 0.37), which illustrates the high quality fit obtained from the MCMCs. From Fig. 5 we can appreciate that the neutrino-majoron interactions act as to partially cancel the background contribution at low multipoles.
To analyze in detail the cosmological implications of our scenario, we perform a MCMC using MontePython [124,125] including only Planck2018+BAO data [1,2] 67 . For each of our analyses we vary the standard cosmological parameters and nuisance parameters 6 Planck2018 data [1] includes the high-and low-(temperature and polarization) and lensing likelihoods [2]. In this work we use BAO data that includes the 6DF galaxy survey [126], the MGS galaxy sample of SDSS [127], and the CMASS and LOWZ galaxy samples of BOSS DR12 [128][129][130][131], as used in the fiducial Planck analysis. 7 We have also run an analysis including Pantheon type Ia supernova data, and have found that the preferred value of H 0 is surprisingly slightly shifted to yet higher values, and the ∆χ 2 with respect to LCDM is reduced. This emphasizes the robust nature of the obtained fit.
in the same way as the Planck collaboration in their legacy analysis [1]. For the majoron mass and interaction we use log-scale priors with the following ranges 8 : We note that the lower limit on m φ is imposed on the physicality of our model. It corresponds to the minimal mass for which the majoron can decay into the most massive neutrino m φ ≥ 2 |∆m 2 atm | = 0.1 eV. For smaller masses the majoron could potentially participate in the process of neutrino decay leading to a very different phenomenology [121,[132][133][134][135][136]. This, however, will in fact not occur within the framework of the singlet majoron model considered here, see [59].
In our runs we do not vary the initial primordial majoron population, i.e. ∆N BBN eff or equivalently T d , but instead we explore some representative values expected to arise from low-scale leptogenesis; this is because of the difficulty in computing the thermodynamic evolution on-the-fly, something the authors hope to improve upon in future work. We also explore scenarios with N int = 1 , 2 , 3 interacting neutrino species. Nevertheless, in what follows, we concentrate on N int = 2 for concreteness because our runs show that the posteriors are fairly independent of whether N int = 1 , 2 , 3.
In Fig. 6 we show the resulting posterior (see Appendix D for the full results) for analyses with N int = 2 and ∆N BBN eff = 0.37 in green, and ΛCDM in grey. For comparison, we also plot the one and two sigma posterior from local measurements of H 0 performed by the SH 0 ES collaboration, which find H 0 = 73.2±1.3 km/s/Mpc [15].
2. With respect to ΛCDM, there are ∼ 1σ upward shifts on other cosmological parameters such as n s , Ω b h 2 , and Ω cdm h 2 . Importantly, the observed shift in the preferred value of Ω b h 2 is in agreement with observations from BBN (see Appendix C for additional details).
3. The majoron mass needed to obtain these rather large values of H 0 is bounded from above to be: and is bounded from below at m φ > 2 |∆m 2 atm | = 0.1 eV. We note that these masses are in agreement with the expectations from Fig. 4. In addition, we note that this upper limit highlights that primordial populations of majorons with larger than eVscale masses that decay well before recombination are not favored by Planck legacy data. 4. The preferred region of parameter space for majoron interactions is Γ eff 0.03 at 2σ .
5. The Planck legacy data points to scales of spontaneous lepton number breaking: where this number is obtained by taking Eq. (10) with m ν = 0.05 eV, but we note that the bound on v L could be as large as 2 − 3 TeV for more massive neutrinos.
6. The best-fit point for of the fiducial model has a χ 2 of: and corresponds to These ∆χ 2 values highlight that for moderately large ∆N BBN eff the fit to Planck+BAO data can be improved relative to ΛCDM. However, our other analyses show that the fit becomes degraded for primordial populations with ∆N BBN eff 0.5. Notice that these values of ∆N eff are excluded at more than 2σ if the damping of neutrino free streaming is not included [1].

C. Implication for the Hubble Tension
In order to determine the significance with which the leptogenesis motivated majoron resolves the Hubble tension, we use the so-called posterior predictive distribution [137] (PPD) as proposed in [138]. Simply put, the idea is to determine the distribution of expected measurements obtained from new data (d ) given a set of observed data (d) and a model (I); in the context of the Hubble tension, this is equivalent to asking the question: given the value of H 0 observed by CMB and BAO data within a particular cosmology (taken here to be either ΛCDM or the primordial majoron cosmology), how likely is it for the SH 0 ES collaboration to measure some other valueĤ 0 CDL, obs ? The PPD is obtained by averaging the likelihood of the new data over the posterior of the parameters (θ) of the existing data, i.e. P r(d |d, I) = dθ P r(d |θ, I) P r(θ|d, I) . Importantly, there is no assumption that the two datasets under consideration are consistent. Rather, the PPD allows one to determine whether the second dataset is a likely draw from the first. We apply this method using the measurement of H 0 by the SH 0 ES collaboration, for which we adopt a Gaussian likelihood with central value 73.2 km/s/Mpc and uncertainty σ = 1.3 km/s/Mpc [15]. We plot the PPD for ΛCDM and the leptogenesis-inspired majoron model with N int = 2 in Fig. 7 (shaded regions). The central value of the SH 0 ES measurement is also shown for comparison (vertical yellow line). The potential tension between the CMB+BAO measurement and that of SH 0 ES can then be obtained using the PPD ratio, defined as whereĤ CDL, obs 0 is the value of the Hubble constant observed by SH 0 ES by using the Cepheid Distance Ladder (CDL). The PPD ratio can be interpreted as a lower bound on posterior probability of the hypothesis that the CMB+BAO data and that of SH 0 ES arise from the same value of H 0 without unaccounted for systematics. In the case of ΛCDM, we find a value of ρ PPD ∼ 10 −3 , implying the probability that the SH 0 ES measurement is unaffected by systematics and appears only from statistical fluctuations, is roughly ∼ 0.1%. In the majoron model, on the other hand, we find a value of ρ PPD ∼ 10%. Ad- mittedly, this metric does not account for the fact that additional model parameters have been introduced; nevertheless, the statistical tension between these datasets in the context of ΛCDM is sufficiently large that novel physics must now be considered, and thus it is perhaps more appropriate to compare the number of model parameters introduced only between proposed solutions to the Hubble tension.
A potential criticism of this treatment is that the decoupling temperature is fixed, rather than scanned, and thus has preferentially selected a large value of ∆N eff 9 ; the current treatment can be interpreted as adopting a strong prior on the ARS leptogenesis-motivated decoupling temperature. As shown in Appendix D, the ∆χ 2 of the best-fit is approximately equivalent to, and actually slightly below, that of ΛCDM. It is thus straightforward to understand that the effect of scanning T d would simply be to broaden the posterior on H 0 about a value of ∼ 70 km/s/Mpc, extending at lower values to what is found in ΛCDM 10 and potentially extending in the low decoupling temperature limit to a value H 0 ∼ 73 km/s/Mpc (as shown in [54]). Including a prior on H 0 from e.g. the SH 0 ES measurement, however, would shift the posterior to values of H 0 larger than that obtained here, and remove the part of the posterior near the ΛCDM value. 9 Note that implementing a scan of T d is quite complicated and time consuming, as the out-of-equilibrium evolution of the majoron must be solved on the fly in CLASS [119,120]. 10 The limit in which T d → inf, Γ 1, and m φ eV is as close as this model can be to recovering ΛCDM, however even in this limit energy is injected at the level of ∆N eff 0.12.

V. SUMMARY AND CONCLUSIONS
Majorons with ∼ eV masses can thermalize with neutrinos prior and close to recombination, damping neutrino free streaming, altering the energy density stored in the neutrinos themselves. These effects manifest in the CMB power spectrum in a manner that partially cancels, however there is a residual scale-dependent phase shift that leads to a preference for larger values of H 0 , n s , Ω b h 2 and Ω cdm h 2 . This motivated the authors' previous study [54], where it was shown that a majoron with a mass m φ ∼ eV arising from a the spontaneous breaking of a lepton number symmetry at scales v L ∼ O(1) TeV could reduce the Hubble tension to the ∼ 2σ level, however only if additional dark radiation was present at the level of ∆N eff ∼ 0.5.
In this work we have investigated the extent to which a primordial population of majorons arising from low-scale models of leptogenesis can source the additional radiation required to ameliorate the H 0 tension. We have found that so long as the Higgs' portal coupling is sufficiently small (|λ φH | 10 −7 ) so as to avoid thermalizing the new scalars at high temperatures and that the lepton number phase transition occurs at T > 10 4 − 10 6 GeV, ARS leptogenesis in the singlet majoron model will proceed as normal. Given the sterile neutrino masses required for ARS leptogensis and the range of interaction strengths required for majorons to thermalize at late times, sterile neutrinos will inevitably create a thermal majoron population at temperatures T M N . This primordial majoron population will eventually decouple as sterile neutrinos decay T ∼ M N /10, naturally sourcing a contribution of ∆N eff ∼ 0.4 (at BBN), see Fig. 3. Given that a large population of primordial majorons can be produced in the early Universe, the question becomes whether the late-time evolution of the majoron-neutrino system behaves similarly to the analysis of [54] (since the time-dependence of the interaction rate has been modified), and whether the result is sensitive to the number of interacting neutrinos. Our results indicate that all three scenarios (i.e. that in which one, two, or all three neutrinos interact) favor similar parameter space, provide a similar shift in H 0 , and yield comparable χ 2 values to CMB and BAO data. The preferred value of H 0 of these models using only Planck2018+BAO data is found to be ∼ 70.2 ± 0.6 km/s/Mpc, still below the value reported by the SH 0 ES collaboration 73.2 ± 1.3 km/s/Mpc [15], but compatible at the ∼ 10% level (to be compared with ∼ 0.1% in ΛCDM). It remains to be seen whether a more robust analysis including non-zero neutrino masses, both in the evolution of the background and at the level of the perturbations, would alter the results in a significant manner.
To summarize, we have shown that the Hubble tension could be a signal of low-scale leptogenesis and neutrino mass generation within a simple and well-motivated neutrino mass model. In this set up, lepton number is a global symmetry that is spontaneously broken at an energy scale v L . This generates Majorana masses for sterile neutrinos, which in turn via the type-I seesaw lead to an understanding of the small neutrino masses. In addition, upon breaking of U (1) L a pseudo-Goldstone boson appears on the spectrum, the majoron φ. This particle is naturally very light and eV masses for it are motivated from Planck-scale suppressed operators that explicitly break lepton number, m φ v L v L /M Pl . Furthermore, in this set up, the interactions between neutrinos and majorons are extremely feeble, λ ∼ m ν /v L ∼ 10 −13 . Intriguingly, these very small couplings and masses precisely correspond to φ ↔νν processes turning on right before recombination for v L ∼ v H . Indeed, our Planck legacy data analysis shows that the preferred region of parameter space within a primordial majoron cosmology is: which point towards Finally, given that the preferred scale of lepton number breaking is 2 TeV this naturally points to GeV-scale sterile neutrinos. In particular, these sterile neutrinos can produce the right primordial majoron population (see Fig. 3) and also generate via their oscillations the sufficient leptonic CP asymmetry in the early Universe. Therefore providing an understanding of the observed asymmetry between matter and antimatter.

VI. OUTLOOK
We have presented a connection between the Hubble tension, leptogenesis and the origin of neutrino masses. In what follows, we enumerate a list of aspects that fall beyond the scope of this study but that we believe deserve future attention: 1. Refined cosmological analysis. As discussed in the main text, we have made three approximations which simplify the treatment of the majoron cosmology: i) We have neglected the effect of neutrino masses in the evolution of both the background and the perturbations of the Universe. This is well-justified for m ν 0.2 eV as Planck CMB observations are not sensitive to such small neutrino masses. However, an actual analysis could potentially reveal the preference of higher neutrino masses which will yield a smaller value of S 8 as relevant for weak lensing measurements of this parameter, see e.g. [139]. ii) We have modeled the majoron-neutrino perturbations by making use of the relaxation time approximation for the collision term. As argued in Section IV we believe that this should be a good approximation because in our parameter space of interest majorons interact with neutrinos once they become mildly nonrelativistic, implying that the final state particles will be roughly isotropically produced (as assumed by the relaxation approximation to capture the effect on the anisotropic stress). It would, however, be very interesting to explore the phenomenology accounting for the exact collision term. iii) We have treated the neutrino-majoron system as a single coupled fluid. While this assumption holds in the relativistic limit, the isotropisation of the fluid decreases as majorons become non-relativistic and decay. It is thus possible that relaxing this assumption, and properly incorporating the decay hierarchy, could potentially reduce the maximum allowed value of ∆N eff . We hope to refine this treatment, and properly assess the importance of this effect, in future work.
2. ARS Leptogenesis in the singlet majoron model. In this work (see also [105]) we have derived a minimal set of conditions necessary for the singlet majoron model not to spoil successful ARS leptogenesis. In order to robustly determine the parameter space of interest more rigorous calculations would be necessary 11 . We note that given the temperature dependence of M N is this model, such analysis may yield even more favorable conditions for the production of a lepton asymmetry from sterile neutrino oscillations.
3. Collider Detectability. Although in this study we have exhausted the cosmological implications of majorons and their companion neutrinos and sterile neutrinos, we have not discussed potential signals at laboratory experiments. From the collider perspective, it appears quite difficult to test the scalar sector of the theory given the smallness of the Higgs portal coupling and the smallness of the neutrinomajoron couplings. However, in the context of sterile neutrinos there are some possibilities. In our set up, sterile neutrinos decay invisibly into a neutrino and a majoron (see Eq. (13)). This means that 11 We note that baryogenesis/leptogenesis has been explored in the majoron model in the context of electroweak baryogenesis [140,141], thermal leptogenesis [142], and resonant leptogenesis [143]. In addition, we note that the CP violating decays of Higgs doublets into sterile neutrinos in the early Universe can yield relevant lepton asymmetries that could (depending upon the mass degeneracy) dominate over the contribution arising from oscillations [144,145].
typical searches at beam-dump experiments (see e.g. [146]) will not have sensitivity to this model. The best avenue to detect these GeV-scale sterile neutrinos may be to look for K/π → N [147,148] decays, where the N particles appears in the form of missing energy. While at the moment current experiments are sensitive to active-sterile neutrino mixings ∼ (1 − 2) orders of magnitude larger than in the naive seesaw limit (see Eq. (5)), ongoing and upcoming experiments looking for these decays [149] may be sensitive to the minimum mixing expected from the seesaw for m N < m K/π − m µ, e .
In this work we have shown that the Hubble tension can be largely ameliorated in a simple framework that explains both the origin of the active neutrino masses (via the seesaw mechanism) and the baryon asymmetry of the Universe (via the ARS leptogenesis mechanism). This proposal may be exhaustively tested by both future cosmological observations and by looking in terrestrial experiments for the presence of GeV-scale sterile neutrinos, which are necessary to source both the primordial majoron population and a primordial lepton asymmetry. where once again we have made use of the fact that dξ dT = − ξ T Γ H , where in this case Γ 5 × 10 −3 |h| 2 T 12 where |h| 2 ≡ α |h αi | 2 , and we have normalized it to the seesaw limit case with M N = 1 GeV, see Eq. (6).
It is important to note that Eqs. (A1) and (A3) do not necessarily give the true evolution of the number densities of N , ρ and φ at very high temperatures. The reason being that we have not accounted for the possibility that processes such as ρ ↔ φφ,N N ↔ ρρ, etc. change the relative number densities of each of these species. Nevertheless, the total number densities of all of these states will remain unchanged, and thus Eqs. (A1) and (A3) represent upper bounds on the number densities. If interactions between these particles are efficient, all species will have roughly equivalent abundances, with the value give by the maximum of Eq. (A1) and Eq. (A3).
Since we know how ξ N and ξ ρ scale as a function of temperature, we can now study the evolution of the U (1) L phase transition. The 1-loop effective potential for the U (1) L sector is given by [151] V =m 2 ρ | eff (T ) where we have defined ρ ≡ ρ + v L . In writing this expression we have neglected logarithmic contributions and terms linear in T , which are inherently small for the parameter space of interest -as such, one can see that the U (1) L phase transition is of 2nd order. Here m 2 ρ | eff (T ) is the effective thermal mass of the ρ scalar which receives contributions from the three self-energy diagrams in Fig. A1. The contribution of each to the effective mass is Summing over all the contributions and including the tree-level value we obtain In order to ensure sterile neutrinos oscillations proceed as normal we must ensure M N = 0 at T lepto ∼ 10 4 −10 5 GeV. This amounts to requiring the phase transition to occur at T c > T lepto , where T c defines the temperature at which thermal corrections are subdominant to the tree-level value (such that the symmetry is broken by the vacuum). For each term in the thermal potential, we then derive constraints ensuring the aforementioned condition is met; this 12 This rate is calculated using [150], and taking the typical neutrino energy to be E ∼ 2T ; collectively this implies process gives We have differentiated here between the case in which the number density of N s is equal to that expected from the Dirac Yukawas and from the Higgs portal coupling. A quick inspection shows that Eqs. (A7d) and (A7e) are trivially satisfied across all of the parameter space of interest (see Fig. 3), that Eq. (A7c) is redundant with the condition required to avoid thermalizing the scalar sector, that Eq. (A7b) is satisfied by the mixings in the seesaw limit (see Eq. (6)), an that Eq. (A7a) poses the only relevant requirement, which matches our Eq. (17) in the main text.
Finally, in order to asses the 3rd requirement, we must ensure that processes of the type N φ ↔ N φ are not efficient at T lepto in order to ensure that the CP violating oscillations of sterile neutrinos are coherent. The rate for N φ ↔ N φ processes can be estimated to be Γ = n N σv λ 4 N /(576π) T ξ N . Comparing this rate with H we obtain the following bound on the λ N coupling: This requirement, thus, implies a mild hierarchy between M N and v L . Therefore, if sterile neutrinos, ρ's and φ's are not populated during reheating, and |λ ΦH | < 10 −7 and Λ N 0.07, we can conclude that i) the U (1) L symmetry is spontaneously broken before the onset of sterile neutrino oscillations, and ii) ρ's, φ's and N 's do not thermalize with the SM plasma prior to the electroweak phase transition, and iii) the CP violating oscillations will be coherent at the time at which the lepton asymmetry is generated. This confirms that ARS leptogenesis can remain viable in the singlet majoron model with symmetry breaking scales v L 1 TeV.
Admittedly, we have not carried out a rigorous calculation of leptogenesis within this framework. What our discussion shows is that in the parameter space of interest, our model converges to the conventional models in which ARS leptogenesis has been shown to be successful. It may be possible to further relax some of the identified requirements and still generate the baryon asymmetry of the Universe, however this requires detailed calculations which are beyond the scope of this work. We note that since the temperature of the phase transition can be adjusted, it could be chosen so as to potentially enhance the CP oscillation rate over some finite range of temperatures. Naively, this may enhance the primordial lepton asymmetry and broaden the parameter space in which ARS leptogenesis can be successful.

Appendix B: Primordial Majorons from GeV-scale Sterile Neutrinos
In this appendix, we outline the calculation detailing the creation and decoupling of the primordial majoron population. As discussed in the main text, sterile neutrinos with 0.1 GeV M N 10 GeV responsible for generating the active neutrino masses generically thermalize with the SM plasma after the electroweak phase transition [111]. Let us begin with a rough estimation to identify the relevant evolution after sterile neutrino thermalization, and then we will return to a more quantitive assessment.
At high temperatures sterile neutrinos annihilate efficiently to majorons, and will quickly generate a thermal population of them. The annihilation cross section (N N ↔ φφ) responsible for thermalizing majorons at high temperatures is given by where we have consider the limit m ρ M N for simplicity. By comparing the rate Γ ∼ n N σv with H, one can see that at temperatures near the sterile neutrino mass, Γ(N N ↔ φφ)/H ∼ 300 (M N /GeV) 2 (2 TeV/v L ) 4 . This implies thermal equilibrium will be achieved if v L < 8 TeV M N GeV 3/4 , which is valid in all of the parameter space of interest. At lower temperatures this rate will fall below the rate for sterile neutrino decay N → φ ν. The sterile neutrino abundance will be strongly depleted by T M N /10 (since inverse decays are no longer efficient in restoring the population), and thus the primordial majoron population will decouple near this epoch. We can estimate the approximate energy density stored in the majoron population at the time of BBN by identifying the temperature at which inverse decays ν φ → N are no longer efficient in altering the majoron distribution function. Since the typical energy exchanged by this process is E ∼ M N , this will roughly occur when ( , which yields an estimate of the decoupling temperature given by: Thus for v L v H , the decoupling of majorons will occur at T d M N /13. Provided that M N is sufficiently small (M N 1 GeV), the decoupling will occur after the QCD phase transition, which implies the a sizable primordial majoron population will be generated. In order to highlight the relevant phenomenology, in Fig. A2 we show the rates for each of the relevant processes as a function of temperature for M N = 1 GeV and v L = v H . We can clearly appreciate that the last rate to drop below Hubble in this example is ν φ → N at T ∼ 100 MeV. This correspondingly predicts 13 a value of ∆N BBN eff 0.34, which is in excellent agreement with the result from our computation (described below), which gives ∆N BBN eff = 0.37. In order to quantitatively study this evolution of the majoron and sterile neutrino populations in more detail, we again have chosen to model the thermodynamic evolution following the formalism developed in [112,117]. This formalism assumes that sterile neutrinos and majorons are both described by thermal equilibrium distributions with evolving temperatures and chemical potentials. This is a good approximation in this context of this problem because: i) majorons are massless, ii) sterile neutrinos have thermal abundances, and iii) sterile neutrinos start to decay while relativistic provided that v L < 2 TeV. The last requirement can be clearly seen by comparing the decay rate with the Hubble at T = M N /3: where we have normalized the number with respect to g ∼ 50, as would be relevant for T ∼ 200 MeV. In this framework, the time evolution of the temperature and chemical potentials of the relevant species reads: where in these expressions ρ i , n i , and p i are the energy density, number density, and pressure of the given species i. In addition, c SM ≡ dρ SM /dT is the heat capacity of the SM plasma that, together with ρ SM and p SM , we take from [108].
In the Maxwell-Boltzmann approximation, we have analytic expressions for the energy and number density exchange rates for decay N ↔ ν φ processes which read: where energy and number density conservation in the decay process implies that: In addition, in the Maxwell-Boltzmann approximation, the rates for annihilationsN N ↔ φφ given Eq. (A1) read: where for the sake of simplicity we have taken σ(s) M 2 N /(64πv 4 L ). This is a good approximation for all temperatures of interest.
The final ingredients needed to obtain an estimation of the size of the primordial majoron population are the initial conditions. Given that i) sterile neutrinos thermalize after the electroweak phase transition with the SM plasma, and that ii) annihilations between sterile neutrinos and majorons are highly efficient across our relevant parameter space, we start with initial conditions corresponding to all species in thermal equilibrium. In particular, we use: Contours correspond to 1 and 2σ CL. We have calculated the predictions using the results of [156], the deuterium measurements from [159], and performed two analyses for the helium measurements of [110] and [109]. In red we show the 1 and 2σ CL posterior from our CMB analyses in Section IV. We can appreciate that our region of interest is in agreement within 2σ with successful BBN irrespectively of the adopted YP value.
with t 0 = 1/(2H(T γ )). Note that for numerical stability we choose these numerical values for the chemical potentials, but that the results are equivalent to choosing µ N = µ φ = 0 since n ∝ e µ/T and therefore a ratio 1 : 100 has a negligible impact on any relevant thermodynamic quantity. The results of solving these equations are shown in Fig. 3 where we display ∆N BBN eff as a function of M N and v L .

Appendix C: Big Bang Nucleosynthesis
The requirement of successful BBN yields relevant constraints on primordial populations of majorons as parametrized by ∆N eff . At present, the two primordial abundances used to constrain non-standard expansion histories at the time of nucleosynthesis are [152,153]: Helium-4, Y p , and Deuterium, D/H. On the one hand, Y p is very sensitive to the expansion history of the Universe because its abundance is mainly controlled by the time at which deuterium starts to form, which corresponds to T D 0.073 MeV [154]. In addition, Y P is only logarithmically sensitive to the baryon abundance, Ω b h 2 , and its prediction has a negligible theoretical uncertainty. On the other hand, the deuterium abundance is strongly dependent upon the baryon energy density while only moderately dependent upon the expansion rate, potentially modified by ∆N eff . Importantly, although the recent results from the LUNA collaboration [155] have reduced the theoretical prediction for deuterium, it is still at the ∼ 2.8% level [156] -see also [157] and [158] which report 1.5% and 4.4% uncertainties, respectively. In order to understand the effect of these constraints in our parameter space we have used the predictions and theoretical uncertainties from the recent analysis of [156]. We contrast these predictions to the measured deuterium abundance [159]: D/H = (2.527 ± 0.030) × 10 −5 , and choose to do an analysis for two values of Y P . One from the recent analysis of [110] that yields Y P = 0.2453±0.0034, and the one from [109] that yields Y P = 0.2551 ± 0.0022. The reason we choose these two values is because although most recent determinations of Y P agree within error bars with [110], the determination of Y P is far from trivial and could have systematic uncertainties. Therefore, we consider the two to highlight the size of potential systematic uncertainties. This being said, we believe that it is likely that the real value is closer to that of [110].
In Fig. A3 we show the resulting constraints from successful BBN on ∆N eff as a function of Ω b h 2 when considering the determination of Y p from [110] (blue) and [109] (purple). In red we show the 1-2σ CL posterior from our Planck+BAO analysis within our benchmark majoron cosmology with a fixed ∆N BBN eff = 0.37 (see Section IV). We can clearly appreciate that the region of parameter space in our benchmark is compatible within 2σ with both values of Y P from [110] and [109]. We note that the fact that Ω b h 2 within the majoron cosmology is shifted upwards with respect to ΛCDM is relevant since it leads to a better agreement with the measured deuterium abundance. Finally, from this figure we can appreciate that given that Ω b h 2 cannot be too different from the ΛCDM value, a very conservative constraint within our cosmology would be ∆N BBN eff < 0.7.
Here, derivatives are taken with respect to conformal time, h and η represent for the metric perturbations, k is defines the given Fourier mode, F νφ represents the th multipole, a is the scale factor, and where Γ φ is the decay width at rest of the majoron, K 2 is a modified Bessel function of order 2, and since chemical potentials are small we effectively approximate e µν Tν 1. We implement these perturbations and the background evolution in the Boltzmann code CLASS [119,120], and we run an MCMC using Montepython [125] using Planck2018+BAO data on the leptogenesis inspired models with N int = 1, 2, 3 and T d = 50, 30 MeV (which correspond to ∆N BBN eff = 0.37, 0.48, respectively). Chains are run until fully converged and all Gelman-Rubin coefficients are ≤ 0.08. We present here triangle plots in the full parameter space (including derived parameters such as σ 8 ); the fiducial model, i.e. the one in which we take T d = 50 MeV to be fully consistent with all constraints from BBN, is shown in Fig. A4, and the result for N int = 3 and T d = 30 MeV is shown in Fig. A5. Both plots contain one and two sigma contours from the SH 0 ES measurement in grey. For completeness, we also include a table (Tab. A1) describing the best-fit values and the one sigma uncertainties.