Neutrino Lines from Majoron Dark Matter

Models with spontaneously broken global lepton number can lead to a pseudo-Goldstone boson as a long-lived dark matter candidate. Here we revisit the case of singlet majoron dark matter and discuss multiple constraints. For masses above MeV, this model could lead to a detectable flux of monochromatic mass-eigenstate neutrinos, which have flavor ratios that depend strongly on the neutrino mass hierarchy. We provide a convenient parametrization for the loop-induced majoron couplings to charged fermions that allows us to discuss three-generation effects such as lepton flavor violation. These couplings are independent of the low-energy neutrino parameters but can be constrained by the majoron decays into charged fermions.


I. INTRODUCTION
The observation of neutrino oscillations has raised the question why neutrino masses are so much smaller than all other known masses. The most-studied solution to this puzzle comes in the form of the seesaw mechanism [1], where heavy right-handed neutrinos suppress active-neutrino masses with respect to the electroweak scale. As a bonus, these heavy neutrinos can dynamically generate a baryon asymmetry in the early Universe via leptogenesis [2], thus solving a further problem of the Standard Model (SM). An inherent feature of the seesaw mechanism is the self-conjugate Majorana nature of neutrinos, which implies that the anomaly-free global U (1) B−L symmetry of the SM has to be broken by two units. Following the success of spontaneous symmetry breaking in particle physics, one can easily imagine that also this B −L symmetry is broken spontaneously, resulting in a Nambu-Goldstone boson named majoron [3,4]. Gravitational or explicit breaking terms then typically generate a mass term, making the majoron a pseudo-Goldstone boson. Since the majoron has couplings that are suppressed by the B − L breaking scale, i.e. the seesaw scale, it can easily be long-lived enough to form the dark matter (DM) of our Universe [5][6][7][8][9][10][11][12]. The CP-even partner of the majoron can on the other hand be used to drive inflation [13,14], although this is not the focus here. Since the seesaw scale is far above the electroweak scale for DM stability reasons, leptogenesis will be hardly modified by the majoron [15].
The most salient and well-studied indirect-detection signature of majoron DM is its decay into two photons, which most prominently arises in cases where one identifies the majoron with the axion [16][17][18][19], for which new particles have to be introduced to create a color anomaly of the U (1), often accompanied by an electromagnetic anomaly as well [20]. In this article we focus on other possible signatures, most notably from the tree-level decays into neutrinos and from the one-loop decays into charged fermions [7][8][9]. To this effect we provide a simple parametrization of the majoron couplings to the three * Electronic address: Camilo.Alfredo.Garcia.Cely@ulb.ac.be † Electronic address: Julian.Heeck@ulb.ac.be generations of fermions that also allows us to discuss lepton flavor violation (LFV) and perturbativity. Owing to its tree-level coupling, the key feature of majoron DM is arguably its two-body decay into monochromatic neutrinos, a topic that has received a lot of attention in recent years in its own right [21][22][23][24][25]. In fact, since neutrinos are the least-detectable SM particles, any limit on their flux automatically provides a model-independent lower bound on the DM lifetime [26]. (The same arguments apply to DM annihilations [27][28][29].) Majorons are a well-motivated DM candidate that can lead to observable monochromatic neutrino fluxes for energies between MeV and 10 TeV. For energies below the electroweak scale, these neutrino lines do not receive Bremsstrahlung corrections that could otherwise lead to observable gamma-ray fluxes [25,30,31], so neutrino detectors have unique detection possibilities. Experiments that are sensitive to MeV-scale supernova neutrinos, most prominently Borexino [32], KamLAND [33], and Super-Kamiokande (SK) [34,35], can thus be used as DM detectors as well.
This article is organized as follows: in Sec. II we provide an introduction to the singlet majoron model, its couplings and decay modes. In particular, we introduce a compact parametrization for the one-loop induced majoron couplings to charged fermions that is invaluable to study majorons. In Sec. III we discuss the signatures of majoron DM, split into neutrino signatures (Sec. III 1) and visible decay modes (Sec. III 2). Low-energy constraints from e.g. LFV that are also relevant if the majoron is not DM are presented in Sec. IV. Finally, we conclude in Sec. V. Appendix A is devoted to a discussion of neutrino flavor ratios after propagation, highlighting the differences between neutrino production from electroweak and majoron interactions.

