---
abstract: '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 charachterizing 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.  The spectral shape of high-energy astrophysical neutrinos could exhibit features due to resonant $\nu\nu\to\phi\to\nu\nu$ scattering.  Such features fall within the sensitivity range of IceCube-like experiments, if $m_\phi$ is of order MeV, making $0\nu2\beta\phi$ a key complimentary laboratory constraint on the scenario.  Our results motivate a dedicated analysis by $0\nu2\beta$ collaborations, analogous to the dedicated analyses targeting massless majoron models.' 
author:
- Kfir Blum
- Yosef Nir
- Michal Shavit
bibliography:
- 'ref.bib'
title: 'Neutrinoless double-beta decay with massive scalar emission'
---

Introduction
============

Neutrinoless double beta ($0\nu2\beta$) decay [@Doi:1985dx; @Rodejohann:2011mu; @Vergados:2012xy; @Pas:2015eia], (A,Z)(A,Z+2)+2e\^[-]{}, is a lepton number violating process.  It is sensitive to the neutrino mass parameter m\_[ee]{}=|\_i m\_[\_i]{} U\_[ei]{}\^2|, where $m_{\nu_i}$ ($i=1,2,3$) are the neutrino masses and $U$ is the lepton mixing matrix [@Patrignani:2016xqp].  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 [@Weinberg:1979sa][^1][^2], \[eq:dimfive\] [L]{}\_[d=5]{}=-(HL\_) (HL\_), where $H$ is the Higgs doublet field and $L_\alpha$ ($\alpha=e,\mu,\tau$) are the lepton doublet fields, leads to neutrino masses, m\_= , with $v=246$ GeV.

We do not know the beyond-SM origin of the dimension-five terms in Eq.  (\[eq:dimfive\]).  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 $\phi$, then the dimension-six terms \[eq:dimsix\] [L]{}\_[d=6]{}=-(HL\_) (HL\_) are possible.  The dimension-six terms lead to Yukawa couplings of $\phi$ to neutrinos, \[eq:Lint\]Ł&&-\_\_\_+[h.c.]{},\
\[eq:Gab\]&=&.  If the mass of the $\phi$ particle is less than the $Q$-value of the $(A,Z)\to (A,Z+2)$ transition, $m_\phi<Q$, then the ${\mathcal{G}}_{ee}$ coupling leads to a decay where $0\nu2\beta$ is accompanied by on-shell $\phi$ emission ($0\nu2\beta\phi$), (A,Z)(A,Z+2)+2e\^-+.  A well known framework that leads to Eq.  (\[eq:dimsix\]) and to the decay mode $0\nu2\beta\phi$ is that of majoron models [@PhysRevD.25.774; @Georgi:1981pg; @Gelmini:1980re; @Vergados:1981gk; @Chikashige:1980ui; @Burgess:1992dt], where $\phi$ 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 $\Lambda$ appearing in Eq.  (\[eq:dimfive\]) is promoted to a dynamical field, and the phase of this field is associated with $\phi$. In such models, (i) the $\phi$ particle is massless, and (ii) the terms of Eq.  (\[eq:dimfive\]) and Eq.  (\[eq:dimsix\]) are related, leading to ${\mathcal{G}}= m_\nu/\Lambda$. For high-scale seesaw models, with ${\mathcal{Z}}=\mathcal{O}(1)$, the seesaw scale is $\Lambda\sim10^{14}$ GeV, leading to ${\mathcal{G}}\sim10^{-24}$. As we review in Sec.  \[sec:0n2bphi\], such tiny coupling is some 20 orders of magnitude below the reach of $0\nu2\beta\phi$ searches. 

