Neutrinoless double-beta decay with massive scalar emission

Searches for neutrino-less double-beta decay ($0\nu2\beta$) place an important constraint on models where light fields beyond the Standard Model participate in the neutrino mass mechanism. While $0\nu2\beta$ experimental collaborations often consider various massless majoron models, including various forms of majoron couplings and multi-majoron final-state processes, none of these searches considered the scenario where the"majoron"$\phi$ is not massless, $m_\phi\sim$~MeV, of the same order as the $Q$-value of the $0\nu2\beta$ reaction. We consider this parameter region and estimate $0\nu2\beta\phi$ constraints for $m_\phi$ of order MeV. The constraints are affected not only by kinematical phase space suppression but also by a change in the signal to background ratio characterizing the search. As a result, $0\nu2\beta\phi$ constraints for $m_\phi>0$ diminish significantly below the reaction threshold. This has phenomenological implications, which we illustrate focusing on high-energy neutrino telescopes. Our results motivate a dedicated analysis by $0\nu2\beta$ collaborations, analogous to the dedicated analyses targeting massless majoron models.


I. INTRODUCTION
Neutrinoless double beta (0ν2β) decay [1][2][3][4], is a lepton number violating process. It is sensitive to the neutrino mass parameter where m νi (i = 1, 2, 3) are the neutrino masses and U is the lepton mixing matrix [5]. While the renormalizable Standard Model (SM) has lepton number as an accidental symmetry and, consequently, predicts that the neutrinos are massless, adding dimension-five terms [6] 12 , where H is the Higgs doublet field and L α (α = e, µ, τ ) are the lepton doublet fields, leads to neutrino masses, with v = 246 GeV. We do not know the beyond-SM origin of the dimension-five terms in Eq. (3). It is possible that additional light particles accompany the neutrino mass mechanism and interact with SM fields in various ways. If there exists a light gauge-singlet scalar φ, then the dimension-six terms are possible. The dimension-six terms lead to Yukawa couplings of φ to neutrinos, If the mass of the φ particle is less than the Q-value of the (A, Z) → (A, Z + 2) transition, m φ < Q, then the G ee coupling leads to a decay where 0ν2β is accompanied by on-shell φ emission (0ν2βφ), A well known framework that leads to Eq. (5) and to the decay mode 0ν2βφ is that of majoron models [8][9][10][11][12][13], where φ is the Goldstone boson related to the spontaneous breaking of the lepton number symmetry. Many variants of the majoron model have been studied in the literature. In the simplest realisations, the seesaw scale Λ appearing in Eq. (3) is promoted to a dynamical field, and the phase of this field is associated with φ. In such models, (i) the φ particle is massless, and (ii) the terms of Eq. (3) and Eq. (5) are related, leading to G = m ν /Λ. For high-scale seesaw models, with Z = O(1), the seesaw scale is Λ ∼ 10 14 GeV, leading to G ∼ 10 −24 . As we review in Sec. II, such tiny coupling is some 20 orders of magnitude below the reach of 0ν2βφ searches.
In other scenarios, like the inverse-seesaw models of Ref. [14][15][16] (see Ref. [17] for a review), neutrino masses arise from effective dimension six terms. Namely, instead of 1/Λ in Eq. (3) we have µ/Λ 2 , where a technically natural hierarchy µ Λ is responsible, at least in part, for the smallness of the neutrino mass. In such case, a light scalar field could arise if we promote the inverse-seesaw parameter µ to a field φ with µ = φ . In this case, G = m ν /µ and if µ is small enough, 0ν2βφ could be observable. If lepton number is broken spontaneously by φ , then the φ particle is still massless.
Global symmetries, however, are not expected to be exact. If lepton number is broken not only spontaneously but also explicitly, by some small parameter, then φ could be light but not massless [18,19]. In addition, if the explicit lepton number violation (LNV) dominates the neutrino mass, then also the relation between G and m ν is modified. Yet another framework that can accommodate this situation is if neutrinos are Dirac particles, in which case lepton number (more precisely some non-anomalous symmetry group containing it, e.g. B − L) may be exact; see [20] for a recent study. While 0ν2β experimental collaborations often consider various massless majoron models, such as different forms of the majoronneutrino couplings and multi-majoron final-state processes, none of these searches considered the scenario of a massive majoron, m φ ∼ MeV, of the same order as the Q-value of the 0ν2β reaction. In this paper we consider this parameter region 3 and estimate 0ν2βφ constraints for the case of m φ of order MeV. As we show, the constraints are affected not only by kinematical phase space suppression near m φ ∼ Q, but also by a change in the signal to background ratio characterising the search. As a result, 0ν2βφ constraints for m φ > 0 diminish significantly below the reaction threshold. Our results motivate a dedicated analysis by 0ν2β collaborations, analogous to the dedicated analyses targeting different massless majoron models.
The constraint on massive φ emission in 0ν2βφ has phenomenological implications, which we illustrate focusing on highenergy neutrino telescopes. Light scalar fields coupled to neutrinos were considered as mediators of anomalous neutrino self-interactions in many other works. Refs. [22,23] studied the effect of light scalar exchange on the energy spectrum of ∼10 MeV neutrinos from core-collapse supernovae (see also [24] where supernovae neutrinos scatter on dark matter). Vector boson or massless majoron exchange were considered in [25][26][27]. Refs. [28][29][30][31] discussed the relation of anomalous neutrino interactions to low-scale neutrino mass generation, focusing on spontaneously broken global and gauged lepton number. Ref. [32] extended the discussion to the technically natural possibility of small explicit LNV, and made a connection to phenomenology at high-energy neutrino telescopes. Recently, Ref. [33] considered light scalar exchange in coherent neutrino-nucleus scattering.
Before we turn into concrete calculations, let us emphasize that while 0ν2β is a LNV process, 0ν2βφ could be lepton number conserving (LNC). It could be therefore that the latter is strongly enhanced compared to the former. Explicitly, for where we used the fact that m ee ∼ < 0.1 eV [34], with Q ∼MeV. As we review in the next section, 0ν2βφ searches have reached a limit |G ee | 10 −5 (for massless φ). The reason that Γ 0νββφ Γ 0νββ , as shown by Eq. (9), is consistent with these limits, is related to the difference in the visible electron energy spectrum between these decay modes, which reduces the signal to background ratio for 0ν2βφ compared with 0ν2β. In what follows we will see an anaolgous effect deteriorating the sensitivity to m φ > 0 compared to the m φ = 0 case. Finally, note that when m φ > Q, the on-shell process 0ν2βφ is kinematically blocked, but the off-shell process 0ν2β(φ * → 2ν), where a virtual φ is emitted and decays to two neutrinos, is always allowed. However, compared to the on-shell process 0ν2βφ (when allowed), the off-shell φ process is strongly suppressed by a factor ∼ 2|G| 2 2 . In addition, the spectral shape with respect to the outgoing electron energy is similar to that of the standard background process 2ν2β. These features are explained in App. A. As a result of these features, the virtual φ decay mode 0ν2β(φ * → 2ν) cannot be constrained with current experiments, and we limit our attention to on-shell 0ν2βφ.