II. SINGLET MAJORON MODEL
We know from neutrino-oscillation experiments that at least two neutrinos are massive, with sub-eV mass splitting. If neutrinos are Majorana particles, this implies that lepton number U (1) L (or U (1) B−L ) is broken by two units, and if this breaking is spontaneous we expect a Goldstone boson, the majoron J [3,4]. We will restrict ourselves to the singlet majoron model, where an SM-singlet complex scalar σ = (f + σ 0 + iJ)/ √ 2 with L(σ) = −2 couples to three right-handed neutrinos N R , with the lepton (scalar) doublet L (H) and the Yukawa matrices y and λ. A generalization to arbitrarily many right-handed neutrinos is straightforward and will not change the discussion [36]. Spontaneous symmetry breaking at the scale f gives rise to the right-handed Majorana mass matrix M R = f λ/ √ 2, diagonal without loss of generality. Electroweak symmetry breaking, H = (v/ √ 2, 0) T , introduces a mixing between the left and right-handed neutrinos via the Dirac mass matrix where V is the 6 × 6 mixing matrix to the states n R , which form the Majorana mass eigenstates n = n R + n c R . The relevant couplings of J, Z, and W − can be rewritten in terms of these mass eigenstates as [37] where Here, U is a unitary mixing matrix from the diagonalization of the charged-lepton mass matrix which we can assume to be the identity matrix without loss of generality. The neutrino couplings to the CP-even scalars can be found in Ref. [37] but are of no importance here. We will assume σ 0 to be very heavy, m σ 0 ∼ f v, and essentially decoupled from the SM to simplify the discussion. It could however be used as an inflaton [13,14], which has little impact on the discussion in this article.
We will further assume the majoron to be massive, i.e. a pseudo-Goldstone boson. The mass could arise because of (quantum-)gravity effects [5,38], heavy chiral fields that render B − L anomalous (which could make J an axion and identify the seesaw scale with the Peccei-Quinn scale [16][17][18][19]) or simply because of explicit breaking in the Lagrangian [9,39]. The actual mechanism is not important for our analysis, its main impact will be on the production mechanism for majoron DM and on its decay into two photons, to be discussed below. It must be mentioned, however, that the existence of U (1) B−L breaking in the Lagrangian, as required for a non-zero majoron mass, can lead to severe fine-tuning issues. Some U (1) breaking terms, such as σ 3 or even some Planck-scale suppressed operators, need to be heavily suppressed in order to keep the majoron mass small [5]. Similar issues arise in axion models [40], but can often be solved by means of additional particles and symmetries.
The limit of interest in this article is the seesaw [1] relation m D M R in Eq. (2), which leads to light neutrino masses of order m 2 D /M R , automatically suppressed with respect to the electroweak scale. This allows for a block-diagonalization and expansion in the small ratio where and describes the mixing between light and heavy neutrinos [41]. Since the mixing angles of the Pontecorvo-Maki-Nakagawa-Sakata matrix U and the mass splittings ∆m 2 21 and |∆m 2 32 | are known [42], the free parameters in the seesaw limit are e.g. m 1 , d h , R, f , m J , and the three CP-violating phases in U . It will prove useful to distinguish three different extreme hierarchies of light neutrinos: Normal Hierarchy (NH): Assuming m 1,2,3 m J m 4,5,6 , the majoron can decay into the light neutrinos, with partial widths proportional to m 2 j due to the diagonal Jνν couplings in the seesaw limit: Neutrino oscillations give a lower bound on j m 2 j of 2.6×10 −3 eV 2 (4.9×10 −3 eV 2 ) for normal (inverted) mass ordering; Cosmology gives a conservative upper limit of 0.17 eV 2 [43], to be used below for the QD regime, although much stronger limits even below 10 −2 eV 2 are possible for certain combinations of datasets [43][44][45]. We see that a majoron can easily be long-lived enough to be DM for typical seesaw scales, assuming J → νν to be the main decay channel. At the one-loop level one obtains a coupling of J to charged fermions [3,37], which is crucial for majoron phenomenology. The Feynman diagrams are shown in Fig. 1 and give rise to the effective on-shell couplings with flavor-diagonal pseudoscalar quark couplings and more involved lepton couplings, to lowest order in the seesaw limit, where T d, The partial width for the charged-fermion modes J →f f is then given by working again in the limit of small fermion masses. A couple of remarks are in order: • All couplings g Jf1f2 are proportional to the corresponding fermion masses as required for derivative couplings of Goldstone bosons. This in turn implies that the processes J → f 1 f 2 are helicity suppressed as expected for a neutral spin-zero particle decaying into SM fermions.
• The diagonal fermion couplings g Jf f are of pure pseudoscalar nature [37].
• The off-diagonal lepton couplings g J are approximately chiral due to the hierarchy of charged-lepton masses, for m m .
• The matrix K is positive semi-definite if the lightest neutrino mass is zero and positive definite otherwise, with determinant and non-negative trace, trK ≥ 3(det K) 1/3 . All diagonal entries K are real and non-negative. Since m lightest = 0 is unstable under renormalization group evolution [46], we can take K to be strictly positive definite, which gives Schwarz inequalities on the off-diagonal entries, As a result, constraints on trK, e.g. from J →qq, constrain all entries of K, courtesy of its positivedefinite nature.
• From the definition K = m D m † D /(vf ) we can estimate a simple perturbativity condition by demanding m D /v < √ 4π (see also Ref. [47]), Typical values for K -without fine-tuned matrix cancellations, i.e. imaginary R -can on the other hand be estimated as with the Yukawa coupling λ from Eq. (1). This is of course nothing but the result one obtains for one fermion generation, as calculated in Ref. [3].
• As shown in Ref. [48], the matrix m D m † D (or K in our case) can be used to replace R and d h in the seesaw parametrization. In other words, the entire seesaw matrix from Eq. (2) can be reconstructed using low-energy neutrino parameters (d l and U ) as well as K and f . This is the parametrization of choice in this article, seeing as K describes the physical couplings of the majoron to charged fermions and furthermore fulfills a number of useful inequalities that would be tedious to translate to e.g. R. It is quite remarkable that the seemingly lost high-energy seesaw parameters encoded in m D m † D become available in the form of majoron couplings, allowing in principle to reconstruct the seesaw mechanism with low-energy data.
Due to the proportionality Γ(J →f f ) ∝ m 2 f , the dominant decay channel of J is typically into the heaviest kinematically available fermion, but there are some notable loopholes: 1) the decay rates into charged fermions all scale with K 2 and can be made small compared to J → νν in the limit λ √ 2d h /f 1; 2) the diagonal lepton couplings J are proportional to trK − 2K , which could be highly suppressed for up to two leptons despite K being large [37,49]. For example, the pattern K ee = K µµ K τ τ 0 turns off the majoron couplings to ee and µµ.
Since we are interested in majoron masses in the MeV-GeV range, the decays J →ūu,dd,ss,cc should be replaced by appropriate decays into hadrons, which in particular moves the kinematic threshold from m J 2m u to m π , with first allowed channel J → πγγ [50], albeit heavily suppressed. Note that J decays into pairs of pseudoscalar mesons are forbidden by CP, so the next threshold is 3m π [51]. Seeing as not even the hadronic decay modes of a CP-even Higgs-like scalar with mass between 0.1-10 GeV have been agreed-on in the literature (see e.g. Ref. [52]), we will not attempt to derive the J → hadron decay rates here, but leave them for future work. Estimates for a pseudoscalar's decay into three mesons can be found in Ref. [53], assuming Higgslike couplings. In the majoron model we have instead a Higgs-like coupling with additional sign-flip for up-and down quarks, just like in two-Higgs-doublet models of type I and X. The only hadronic decay used in the following is J →bb, which can be calculated reliably and will provide the best constraints on K for m J > ∼ 10 GeV. Let us continue our discussion of majoron decay modes. Still at the one-loop level one has virtual internal Bremsstrahlung, J →f f γ, simply by attaching photons to the diagrams in Fig. 1. For quarks this merely gives the well-known final-state radiation spectrum, but the additional diagram with a W boson gives a more interesting result for leptons. The extra photon removes the helicity suppression of the amplitude and leads to a photon spectrum similar to the s-wave Majorana DM annihilation intof f γ [54], with characteristic shape for sizable photon energy In our case, the helicity suppression of the amplitude A ∝ m is however replaced by an additional heavy-neutrino propagator, A ∝ e m 3 J /d 2 h , so the rate is of higher order in the seesaw expansion and hence strongly suppressed. Bremsstrahlung will therefore not give testable signatures and will not be discussed further.
Lastly, let us mention the possible decay mode J → γγ, which could be the prime discovery channel and has been discussed extensively in the literature for other models. For a massless majoron, the coupling to photons vanishes because the global U (1) B−L symmetry is anomaly free [9,37]. The coupling for a pseudo-Goldstone boson then depends on the UV completion of the theory, i.e. the details of how m J is generated (and whether the singlet has some admixture of a triplet majoron [4]). In the absence of U (1) B−L -anomaly-inducing heavy fermions, our singlet-majoron coupling to photons will be generated first at two loops. One contribution comes from the majoron mixing with the longitudinal component of the Z boson, which then decays into two photons, see Fig. 2. Notice that only fermion loops contribute to this piece of the amplitude, because similar diagrams with the W boson and its Faddeev-Popov ghosts cancel each other [55]. The additional diagrams that arise from closing the leptonic lines in the W -boson loop of Fig. 1 b) are much more complicated to calculate, but we expect them to be further suppressed by the W mass or even the heavy neutrino mass, so we will neglect them for now. Notice that such a separation of the diagrams is gauge invariant, as the corresponding amplitudes satisfy the Ward identities separately. Note also that the Z-boson contribution depends on different parameters than the W -boson part (e.g. quark masses), so it is not possible for the neglected diagrams to cancel the entire amplitude; a partial destructive interference could, of course, be possible. Focusing only on the gauge-invariant part of the amplitude induced by J-Z mixing, i.e. Fig. 2, the two-loop rate takes the simple form with the color factor N q c = 3, N c = 1 and the loop function For m J m e , the fermion-mass independent contributions cancel due to anomaly freedom, leading to a rate that is dominated by the lightest fermion, In particular, the coupling Jγγ vanishes for m J = 0 as expected. Up to a prefactor, the rate of Eq. (27) is equivalent to the singlet-triplet majoron case, where the majoron-Z mixing is induced already at tree level by the vacuum expectation value of an SU (2) L triplet ∆ → v T / √ 2 [8]. The singlet-triplet-majoron rate then follows from Eq. (27) We stress once more that the above diphoton rate was obtained by considering only a (gauge-invariant) subset of two-loop diagrams. While we expect the remaining diagrams to be suppressed by m W or d h or even cancel completely, a full calculation is beyond the scope of this article. Furthermore, the rate can be modified by the details of the scalar (admixture of triplets or CP-violating mixing with the Higgs) and fermion sector (B −L anomalous fermions that create a Jγγ coupling at one-loop). The reader should therefore be careful when interpreting the diphoton rate used here.