In other scenarios, like the inverse-seesaw models of Ref.  [@PhysRevLett.56.561; @PhysRevD.34.1642; @BernabŽu1987303] (see Ref.  [@Weiland:2013wha] for a review), neutrino masses arise from effective dimension six terms.  Namely, instead of $1/\Lambda$ in Eq.  (\[eq:dimfive\]) we have $\mu/\Lambda^2$, where a technically natural hierarchy $\mu\ll\Lambda$ 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 $\mu$ to a field $\phi$ with $\mu=\langle\phi\rangle$. In this case, ${\mathcal{G}}= m_\nu/\mu$ and if $\mu$ is small enough, $0\nu2\beta\phi$ could be observable.  If lepton number is broken spontaneously by $\langle\phi\rangle$, then the $\phi$ 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 $\phi$ could be light but not massless [@Rothstein:1992rh; @Akhmedov:1992hi].  In addition, if the explicit lepton number violation (LNV) dominates the neutrino mass, then also the relation between ${\mathcal{G}}$ and $m_\nu$ 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 [@Berryman:2018ogk] for a recent study.  While $0\nu2\beta$ experimental collaborations often consider various massless majoron models, such as different forms of the majoron-neutrino couplings and multi-majoron final-state processes, none of these searches considered the scenario of a massive majoron, $m_\phi\sim$ MeV, of the same order as the $Q$-value of the $0\nu2\beta$ reaction.  In this paper we consider this parameter region[^3] and estimate $0\nu2\beta\phi$ constraints for the case of $m_\phi$ of order MeV. As we show, the constraints are affected not only by kinematical phase space suppression near $m_\phi\sim Q$, but also by a change in the signal to background ratio characterising the search.  As a result, $0\nu2\beta\phi$ constraints for $m_\phi>0$ diminish significantly below the reaction threshold.  Our results motivate a dedicated analysis by $0\nu2\beta$ collaborations, analogous to the dedicated analyses targeting different massless majoron models. 

The constraint on massive $\phi$ emission in $0\nu2\beta\phi$ has phenomenological implications, which we illustrate focusing on high-energy neutrino telescopes.  Light scalar fields coupled to neutrinos were considered as mediators of anomalous neutrino self-interactions in many other works.  Refs.  [@Goldberg:2005yw; @Baker:2006gm] studied the effect of light scalar exchange on the energy spectrum of $\sim$10 MeV neutrinos from core-collapse supernovae (see also [@Farzan:2014gza] where supernovae neutrinos scatter on dark matter).  Vector boson or massless majoron exchange were considered in [@Kolb:1987qy; @Keranen:1997gz; @Hooper:2007jr].  Refs.  [@Chacko:2003dt; @Hall:2004yg; @Davoudiasl:2005ks; @Friedland:2007vv] discussed the relation of anomalous neutrino interactions to low-scale neutrino mass generation, focusing on spontaneously broken global and gauged lepton number.  Ref.  [@Blum:2014ewa] 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.  [@Farzan:2018gtr] considered light scalar exchange in coherent neutrino-nucleus scattering. 

Before we turn into concrete calculations, let us emphasize that while $0\nu2\beta$ is a LNV process, $0\nu2\beta\phi$ could be lepton number conserving (LNC).  It could be therefore that the latter is strongly enhanced compared to the former.  Explicitly, for $m_\phi\ll Q$, we have \[eq:0n2bphi20n2bphi\] \~60()\^2, where we used the fact that $m_{ee}\lsim0.1$ eV [@Barabash:2014uqa], with $Q\sim$MeV. As we review in the next section, $0\nu2\beta\phi$ searches have reached a limit $|{\mathcal{G}}_{ee}|\lesssim10^{-5}$ (for massless $\phi$).  The reason that $\Gamma_{0\nu\beta\beta\phi}\gg\Gamma_{0\nu\beta\beta}$, as shown by Eq.  (\[eq:0n2bphi20n2bphi\]), 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\nu2\beta\phi$ compared with $0\nu2\beta$. In what follows we will see an anaolgous effect deteriorating the sensitivity to $m_\phi>0$ compared to the $m_\phi=0$ case. 