II. NEUTRINOLESS DOUBLE-BETA DECAY WITH MASSIVE SCALAR EMISSION
From the list of 0ν2β experiments surveyed in [34], NEMO-3 [35] using 100 Mo has the highest Q-value, Q ≈ 3.03 MeV. Recent work by NEMO-3 allowed them to surpass this record using 150 Nd [36], with Q ≈ 3.37 MeV, albeit with lower exposure. Thus in principle 100 Mo and 150 Nd experiments probe the highest scalar mass. The strongest constraint on the massless majoron case is from KamLAND-Zen [37] using 136 Xe, which has a somewhat lower value of Q ≈ 2.5 MeV. For these reasons -sensitivity to the highest m φ , and current best sensitivity to massless φ -we focus on 100 Mo, 150 Nd, and 136 Xe in our numerical analysis below. It is straightforward to extend our analysis to other isotopes of common use, like 76 Ge [38,39], 82 Se [40], and 130 Te [41,42]. These isotopes yield comparable, although (currently) somewhat weaker constraints.
Refs. [35], [36], and [37] provide 90%CL bounds on the massless majoron scenario, equivalent to |G ee | < (1. respectively. These bounds are stronger than other constraints in the literature such as those arising from light meson decay (see, e.g. [43][44][45]) and from cosmological and astrophysical considerations [32]. While the massless majoron bounds [46][47][48][49][50] coincide with our model for m φ Q, to our knowledge a study of the kinematical region m φ ∼ Q ∼ MeV has not been done and we consider this region in what follows.
Searches for 0ν2βφ constrain the half-life time for the decay, T 1 2 , which in our model can be approximately decomposed as [1,46,51] (see, e.g. [52] for a recent account) where G ee is defined in Eq. (7), M is a dimensionless nuclear matrix element (NME), and G(m φ ) is a kinematical phase space factor, conventionally expressed in units of yr −1 . The main effect of massive φ emission is to modify the phase space factor [13], The outgoing electrons Coulomb factor, encoded in a(E 1 , E 2 ), is given in Refs. [1,46]. For the range of electron momenta of interest, . For simplicity, for most of the numerical results in this work we use the non-relativistic approximation, We have checked that our results are not affected significantly when using the relativistic expressions for the electron wave function.
The 0ν2βφ constraint on G ee deteriorates as the scalar mass m φ approaches the kinematical limit for the decay. The constraints are affected in two ways: 1. The phase space factor diminishes close to the kinematical limit, G(m φ )/G(0) → 0 as m φ → Q.
2. As m φ is increased, the visible final state electrons kinetic energy is pushed to lower values, because the usual massless majoron phase space factor This brings the distribution of T 2β to overlap more with the T 2β distribution of the irreducible Standard Model 2ν2β background, leading to smaller signal to background ratio in the experimentally relevant energy range.  Proceeding further, we consider first the simple phase space suppression encoded in G(m φ )/G(0). In the Primakoff-Rosen approximation (PRA) [53], we can calculate this ratio analytically. We find: Taking the limitm φ → 0, we recover the phase space factor for massless majoron (see, e.g. [9]). The phase space suppression factor is shown in Fig. 2 for 100 Mo (blue) and 136 Xe (red). The analytical PRA result Eq. (14) is shown in solid lines. Numerical computation using the full relativistic electron wave function is shown by dots. We see that phase space suppression is appreciable already for m φ ∼ 1 MeV. In the case of 136 Xe ( 100 Mo), for m φ ≈ 2 MeV (2.5 MeV), the decay width drops to ∼1% of its value for the massless majoron case. This means that the limit on G ee becomes weaker than the massless majoron limit by a factor of at least ∼10 at that point, even ignoring the signal to background ratio deterioration effect that we discuss later on.
We next consider the varying signal to background ratio at varying m φ . Accounting for this effect properly is difficult outside of the experimental collaboration, due, among other factors, to different experiment-dependent sources of background, which make the derived limit sensitive to the signal spectral shape. Our discussion here provides a crude approximation of the limits, and motivates a dedicated analysis by the experimental collaborations, analogous to the analysis done when setting limits on various massless majoron models with different spectral indices.
We can estimate the spectral effect on the limits by considering s/ √ b, where for background (b) we take the 2ν2β spectrum and for signal (s) the spectrum of 0ν2βφ with the selected value of m φ . The limit would deteriorate approximately We stress that the exercise in Fig. 3 is a rough approximation only: the official experimental analyses contain additional important sources of background from various radioactive contaminants, that are typically fit alongside with the signal. Nevertheless, in what follows we set estimated limits in the parameter space of G ee and m φ , using the max s/ √ b information from Fig. 3 together with the phase space suppression factor G(m φ )/G(0). The lower bound we take for T 1 2 is where T limit 1 2 (0) is the limit placed by the collaboration for the massless φ case. This translates to an upper bound on G ee that reads (here and elsewhere, where experimental limits are considered they always refer to the absolute value |G ee |) In Fig. 4 we plot the 90%CL upper bound on G ee , evaluated using Eq. (16). The region above the shaded bands is excluded by KamLAND-Zen [37] (blue, 136 Xe), NEMO-3 [35] (orange, 100 Mo), and NEMO-3 [36] (green, 150 Nd). The width of the band represents the uncertainties quoted by the collaborations for the massless φ limits. We stress that our approximate signal to background analysis implies larger uncertainty at m φ > 0. For comparison with other constraints, the dark shaded region above the horizontal black line shows the constraint from light meson decays [45]. The region above the coloured shaded bands is excluded by KamLAND-Zen [37] (blue, 136 Xe), NEMO-3 [35] (orange, 100 Mo), and NEMO-3 [36] (green, 150 Nd). The width of the band represents the uncertainties quoted by the collaborations regarding the massless φ limits; we stress that our approximate signal to background analysis implies larger uncertainty at m φ > 0. For comparison with other constraints, the dark shaded region above the horizontal black line shows the constraint from light meson decays [45].

III. IMPLICATIONS FOR SCALAR-MEDIATED NEUTRINO SELF-INTERACTIONS
Neutrino-neutrino interactions through light mediator exchange can cause observable features in the diffuse high-energy neutrino flux seen by IceCube and future neutrino telescopes [32,54,55]. This can occur if resonant s-channel scattering of a high-energy astrophysical neutrino with energy ν off the cosmic neutrino background (CνB) is possible, which in turn requires the mediator mass to match the center of mass energy (CME) of the collision, where m ν is the mass of the CνB neutrino participating in the collision. High-energy neutrino telescopes like IceCube and its planned Gen2 upgrade would be most sensitive to features in the astrophysical neutrino flux in the energy range from a few tens of TeV (below which the atmospheric background kicks in) up to around PeV (above which statistics fall off). The MeV mediator mass range is therefore of particular interest to this phenomenology. The detailed connection between 0ν2βφ limits and high-energy neutrino phenomenology is model-dependent. For concreteness, in the rest of this section we consider the framework of Ref. [32], where the coupling G αβ was proportional to the neutrino mass matrix. In particular, each real non-negative neutrino mass eigenvalue m νi is accompanied by a real non-negative value of G i = G αβ U αi U βi , while off-diagonal terms vanish in the mass basis, G αβ U αi U βj = 0 for i = j. The optical depth for resonant scattering, ν i ν i → φ → νν, proceeding through the scalar Φ in Eq. (6), considering neutrinos with observed energy ν , is bounded approximately by 5 In Eq. (18), the high-energy neutrino is assumed to have been emitted at redshift z. The relevant astrophysical emission is typically thought to be dominated around z ∼ 2 or so (see, e.g. [56][57][58][59][60]). Detectable effects at neutrino telescopes require If we are given the neutrino mass hierarchy (e.g. by upcoming neutrino oscillation experiments [61]) and the sum of neutrino masses (e.g. by cosmology [62,63]), then, using the measured PMNS values, we can relate the bound on G ee to a bound on G i for any i. In Fig. 5 we do this exercise, using neutrino oscillation parameters from Ref. [5]. Finally, in Fig. 6 we present our result for the 0ν2βφ in the parameter space relevant for high-energy neutrino telescopes. In and to the left of the red (blue) shaded area, G ee < 10 −4 (G ee < 10 −3 ) at 90%CL, where we have used the upper side (more conservative) of the KamLAND-Zen limit, rescaled from the m φ = 0 case as shown in Fig. 4. As explained above, the constraint on G ee can be readily converted to a constraint on neutrino optical depth, given information about the neutrino mass hierarchy and total mass. Diagonal lines show the observer frame neutrino energy which enters resonance s-channel scattering for scalar exchange with a CνB neutrino of mass m ν . In and above the green shaded region, the sum of neutrino masses exceeds 0.3 eV and is excluded by cosmological observations. The Q-values for 0ν2β in 100 Mo and 136 Xe are indicated by black dots at the bottom of the plot. We add these indicators to signify the effect of the phase space and signal to background considerations, which cause the naive 0ν2βφ constraint for a massless φ, |G ee | 10 −5 , to deteriorate in the massive φ case.

IV. CONCLUSIONS
A gauge-singlet scalar φ is expected to couple to two neutrinos, ν α ν β (α, β = e, µ, τ ) with couplings G αβ suppressed by v 2 /Λ 2 , where Λ is a scale of additional new physics. The G ee coupling can lead to neutrino-less double beta decay accompanied by scalar emission, 0ν2βφ. Experimental searches for 0ν2βφ have been conducted under the assumption that φ is the majoron, that is the massless Goldstone boson related to the spontaneous breaking of lepton number symmetry. It could, however, be the case that lepton number is explicitly broken, and φ is a massive scalar. The 0ν2βφ decay will proceed if G ee = 0 and the decay is kinematically allowed, m φ < Q.
If m φ is not much smaller than Q, then the bound on G ee extracted from the experimental upper bound on T −1 1 2 , the 0ν2βφ decay rate, is weakened compared to the massless majoron case. In this work, we obtained these bounds by considering the two main relevant effects: • The phase space factor G(m φ ) (see Fig. 2), which is suppressed compared to G(0); • The reduction in s/ √ b (see Fig. 3), the signal to root-background ratio, which is a consequence of the modification of the T 2β spectrum (T 2β is the sum of the kinetic energies of the two electrons).
The bounds on G ee for the massive scalar case are presented in Fig. 4.
The modification of the bounds from the massless majoron case to the massive scalar case are relevant for m φ = O(MeV). A scalar in this mass range which couples to neutrinos can have a strong effect on high energy astrophysical neutrinos observed by IceCube, as it mediates resonant scattering of these neutrinos on the cosmic neutrino background. Thus, 0ν2βφ constraints on massive scalars exclude part of the parameter space where the relevant features in the neutrino spectrum measured by IceCube may appear. This relation between 0ν2βφ and high energy neutrino phenomenology is presented in Fig. 6.
The exciting possibility to discover gauge-singlet scalars with mass in the MeV range via 0ν2β experiments (and the possible relation with high energy astrophysical neutrino observations) call for dedicated analyses by the experiments, where the effects of m φ = 0 are carefully taken into account.
Finally, for 0ν2β(φ * → 2ν) we have Here, the decay width of φ into two neutrinos is given by with |G| 2 = i G 2 i . Considering the virtual φ process, Eq. (A4), we see that: • Compared to the on-shell process 0ν2βφ, the decay rate for the off-shell process 0ν2β(φ * → 2ν) is suppressed by a factor ∼ 2|G| 2 . This suppression can be recognised as the product of (i) an additional final state phase space factor, (ii) an insertion of G 2 , and (iii) an over-all kinematical factor ∼ Q 4 /m 4 φ .
• Besides from the over-all suppression of the 0ν2β(φ * → 2ν) process, the spectral shape with respect to the outgoing electron energy is of the form ∼ (Q + 2m e − 1 − 2 ) 5 , which is, of course, just the spectral shape of the standard background process 2ν2β.
As a result of these features, the virtual φ decay mode 0ν2β(φ * → 2ν) cannot be constrained with current experiments.