III. DARK MATTER
Possible production mechanisms for majoron(-like) DM have been extensively discussed in Ref. [9], assuming a restricted set of couplings in order to obtain predictions. For example, taking to be the only explicit U (1) breaking term in the scalar potential and neglecting the U (1)-invariant portal |σ| 2 H † H, the relic density Ω J of J is completely fixed for a given mass m J (assuming small Yukawa couplings λ to the heavy neutrinos). For sufficiently large λ h = m 2 J /v 2 , a thermal population of majorons is produced in the Early Universe from annihilations and the (inverse) decays of the Higgs boson; after the Higgs disappears from the thermal plasma, the DM density eventually freezes out. The required value for λ h in this scenario is typically incompatible with constraints from direct detection or h → invisible, at least in the mass range of interest here [56]. Another possible situation is to assume that the number of DM particles was negligible with respect to those of the SM after reheating. If the portal interaction coupling λ h has small values, the population of majorons never reaches thermal equilibrium; for temperatures much smaller than the Higgs mass -when the majoron decouples from the SM plasma -its abundance approaches a constant value. This leads to [57] where g s * and g ρ * are the number of degrees of freedom contributing to the entropy and energy density when the majoron decouples. This is the freeze-in mechanism, which obviously requires m J < m h /2 and a very small decay rate of the Higgs boson into majorons (automatically satisfying LHC constraints on h → invisible). From the observed DM density, and taking g s * ∼ g ρ * ∼ 100, we obtain m J 2.7 MeV for the λ h = m 2 J /v 2 case described above [9].
In a more general case, one can consider separate U (1) breaking terms for the majoron mass and the Higgs portal, disentangling relic density and DM mass. For the freeze-out production mechanism, this is just the singlet DM scenario, heavily constrained and only viable around the Higgs resonance [56]. For the production via freezein, Eq. (31) leads to Freeze-in is thus a viable mechanism to produce majoron DM in the MeV and GeV range. Other production mechanisms exist, see e.g. Ref. [9] and references therein. Nevertheless, from now on we will remain agnostic about how DM was produced in the Early Universe and only assume that (cold) majorons constitute all the DM and that its mass lies below the electroweak scale. In any case, the specific indirect detection signatures discussed below do not depend on the details of the majoron mass generation or its production mechanism.