Finally, note that when $m_\phi>Q$, the on-shell process $0\nu2\beta\phi$ is kinematically blocked, but the off-shell process $0\nu2\beta(\phi^*\to2\nu)$, where a virtual $\phi$ is emitted and decays to two neutrinos, is always allowed.  However, compared to the on-shell process $0\nu2\beta\phi$ (when allowed), the off-shell $\phi$ process is strongly suppressed by a factor $\sim\frac{2|\mathcal{G}|^2}{15(4\pi)^2}\,\frac{Q^4}{m_\phi^4}\lesssim10^{-7}\left(\frac{\mathcal{G}}{10^{-2}}\right)^2$. In addition, the spectral shape with respect to the outgoing electron energy is similar to that of the standard background process $2\nu2\beta$. These features are explained in App.  \[app:virt\].  As a result of these features, the virtual $\phi$ decay mode $0\nu2\beta(\phi*\to2\nu)$ cannot be constrained with current experiments, and we limit our attention to on-shell $0\nu2\beta\phi$.

Neutrinoless double-beta decay with massive scalar emission {#sec:0n2bphi}
===========================================================

From the list of $0\nu2\beta$ experiments surveyed in [@Barabash:2014uqa], NEMO-3 [@Arnold:2013dha] using $^{100}$Mo has the highest $Q$-value, $Q\approx3.03$ MeV. Recent work by NEMO-3 allowed them to surpass this record using $^{150}$Nd [@Arnold:2016qyg], with $Q\approx3.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 [@Gando:2012pj] using $^{136}$Xe, which has a somewhat lower value of $Q\approx2.5$ MeV. For these reasons – sensitivity to the highest $m_\phi$, and current best sensitivity to massless $\phi$ – 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 [@GERDA:2018zzh; @Aalseth:2017btx], $^{82}$Se [@Azzolini:2018dyb], and $^{130}$Te [@Andringa:2015tza; @Alduino:2017ehq].  These isotopes yield comparable, although (currently) somewhat weaker constraints. 

Refs.  [@Arnold:2013dha], [@Arnold:2016qyg], and [@Gando:2012pj] provide 90%CL bounds on the massless majoron scenario, equivalent to $|{\mathcal{G}}_{ee}|<(1.6-4.1)\times10^{-5}$, $|{\mathcal{G}}_{ee}|<(3.8-14.4)\times10^{-5}$, and $|{\mathcal{G}}_{ee}|<(0.4-1.0)\times10^{-5}$, respectively.  These bounds are stronger than other constraints in the literature such as those arising from light meson decay (see, e.g. [@Barger:1981vd; @Lessa:2007up; @Pasquini:2015fjv]) and from cosmological and astrophysical considerations [@Blum:2014ewa].  While the massless majoron bounds [@Doi:1987rx; @Rath:2016saq; @Hirsch:1995in; @Pas:1999fc; @Pas:2000vn] coincide with our model for $m_\phi\ll Q$, to our knowledge a study of the kinematical region $m_\phi\sim Q\sim$ MeV has not been done and we consider this region in what follows. 

Searches for $0\nu2\beta\phi$ constrain the half-life time for the decay, $T_{\frac{1}{2}}$, which in our model can be approximately decomposed as [@Doi:1985dx; @Doi:1988zz; @Doi:1987rx] (see, e.g. [@Kotila:2015ata] for a recent account) \[eq:T12\]T\_\^[-1]{}&&|\_[ee]{}|\^2||\^2G(m\_),where ${\mathcal{G}}_{ee}$ is defined in Eq.  (\[eq:Gab\]), $\M$ is a dimensionless nuclear matrix element (NME), and $G(m_\phi)$ is a kinematical phase space factor, conventionally expressed in units of yr$^{-1}$. The main effect of massive $\phi$ emission is to modify the phase space factor [@Burgess:1992dt], $G(0)\to G(m_\phi)$: \[eq:G\] G(m\_)&& dE\_1dE\_2 E\_1p\_1E\_2p\_2((Q+2m\_e-E\_1-E\_2)\^2-m\_\^2)\^a(E\_1,E\_2) .The outgoing electrons Coulomb factor, encoded in $a(E_1,E_2)$, is given in Refs.  [@Doi:1985dx; @Doi:1987rx].  For the range of electron momenta of interest, $p_i>0.5$ MeV, $a(E_1,E_2)$ can be factorized as $a(E_1,E_2)\approx F_0(E_1)F_0(E_2)$. For simplicity, for most of the numerical results in this work we use the non-relativistic approximation, \[eq:PRA\]F(v)&&.  We have checked that our results are not affected significantly when using the relativistic expressions for the electron wave function. 

The $0\nu2\beta\phi$ constraint on ${\mathcal{G}}_{ee}$ deteriorates as the scalar mass $m_\phi$ 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_\phi)/G(0)\to0$ as $m_\phi\to Q$.