Neutrino signatures
The only tree-level decay mode of the singlet majoron J is into neutrinos, Eq. (13), completely specified in terms of neutrino masses and U (1) breaking scale f . An interesting side effect of the majoron coupling to neutrino mass eigenstates is that the emitted neutrinos will not oscillate, resulting in flavor ratios that can be completely different from astrophysical sources [58]; for a detailed discussion using the density-matrix formalism, see Appendix A. The branching ratio of J decaying into ν j is proportional to m 2 j , and ν j contains a fraction |U j | 2 of flavor , so the flavor composition of the majoron-decay neutrino flux is given by α e : α µ : α τ with normalized so that α = 1. The self-conjugate Majorana nature ensures that α = α¯ . Contrary to most other neutrino fluxes, these ratios are the same at the source, where DM decays, and on Earth, so α = α S = α ⊕ , up to matter effects inside the Earth. See Fig. 3 for an illustration using a ternary plot with a scan over the 1σ and 3σ ranges of the oscillation parameters obtained in Ref. [42].
The mixing angles θ 23 π/4 θ 13 result in an almost µ-τ -symmetric mixing matrix, i.e. |U µj | |U τ j |, which ensures α µ α τ independent of the mass ordering. α e on the other hand depends strongly on the neutrino mass hierarchy, with lowest value for NH (α e sin 2 θ 13 ) and largest value for IH (α e 1/2). Using the best-fit values from Ref. [42] for the mixing angles, we obtain the following benchmark values for the flavor ratios in the hierarchical regime,   Fig. 3. These are the values we will use in the following, but most results can be rescaled The 1σ (3σ) ranges of the neutrino-oscillation parameters from Ref. [42] correspond to the green (blue) lines; lighter colors correspond to a larger lightest-neutrino mass, converging to 1 : 1 : 1 for the QD spectrum. The three stars denote the benchmark values of Eq. (34). The expected flavor ratios from realistic astrophysical processes (e.g. pion decay followed by averaged-out neutrino oscillations) fall in the red contour, taking into account the uncertainties in the mixing parameters (95% C.L.) [59].
without much effort. The NH composition with its tiny ν e fraction α e sin 2 θ 13 is particularly interesting, because there is no astrophysical mechanism that would suppress ν e to such a degree without physics beyond the SM [58]. This is illustrated in Fig. 3, where we also show the expected flavor ratios from astrophysical processes (red contour) under the assumption that the neutrino oscillations have been averaged out when the flux arrives at Earth [58,59], see Appendix A for details. As can be seen, the NH region lies outside of the typical astrophysical expectation, making flavor ratios a potential discriminatory tool for DM detection.
Seeing as the majoron itself forms cold DM in our model, the neutrino spectrum with its line-like feature will be a much better discovery tool than the flavor ratios of Fig. 3. Let us mention, however, that the monochromatic signature becomes less important as soon as we consider m J above the electroweak scale; since the coupling to neutrinos of Eq. (3) also induces a coupling to the SM Higgs of the form Jν j ν j (m j /f )(1 + h/v) 2 , the decay modes J → ννh(h) open up for m J > (2)m h , and in fact dominate over J → νν for m J > ∼ 10 TeV [62]. The neutrino spectrum from J → ννh(h) is then obviously no longer monochromatic, but the flavor ratios of the primary neutrinos illustrated in Fig. 3 continue to be valid. In addition, there will be secondary neutrinos with a different spectrum and flavor ratio from the h decay and electroweak Bremsstrahlung. A thorough discussion of these effects will be discussed elsewhere, but we expect the flavor ratios of the secondary neutrinos to fall into the red contour of Fig. 3, because they are created as flavor eigenstates (see Appendix A). Let us also mention that in models with a larger dark sector it is possible to obtain, for example, boosted majorons that decay into a continuous neutrino spectrum, for which the flavor ratios could again be more important.
As mentioned above, the spectral feature of J → νν should serve as a sufficient discriminant from the continuous background. As shown in Ref. [62], this two-body decay mode is suppressed compared to J → ννh(h) for m J > ∼ 10 TeV, which induces a continuous spectrum. We will further restrict ourselves to masses m J < 100 GeV in this analysis, in order to avoid discussing effects from e.g. J → W W, ZZ that could be induced in some UVcompletions of our model. We stress, however, that J → νν could still be an important discovery channel for majoron masses up to 10 TeV, for which IceCube becomes the ideal observatory [63,64]. The neutrino (plus antineutrino) flux per flavor from the J → νν decay in our galaxy is given by [26,65] where J = ∞ 0 ρ Halo ds is the astrophysical factor associated to the DM density ρ Halo in the Milky Way halo. For simplicity we write here the flux associated to the full sky, the general case for an angular signal is a straightforward generalization of this case. The J -factor introduces uncertainties in the determination of the flux because the precise shape of ρ Halo is unknown in the center of the Galaxy. Nevertheless, in contrast to DM annihilations for which the J factor scales quadratically with ρ Halo and thus varies by many orders of magnitude depending on the assumptions on the DM halo, the uncertainty for DM decays is of less than one order of magnitude [26] and the determination of neutrino fluxes or limits on them is more robust. Notice that here we are neglecting the neutrino flux arising from DM decays outside our Galaxy, whose spectrum is in any case red-shifted and not necessarily line-like [66,67].
For the two-body decay J → νν we have dN/dE ν = 2δ(E ν − m J /2), which is smeared out by the velocity distribution and detector resolution. Low-threshold neutrino detectors such as Borexino [32], KamLAND [33], and SK [34,35] give limits on the (monochromatic) flux ofν e from searches for the diffuse supernova neutrino background. Due to the large cross section and tagging possibilities, the detection channel of choice here is inverse beta decayν e p → ne + , which has a kinematic threshold of E ν > 1.8 MeV. This makes it difficult to obtain limits for m J < ∼ 4 MeV, seeing as the background from reactor neutrinos also increases dramatically for such low energies. For 5 MeV < ∼ m J < O(100) MeV on the other hand, searches for supernovaν e neutrinos give useful constraints on DM-induced neutrino fluxes, as can be seen in Fig. 4. Note that in our notation this is a limit on the flux Φν e = 1 2 Φ e , because only half of our electron neutrinos are antineutrinos. A near-future improvement of these limits is realistic, especially with the  [60]. Adopted limits from supernovā νe searches come from Borexino [32] (green), KamLAND [33] (red), SK [34,35] (blue), and reinterpreted SK data (orange) [26]. The black lines for mJ > 0.1 GeV come from a comparison with atmospheric νµ spectra [26], while pink shows the preliminary limit from a designated DM search using angular-anisotropy SK data [61].
proposed Gadolinium-extension of SK [68], which should reduce background and potentially reach the diffuse supernova regime. Even ton-scale liquid-xenon detectors build for the direct detection of DM, such as XENONnT, LZ or DARWIN, could be sensitive to O(10 MeV) neutrino lines and might give useful information about the flavor ratios [69]. In any case, dedicated DM searches by the experimental collaborations are desirable, especially considering the apparent gap of official limits between m J = 60 MeV and GeV. Above GeV, we have preliminary SK limits on DM decay into muon neutrinos [61]. In the gap 60 MeV < m J < GeV, we adapt the limits of Ref. [26], based on a reinterpretation of SKν e data as well as a comparison to the well-understood atmospheric muon neutrino flux (see also Ref. [67]). Here, we strongly urge the SK collaboration to check for neutrino lines, both in electron and muon neutrinos. Hyper-K is expected to further improve on the higher-energy region.
Depending on the neutrino mass hierarchy, these flux limits can be translated into a lower bound on the U (1) B−L breaking scale f , see Fig. 4. The latter is naturally strongest for QD, seeing as Γ(J → νν) ∝ j m 2 j /f 2 scales with the neutrino masses. In contrast, the weakest bounds arise for NH, which is quite obvious for limits that come from bounds on the total lifetime or from the α e sin 2 θ 13 suppressed ν e flux; surprisingly, limits from Φ µ lead to roughly the same bounds on f for NH and IH, because of the accidental numerical relation For the sake of clarity, it is therefore sufficient to discuss the limits for the regimes QD and NH in Fig. 4, as those associated to IH happen to fall in between.
Less direct limits on the J → νν decay come from cosmology. The most conservative bound is surely to demand J to have a lifetime that exceeds the age of our Universe, τ 4 × 10 17 s. Better limits can be obtained by studying the effect that the decay of a nonrelativistic DM particle into relativistic daughter particles has on e.g. the matter power spectrum. A recent analysis provides a 95% C.L. constraint of order τ > 5 × 10 18 s [70]. Future measurements of the cosmic microwave background (CMB), e.g. by CORE, could improve the bound on τ by a factor of 2 [71]. This is currently the only constraint on the J → νν mode for majoron masses below 4 MeV and will be hard to improve with line searches due to the huge neutrino background below 10 MeV from e.g. reactor neutrinos [29].
The limits on f from Fig. 4 can be translated into upper bounds on |K αβ | with the help of the perturbativity constraint of Eq. (24). For m J = 1 MeV (100 GeV) this implies |K| < 5 × 10 −6 (3 × 10 −11 ) for NH and about an order of magnitude stronger for QD. These limits are much weaker then the direct constraints from J →f f derived below (Fig. 5), but are valid even if the J decay is kinematically forbidden.  [72] and indirect DM searches with AMS-02 [65,73]; γ-ray telescope limits on J → γγ and J →bb refer to INTEGRAL [74] for mJ < 7 MeV, to COMPTEL/EGRET [75] for 7 MeV ≤ mJ ≤ 400 MeV, and to Fermi-LAT [76,77] for mJ > 400 MeV. For indirect DM searches, we only show the most constraining limits in a given channel. We remind the reader that K is expected to have an order of magnitude of 10 −13 λ, where λ is the Yukawa in Eq. (1).