2.   As $m_\phi$ is increased, the visible final state electrons kinetic energy T\_[2]{}=E\_1+E\_1-2m\_eis pushed to lower values, because the usual massless majoron phase space factor $(Q-T_{2\beta})$ is replaced by $((Q-T_{2\beta})^2-m_\phi^2)^{\frac{1}{2}}$. This brings the distribution of $T_{2\beta}$ to overlap more with the $T_{2\beta}$ distribution of the irreducible Standard Model $2\nu2\beta$ background, leading to smaller signal to background ratio in the experimentally relevant energy range. 

Fig.  \[fig:spectrum\] illustrates both of the above points, showing the visible electron $T_{2\beta}$ spectrum for $0\nu2\beta\phi$ for different values of $m_\phi$ and comparing to the $2\nu2\beta$ case (in arbitrary normalisation).  This plot considers $^{100}$Mo as an example.  Note that the experimental analyses typically impose a lower energy cutoff of $T_{2\beta}\gtrsim0.5$ MeV.

![ Energy spectrum (sum of kinetic energies of the two electrons) for $0\nu2\beta\phi$. The grey dashed curve (solid vertical line) shows the spectrum for $2\nu2\beta$ ($0\nu2\beta$) decay.  The normalisation is arbitrary, apart from the relative normalisation of the $0\nu2\beta\phi$ curves which follow the phase space factor.[]{data-label="fig:spectrum"}](spectrum){width="80.00000%"}

Proceeding further, we consider first the simple phase space suppression encoded in $G(m_\phi)/G(0)$. In the Primakoff-Rosen approximation (PRA) [@0034-4885-22-1-305], we can calculate this ratio analytically.  We find: \[eq:Gm2G0PRA\] &=&\
&&,where $\bar{m}_{\phi}= m_{\phi}/Q$ and $\bar{m}_{e}= m_{e}/Q$. Taking the limit $\bar{m}_{\phi}\rightarrow0$, we recover the phase space factor for massless majoron (see, e.g. [@Georgi:1981pg]).  The phase space suppression factor is shown in Fig.  \[fig:Gm2G0\] for $^{100}$Mo (blue) and $^{136}$Xe (red).  The analytical PRA result Eq.  (\[eq:Gm2G0PRA\]) 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_\phi\sim1$ MeV. In the case of $^{136}$Xe ($^{100}$Mo), for $m_\phi\approx2$ MeV (2.5 MeV), the decay width drops to $\sim$1% of its value for the massless majoron case.  This means that the limit on ${\mathcal{G}}_{ee}$ becomes weaker than the massless majoron limit by a factor of at least $\sim$10 at that point, even ignoring the signal to background ratio deterioration effect that we discuss later on. 

![ Phase space suppression factor for massive $\phi$, vs. $m_\phi$, for $^{100}$Mo (blue) and $^{136}$Xe (red).  The analytical PRA result Eq.  (\[eq:Gm2G0PRA\]) is shown in solid lines.  Numerical computation using the full relativistic electron wave function is shown by dots.[]{data-label="fig:Gm2G0"}](Gm2G0){width="55.00000%"}

We next consider the varying signal to background ratio at varying $m_\phi$. 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/\sqrt{b}$, where for background ($b$) we take the $2\nu2\beta$ spectrum and for signal ($s$) the spectrum of $0\nu2\beta\phi$ with the selected value of $m_\phi$. The limit would deteriorate approximately in proportion to the maximum value of $s/\sqrt{b}$, attained in the region $T_{2\beta}\gtrsim0.5$ MeV analyzed by the experiments.  This information is presented in Fig.  \[fig:s2rb\].  The x-axis of the plot gives the scalar mass $m_\phi$. The y-axis gives the value of $s/\sqrt{b}$, where we calculate this ratio at the point[^4] in $T_{2\beta}$ that maximizes $s$ within the range $T_{2\beta}>0.5$ MeV. To make the interpretation easier, we normalize max $s/\sqrt{b}$ obtained at any value of $m_\phi$ to the value of max $s/\sqrt{b}$ obtained for $m_\phi=0$. Fig.  \[fig:s2rb\] then implies that the $^{136}$Xe limit on $T_{\frac{1}{2}}$ for $m_\phi=1$ MeV, for example, should be about a factor of 10 weaker than the limit on $T_{\frac{1}{2}}$ for $m_\phi=0$.

![ Normalised signal to root-background, as function of $m_\phi$.[]{data-label="fig:s2rb"}](s2rootb){width="75.00000%"}

We stress that the exercise in Fig.  \[fig:s2rb\] 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 $\mathcal{G}_{ee}$ and $m_\phi$, using the max $s/\sqrt{b}$ information from Fig.  \[fig:s2rb\] together with the phase space suppression factor $G(m_\phi)/G(0)$. The lower bound we take for $T_{\frac{1}{2}}$ is T\^[limit]{}\_(m\_)&=&T\^[limit]{}\_(0),where $T^{\rm limit}_{\frac{1}{2}}(0)$ is the limit placed by the collaboration for the massless $\phi$ case.  This translates to an upper bound on $\mathcal{G}_{ee}$ that reads (here and elsewhere, where experimental limits are considered they always refer to the absolute value $|{\mathcal{G}}_{ee}|$) \[eq:cstrtGee\]\_[ee]{}\^[limit]{}(m\_)&=&\_[ee]{}\^[limit]{}(0).  In Fig.  \[fig:GeeLim\] we plot the 90%CL upper bound on ${\mathcal{G}}_{ee}$, evaluated using Eq.  (\[eq:cstrtGee\]).  The region above the shaded bands is excluded by KamLAND-Zen [@Gando:2012pj] (blue, $^{136}$Xe), NEMO-3 [@Arnold:2013dha] (orange, $^{100}$Mo), and NEMO-3 [@Arnold:2016qyg] (green, $^{150}$Nd).  The width of the band represents the uncertainties quoted by the collaborations for the massless $\phi$ limits.  We stress that our approximate signal to background analysis implies larger uncertainty at $m_\phi>0$. For comparison with other constraints, the dark shaded region above the horizontal black line shows the constraint from light meson decays [@Pasquini:2015fjv]. 