Signatures from visible decay channels
Having identified MeV < ∼ m J < ∼ 100 GeV as the region of interest where majoron DM can lead to a particularly clean observable flux of monochromatic neutrinos, let us discuss the constraints from the visible decay channels, i.e. J →f f at one loop and J → γγ at two loop. There are stringent constraints on DM decays into charged fermions from a wide range of indirect searches, see e.g. Ref. [78] for a review. In our model, the majoron decay modes into charged fermions all depend on the matrix K = m D m † D /(vf ) introduced in Eq. (18). A crucial observation here is that K does not depend on the lowenergy neutrino parameters, but is a completely free parameter matrix in the seesaw limit, up to the inequalities given below Eq. (18). Typical values can be estimated as K ∼ d h d l f v 2 × 10 −13 λ, but it is entirely possible to have values orders of magnitude larger or smaller. While the J → νν modes discussed above gave a direct limit on the seesaw scale f , the charged-fermion decay modes will give limits on the remaining parameters of our model, which are encoded in the elements of K. The decays J → νν and J →f f thus provide completely orthogonal information about majoron DM.
Majoron decays into charged leptons are in particular constrained by the AMS-02 measurements of the positron flux in cosmic rays [65]. The corresponding 95% C.L. upper bounds on the K matrix elements are shown in Fig. 5. For masses above a few GeVs, other limits on the leptonic decay channels arise from the diffuse-gamma-ray obser-vations of the sky [75,77,[79][80][81][82][83], but these are typically less stringent than those of positrons for DM masses below 100 GeV. In addition, for m J > ∼ 10 GeV, the majoron decays dominantly into bottom quarks, which subsequently decay and fragment producing antiprotons. The AMS-02 experiment has also measured the corresponding flux [84], which, within astrophysical uncertainties, can be interpreted as originating from only cosmic ray collisions with the interstellar material [73]. Slightly stronger bounds can be obtained with Fermi [77]. This allows to set a strong upper bound on the decay rate into bottom quarks, trK < ∼ 10 −22 at 95% C.L. for m J > 10 GeV, as shown in Fig. 5.
The strongest of the indirect detection bounds is the one on trK by J →bb. As a matter of fact, this bound also applies to all entries of K due to the inequality of Eq. (23). Thus, majorons with masses greater than ∼ 10 GeV are severely constrained, because such a small K would require tiny Yukawa couplings λ ∼ 10 −9 . We expect constraints on trK from the hadronic decay modes even below 10 GeV, but the corresponding decays into mesons are difficult to calculate reliably. Notice that because of these constraints, it is hopeless to observe Majoron DM in direct detection experiments looking for nuclear recoils.
Below few GeVs, indirect detection bounds become very weak compared to CMB bounds. If DM decays into photons or charged particles during the time between recombination and reionization, when the Universe was transparent and no large-scale structures were formed, it injects energy into the photon-baryon fluid and potentially modifies the anisotropies of the CMB and its blackbody shape. Consequently, the precise measurements of Planck [85] set stringent constraints on majoron decays into charged fermions. We calculate the corresponding constraints 1 on K following Ref. [72], and show them in Fig. 5. These bounds are very important for two reasons. On the one hand, they constrain majoron decays at the MeV scale, where J → e + e − and J → µ + µ − are the dominant decay channels. On the other hand, they do not suffer from astrophysical uncertainties such as those associated to halo DM densities or cosmic-ray propagation parameters.
Finally, let us discuss constraints from J → γγ, arguably the most popular decay channel for majorons [8,10,11]. Using our estimate for this two-loop decay of Eq. (27), we can translate γ-line limits from IN-TEGRAL [74], COMPTEL/EGRET [75], and Fermi-LAT [76,77] into upper bounds on trK (Fig. 5). These γ-ray telescope limits are stronger than the corresponding CMB limits on J → γγ [72], so we will not show them here. Due to the suppression by α 2 and an additional loop compared to J →f f , the limits from J → γγ are for the most part weaker than those from charged fermions. Nevertheless, the diphoton decays probe trK, which in turn limits all entries of K via Eq. (23), whereas the J → decays only probe specific linear combination of K elements. This makes the J → γγ (and J →bb) constraints particularly interesting. Future prospects for this channel are good, with considerable current effort to improve limits in the MeV gap between 7 MeV < ∼ m J < ∼ 400 MeV, for example by AdEPT [86] and e-ASTROGAM [87]. An improvement by several orders of magnitude seems feasible, which could open the door to a double-line observation in the MeV range, both in neutrinos and γ-rays. For m J < MeV, the diphoton decay is the only feasible DM detection channel, seeing as sub-MeV neutrinos are extremely difficult to detect, especially when it comes to their spectral shape and flavor.
In summary, the constraints on majoron DM from its visible decay channels provide information on the model that is complementary to the main decay mode J → νν. In the region of interest for neutrino lines, MeV < ∼ m J < ∼ 100 GeV, the constraints on the elements of K range from 10 −13 to 10 −23 , which translates into typical values for the Yukawa coupling λ of 1 to 10 −10 via Eq. (25). This should not be taken too literally in the three-generation framework, but can give some idea about the level of tuning necessary to evade the bounds. In particular, the region m J > ∼ 10 GeV could be regarded as less motivated, which is however a highly subjective statement.

IV. OTHER CONSTRAINTS
For m J > m f1 + m f2 , the best constraints on the majoron couplings g Jf1f2 come from the decay J → f 1 f 2 or J → γγ, as we have seen in Fig. 5. Let us briefly discuss limits from the production of J, e.g. from f 1 → f 2 J, which gives limits on g Jf1f2 for m J < m f1 − m f2 . For m J > MeV, all these constraints turn out to be weaker than the perturbativity bounds of Eq. (24) in connection with the limits on f from Fig. 4, which imply that |K| can be at most 5 × 10 −6 for m J MeV. Even stronger bounds apply when considering the limits from J → γγ (Fig. 5). We nevertheless list the direct constraints below for completeness, and stress that they can be relevant for m J < MeV or if J makes up only a subcomponent of DM.
The off-diagonal majoron couplings are directly constrained by the lepton flavor violating (LFV) decays → J [37,88,89], with strongest bound in the muon sector, Br(µ → eJ) < 2.6 × 10 −6 [90], and Br(τ → J) < O(10 −3 ) [91]. The strong µ → eJ bound of Ref. [90] rests on the assumption of isotropic electron emission; in our case, however, the emission is maximally anisotropic, see Eq. (21), with dominant emission of the left-handed electron in the direction opposite to the muon polarization. This also happens to be the region where the background from µ → eνν is largest, diminishing the limit by an order of magnitude [92] to |K µe | < ∼ 1 × 10 −5 for m J m µ . An almost identical limit has been obtained long ago by considering µ → eJγ with a massless J, which does not depend on the chirality properties of the Jµe coupling, but is of course further suppressed by α and phase space [93]. We checked explicitly that the rate for µ → eJγ in our model is well described by the effective off-diagonal coupling of Eq. (21) followed by Bremsstrahlung, leading to the same differential distributions given in Refs. [89,93]. 2 Since the Bremsstrahlung rate formally diverges for small photon energies and small electron-photon opening angle, the number of events crucially depends on the detector resolution. It would be interesting to see how current and future experiments such as MEG and Mu3e can improve on these 30-year-old limits with their modern detectors [89], but this will be discussed elsewhere.
For m , m J m , the partial widths are simply which then translate to the bounds neglecting the dependence on the majoron mass for simplicity. Perturbativity plus J → νν limits give stronger limits, unless m J < MeV; J → γγ even requires m J < 10 keV for LFV to be observable, at least if J makes up 100% of DM. Since such low-mass DM is typically not cold, a dedicated analysis is necessary to evaluate its validity. Additional LFV in the form of → γ arises from the right-handed neutrinos, which is independent of the majoron or breaking scale, with the strongest bound coming from Br(µ → eγ) < 4.2 × 10 −13 [94]. In the seesaw limit, m 1,2,3 m W m 4,5,6 , the partial widths take the form [36,95] which has a different matrix structure than K, making it difficult to directly compare limits. In principle one can calculate the above for a given d l , U , K and f [48], but the expression will be far from illuminating. Large rates typically require some fine-tuning, e.g. large Im(R) in the Casas-Ibarra parametrization, or a symmetry-motivated structure in m D [96]. Let us focus on the case of degenerate heavy neutrinos, d h = M × 1, for which the above is proportional to |K | 2 , allowing us to directly compare the two LFV rates, The ratio is heavily suppressed for M ∼ f m , making the majoron final state the prime LFV channel despite its more difficult signature; the photon rate can dominate for small Yukawa coupling, λ = √ 2M/f 1, implying nottoo-heavy right-handed neutrinos. Both channels should hence be searched for experimentally.
The diagonal majoron couplings, i.e. the diagonal K entries, are constrained via the J coupling to electrons and quarks. At low energies, we typically require the couplings to nucleons N = (p, n) T instead of quarks, which can be estimated naively as JN iγ 5 σ 3 N m N trK/(16π 2 v). The coupling to quarks and nucleons is of particular interest, because it depends on trK, which automatically limits all entries in K due to Eq. (23), even the LFV couplings. Limits on (light) pseudoscalars can be readily found in the literature, usually assuming an effective coupling to fermions that is then used to calculate scattering processes etc.; this will be at best an approximately accurate procedure in our model, because our effective Jf f couplings from Eq. (14) are by construction only valid for on-shell particles. As such, scattering processes -which necessarily involve off-shell particles -would have to be calculated from scratch using the loop diagrams to obtain the correct dependence of the cross sections on our parameters.
A full calculation of all the required scattering rates being beyond the scope of this work, let us assume that the Jf f couplings provide a reasonable estimate for majoron scattering. For m J < 10 keV, the best limits then come from astrophysics and imply |K ee − K µµ − K τ τ | < 2 × 10 −5 , trK < 10 −5 , (40) from the electron [97] and nucleon coupling [98], respectively. For m J up to 100 keV one has (slightly weaker) direct-detection bounds on g P Jee from EDELWEISS [99], XENON [100], XMASS [101], and MAJORANA [102], assuming J to be DM; this gives |K ee − K µµ − K τ τ | < ∼ 10 −4 [101] for m J = 100 keV, roughly ten orders of magnitude weaker than the bound at m J = O(1) MeV (Fig. 5). The couplings to quarks are much less constrained for m J > 10 keV; since there are no flavorchanging processes in the quark sector mediated by the majoron at the one-loop level, quark-flavor constraints are suppressed. Going to the next loop level we can estimate constraints from K → πJ etc. following Ref. [53], which give constraints trK < ∼ 2×10 −2 for m J < 100 MeV, much weaker for larger m J . In the region of interest in this article, MeV ≤ m J ≤ 100 GeV, majoron production gives weaker limits on K than perturbativity in combination with the neutrino limits on f , and much weaker than the J → γγ bounds.
Lastly, let us mention another signature of our model: neutrinoless double beta decay 0νββ [103]. In the seesaw limit, the amplitude for this ∆L = 2 process is dominated by light-neutrino exchange, proportional to (U d l U T ) ee = 3 j=1 U 2 ej m j . This is in particular sensitive to the Majorana CP phases in U , which cannot be measured via neutrino oscillations. Current experiments probe the QD regime, with limits of order |(U d l U T ) ee | < 0.2 eV [104]. Future experiments are expected to ultimately reach the IH regime, while NH leads to discouragingly small rates. The observation of 0νββ would be an incredible discovery and prove beyond doubt that neutrinos are Majorana particles, leading further credence to the seesaw mechanism. This would of course be good news for our majoron model at hand, as it would in particular fix the rather strong dependence of e.g. J → νν on the neutrino hierarchy. It should be mentioned, however, that our (sub-MeV) majoron DM gives completely negligible rates for the associated 0νββJ process (A, Z) → (A, Z + 2) + 2e − + J [103], seeing as the majoron couplings to neutrinos m ν /f are minuscule. The discovery of such a mode would therefore strongly hint at a more complicated majoron realization. Due to the small Jνν coupling, supernova constraints are also easily evaded [105].

V. CONCLUSION
In this work, we have revisited the singlet majoron model, in which lepton number is a nearly exact symmetry that is spontaneously broken at the seesaw scale. The corresponding pseudo-Goldstone boson, the majoron, is stable on cosmological scales due to its highly suppressed couplings and can act as DM. One of the most remarkable features of this model is the prediction of monochromatic neutrinos arising from DM decays, practically testable at energies between MeV (e.g. Borexino) and 100 GeV (e.g. Super-K), potentially up to 10 TeV (IceCube). We urge the experimental collaborations to perform designated searches for such low-energy DM-induced neutrino lines. Since the majoron couples directly to the neutrino mass eigenstates, the decay neutrinos do not oscillate and have flavor ratios on Earth that depend strongly on the neutrino mass hierarchy, see Figs. 3 and 4. In particular, the electron-neutrino flux is suppressed compared to the other flavors for the normal mass hierarchy.
Other constraints on the model arise from the majoron couplings to charged fermions, induced at the one-loop level, and the decay into two photons, induced by twoloop diagrams. We have provided a convenient and compact three-generation parametrization of these couplings in terms of the matrix m D m † D , which contains precisely those seesaw parameters that are usually unobservable at low energies. A measurement of the majoron couplings could then in principle complete our knowledge of the seesaw mechanism. In the DM context, majoron decays into charged fermions and diphotons are constrained by CMB observations and indirect DM searches, which put strong limits on m D m † D , especially for m J > 10 GeV, as illustrated in Fig. 5. Our parametrization also allows us to study constraints from lepton flavor violation; the rates for anisotropic → J turn out to be small for m J > ∼ MeV if J makes up all of DM, but → γ can be observable for not-too-heavy right-handed neutrinos. m 2 i 2E )δ ij . In fact, the density matrix describing the flux of neutrinos after a distance L at Earth reads ρ ⊕ = e −iHL ρ S e iHL , or more precisely, For a sufficiently large oscillation length L, neutrino oscillations average out and exp{−i The flavor composition at Earth, given by the diagonal elements of the density matrix, is thus By varying the oscillation angles within the ranges allowed by neutrino experiments and assuming an arbitrary composition of flavors at the source, we obtain the red contour of Fig. 3. The situation is different for neutrinos arising from majoron decay. In this case, the branching ratios associated to J → ν i ν j are proportional to m 2 j δ ij , at least in the lowest seesaw order we consider. Accordingly, the density matrix at the source is diagonal in the mass basis and commutes with the Hamiltonian. As a result, ρ ⊕ = e −iHL ρ S e iHL = ρ S and therefore α ⊕ = α S .