![ 90%CL upper bound on ${\mathcal{G}}_{ee}$. The region above the coloured shaded bands is excluded by KamLAND-Zen [@Gando:2012pj] (blue, $^{136}$Xe), NEMO-3 [@Arnold:2013dha] (orange, $^{100}$Mo), and NEMO-3 [@Arnold:2016qyg] (green, $^{150}$Nd).  The width of the band represents the uncertainties quoted by the collaborations regarding the massless $\phi$ limits; we stress that our approximate signal to background analysis implies larger uncertainty at $m_\phi>0$. For comparison with other constraints, the dark shaded region above the horizontal black line shows the constraint from light meson decays [@Pasquini:2015fjv].[]{data-label="fig:GeeLim"}](GeeLim2){width="75.00000%"}

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 [@Ng:2014pca; @Ioka:2014kca; @Blum:2014ewa].  This can occur if resonant s-channel scattering of a high-energy astrophysical neutrino with energy $\epsilon_\nu$ off the cosmic neutrino background (C$\nu$B) is possible, which in turn requires the mediator mass to match the center of mass energy (CME) of the collision, m\_&=&=2()\^()\^ [MeV]{},where $m_\nu$ is the mass of the C$\nu$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\nu2\beta\phi$ limits and high-energy neutrino phenomenology is model-dependent.  For concreteness, in the rest of this section we consider the framework of Ref.  [@Blum:2014ewa], where the coupling ${\mathcal{G}}_{\alpha\beta}$ was proportional to the neutrino mass matrix.  In particular, each real non-negative neutrino mass eigenvalue $m_{\nu_i}$ is accompanied by a real non-negative value of ${\mathcal{G}}_i={\mathcal{G}}_{\alpha\beta}U_{\alpha i}U_{\beta i}$, while off-diagonal terms vanish in the mass basis, ${\mathcal{G}}_{\alpha\beta}U_{\alpha i}U_{\beta j}=0$ for $i\neq j$. The optical depth for resonant scattering, $\nu_i\nu_i\to\phi\to\nu\nu$, proceeding through the scalar $\Phi$ in Eq.  (\[eq:Lint\]), considering neutrinos with observed energy $\epsilon_\nu$, is bounded approximately by [^5] \[eq:optd\]\_[\_i\_i]{}&&2()\^2()\^[-2]{}()\^3(-\_)((1+z)\_-).In Eq.  (\[eq:optd\]), 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\sim2$ or so (see, e.g. [@Waxman:1998yy; @Bahcall:1999yr; @Waxman:2013zda; @Halzen:2013dva; @Meszaros:2014tta]).  Detectable effects at neutrino telescopes require $\tau=\mathcal{O}(1)$, implying ${\mathcal{G}}_i\gsim10^{-4}$.

In this model, because ${\mathcal{G}}_i\propto m_{\nu_i}$, we have \_i=|\_[ee]{}|||.If we are given the neutrino mass hierarchy (e.g. by upcoming neutrino oscillation experiments [@Capozzi:2017ipn]) and the sum of neutrino masses (e.g. by cosmology [@Lesgourgues:2006nd; @Abazajian:2016yjj]), then, using the measured PMNS values, we can relate the bound on ${\mathcal{G}}_{ee}$ to a bound on ${\mathcal{G}}_i$ for any $i$. In Fig.  \[fig:Gi2Gee\] we do this exercise, using neutrino oscillation parameters from Ref.  [@Patrignani:2016xqp]. 

![ The ratio ${\mathcal{G}}_i/|{\mathcal{G}}_{ee}|$ for the different neutrino mass eigenstates $i=1,2,3$, for normal and inverted hierarchy, as function of the sum of neutrino mass $\sum m_\nu$.[]{data-label="fig:Gi2Gee"}](Gi2Gee){width="60.00000%"}

Finally, in Fig.  \[fig:mnu\] we present our result for the $0\nu2\beta\phi$ in the parameter space relevant for high-energy neutrino telescopes.  In and to the left of the red (blue) shaded area, ${\mathcal{G}}_{ee}<10^{-4}$ (${\mathcal{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_\phi=0$ case as shown in Fig.  \[fig:GeeLim\].  As explained above, the constraint on ${\mathcal{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$\nu$B neutrino of mass $m_\nu$. 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\nu2\beta$ 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\nu2\beta\phi$ constraint for a massless $\phi$, $|{\mathcal{G}}_{ee}|\lesssim10^{-5}$, to deteriorate in the massive $\phi$ case. 

![ $0\nu2\beta\phi$ constraints presented in the $m_\phi-m_\nu$ parameter space, relevant for high-energy neutrino phenomenology.  In and to the left of the red (blue) shaded area, ${\mathcal{G}}_{ee}<10^{-4}$ (${\mathcal{G}}_{ee}<10^{-3}$) at 90%CL.  Diagonal lines show the observer frame neutrino energy which enters resonance s-channel scattering for scalar exchange with a C$\nu$B neutrino of mass $m_\nu$. 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\nu2\beta$ in $^{100}$Mo and $^{136}$Xe are indicated by black dots at the bottom of the plot.[]{data-label="fig:mnu"}](m_nu){width="75.00000%"}

Conclusions
===========

A gauge-singlet scalar $\phi$ is expected to couple to two neutrinos, $\nu_\alpha\nu_\beta$ ($\alpha,\beta=e,\mu,\tau$) with couplings ${\mathcal{G}}_{\alpha\beta}$ suppressed by $v^2/\Lambda^2$, where $\Lambda$ is a scale of additional new physics.  The ${\mathcal{G}}_{ee}$ coupling can lead to neutrino-less double beta decay accompanied by scalar emission, $0\nu2\beta\phi$. Experimental searches for $0\nu2\beta\phi$ have been conducted under the assumption that $\phi$ 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 $\phi$ is a massive scalar.  The $0\nu2\beta\phi$ decay will proceed if ${\mathcal{G}}_{ee}\neq0$ and the decay is kinematically allowed, $m_\phi<Q$.

If $m_\phi$ is not much smaller than $Q$, then the bound on ${\mathcal{G}}_{ee}$ extracted from the experimental upper bound on $T_{\frac{1}{2}}^{-1}$, the $0\nu2\beta\phi$ 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_\phi)$ (see Fig.  \[fig:Gm2G0\]), which is suppressed compared to $G(0)$;

-   The reduction in $s/\sqrt{b}$ (see Fig.  \[fig:s2rb\]), the signal to root-background ratio, which is a consequence of the modification of the $T_{2\beta}$ spectrum ($T_{2\beta}$ is the sum of the kinetic energies of the two electrons). 

The bounds on ${\mathcal{G}}_{ee}$ for the massive scalar case are presented in Fig.  \[fig:GeeLim\]. 

The modification of the bounds from the massless majoron case to the massive scalar case are relevant for $m_\phi={\cal O}({\rm 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\nu2\beta\phi$ 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\nu2\beta\phi$ and high energy neutrino phenomenology is presented in Fig.  \[fig:mnu\]. 

The exciting possibility to discover gauge-singlet scalars with mass in the MeV range via $0\nu2\beta$ experiments (and the possible relation with high energy astrophysical neutrino observations) call for dedicated analyses by the experiments, where the effects of $m_\phi\neq0$ are carefully taken into account. 

We thank Laura Baudis and Jacobo Lopez-Pavon for useful correspondence.  KB is incumbent of the Dewey David Stone and Harry Levine career development chair, and is supported by grant 1937/12 from the I-CORE program of the Planning and Budgeting Committee and the Israel Science Foundation and by grant 1507/16 from the Israel Science Foundation.  YN is the Amos de-Shalit chair of theoretical physics, and is supported by grants from the Israel Science Foundation (grant number 394/16), the United States-Israel Binational Science Foundation (BSF), Jerusalem, Israel (grant number 2014230), the I-CORE program of the Planning and Budgeting Committee and the Israel Science Foundation (grant number 1937/12), and the Minerva Foundation. 

Phase space for different decay modes {#app:virt}
=====================================

In the main text, it was convenient to decompose the inverse decay lifetime for the $0\nu2\beta\phi$ process as (reproducing Eq.  (\[eq:T12\])) \[eq:T12mtxt\]T\_\^[-1]{}&&|\_[ee]{}|\^2||\^2G(m\_),where $G(m_\phi)$ is the kinematical factor, encoding the important features of the $0\nu2\beta\phi$ mode, and $\mathcal{M}$ is the nuclear matrix element. 

Here we reproduce this decomposition in some more detail, in order to compare the different processes $0\nu2\beta$, $0\nu2\beta\phi$, and the off-shell scalar process $0\nu2\beta(\phi*\to2\nu)$. In the latter process, we refer to the case where no neutrinos are emitted from the weak current terms in the direct nuclear decay calculation (as in $0\nu2\beta\phi$), but where the scalar mass $m_\phi$ is larger than the reaction threshold, leading to off-shell scalar diagram that can be decomposed as $0\nu2\beta\phi^*$ with virtual $\phi$ decaying via $\phi^*\to\nu\nu$.

For $0\nu2\beta$, we have \_[MMee]{}&&|\_[MMee]{}|\^2k\_1k\_2(2)\^4\^[(4)]{}(P+k\_1+k\_2-P)a(\_1,\_2)\
&=&|\_[MMee]{}|\^2d\_1k\_1d\_2k\_2\_1\_2(Q+2m\_e-\_1-\_2)a(\_1,\_2).  For $0\nu2\beta\phi$, with on-shell $\phi$, we have \_[MMee]{}&&|\_[MMee]{}|\^2k\_1k\_2(2)\^4\^[(4)]{}(P+k\_1+k\_2+k-P)a(\_1,\_2)\
&=&|\_[MMee]{}|\^2d\_1k\_1d\_2k\_2\_1\_2 ((Q+2m\_e-\_1-\_2)\^2-m\_\^2)\^a(\_1,\_2),in which $G(m_\phi)$ of Eq.  (\[eq:G\]) can be identified. 

Finally, for $0\nu2\beta(\phi*\to2\nu)$ we have \[eq:virt\]\_[MMee(\^\*)]{}&&|\_[MMee]{}|\^2k\_1k\_2\
&&(2)\^4\^[(4)]{}(P+k\_1+k\_2+k\_[\_1]{}+k\_[\_2]{}-P)a(\_1,\_2)\
&=& |\_[MMee]{}|\^2d\_1k\_1d\_2k\_2\_1\_2a(\_1,\_2)\_0\^d(Q+2m\_e-\_1-\_2-)\
&&\_0\^1 dx\
&& |\_[MMee]{}|\^2d\_1k\_1d\_2k\_2\_1\_2a(\_1,\_2)(1+()).\
Here, the decay width of $\phi$ into two neutrinos is given by \_&=&k\_[\_1]{}k\_[\_2]{}(2)\^4\^[(4)]{}(k\_[\_1]{}+k\_[\_2]{}-k\_)= ,with $|\mathcal{G}|^2=\sum_i\mathcal{G}_i^2$.

Considering the virtual $\phi$ process, Eq.  (\[eq:virt\]), we see that:

-   Compared to the on-shell process $0\nu2\beta\phi$, the decay rate for the off-shell process $0\nu2\beta(\phi*\to2\nu)$ is suppressed by a factor $\sim\frac{2|\mathcal{G}|^2}{15(4\pi)^2}\,\frac{Q^4}{m_\phi^4}\sim10^{-7}\left(\frac{\mathcal{G}}{10^{-2}}\right)^2\frac{Q^4}{m_\phi^4}$. This suppression can be recognised as the product of (i) an additional final state phase space factor, (ii) an insertion of $\mathcal{G}^2$, and (iii) an over-all kinematical factor $\sim Q^4/m_\phi^4$.

-   Besides from the over-all suppression of the $0\nu2\beta(\phi*\to2\nu)$ process, the spectral shape with respect to the outgoing electron energy is of the form $\sim \left(Q+2m_e-\epsilon_1-\epsilon_2\right)^5$, which is, of course, just the spectral shape of the standard background process $2\nu2\beta$.

As a result of these features, the virtual $\phi$ decay mode $0\nu2\beta(\phi*\to2\nu)$ cannot be constrained with current experiments. 

[^1]: Repeated flavour indices are summed-over, and the bracket $(HL)$ denotes contraction to an SU(2) singlet. 

[^2]: See Ref.  [@Cirigliano:2017djv] for a recent discussion of $0\nu2\beta$ in the SM effective field theory. 

[^3]: We note that Ref.  [@1993PhLB..308...85C] considered neutrino-less double-beta decay with emission of a massive vector boson. 

[^4]: In this calculation we do not account for the energy resolution of KamLAND-Zen and of NEMO-3 [@Gando:2012pj; @Arnold:2013dha].  We checked that this does not modify the results appreciably. 

[^5]: Eq.  (\[eq:optd\]) is an upper bound to $\tau_{\nu\nu\to\phi\to\nu\nu}$ because it assumes that the scalar $\phi$ can only decay back to neutrinos via Eq.  (\[eq:Lint\]), minimizing its width; if other decay modes are possible for $\phi$, the resonant scattering cross section and the optical depth are suppressed by $\Gamma_{\phi\to\nu\nu}/\Gamma_\phi$.
