Neutrinoless Double Beta Decay with Non-standard Majoron Emission

We present a novel mode of neutrinoless double beta decay with emission of a light Majoron-like scalar particle phi. We assume it couples via an effective 7-dimensional operator with a (V+A) lepton current and (V+-A) quark currents leading to a long-range contribution that is unsuppressed by the light neutrino mass. We calculate the total double beta decay rate and determine the fully differential shape for this mode. We find that future double beta decay searches are sensitive to scales of the order Lambda_NP ~ 1 TeV for the effective operator and a light scalar m_phi<0.2 MeV, based on ordinary double beta decay Majoron searches. The angular and energy distributions can deviate considerably from that of two-neutrino double beta decay, which is the main background. We point out possible ultraviolet completions where such an effective operator can emerge.


I. INTRODUCTION
Double beta decay processes are sensitive probes of physics beyond the Standard Model (SM). The SM process of two-neutrino double beta decay (2νββ) is the rarest process ever observed with half lives of order T 2νββ On the other hand, one or more exotic neutral particles may also be emitted, with a signature of anomalous missing energy beyond that expected in 2νββ decay. A well studied set of theories involve the emission of a scalar particle, called Majoron J. The first such proposed Majoron was a Goldstone boson associated with the spontaneous breaking of lepton number symmetry [1,2], coupling to a neutrino ν as g J ννJ, cf. Fig. 1 (left).
Current searches have a sensitivity of the order T 0νββJ 1/2 ∼ (10 −5 /g J ) 2 × 10 24 y. The term Majoron has been used in a wider sense, implying just a charge-neutral scalar particle (Goldstone boson or not) or vector particle [3]. Originally considered to be massless, it may also be a light particle [4][5][6] that can potentially be a Dark Matter candidate [7][8][9].
Searches for extra particles in double beta decay are crucial in understanding neutrinos. with the electrons and quarks. In this letter, we instead consider 0νββφ decay with emission of a light neutral scalar φ from a single effective dimension-7 operator of the form Λ −3 NP (ūOd)(ēOν)φ, cf. Fig. 1 (center), with the fermion currents having a different chiral structure from that in the SM. In the following, we will refer to the light scalar as "Majoron", independent of its origin. We determine the sensitivity to Λ NP and analyse the effect on the energy and angular distributions in comparison with 2νββ decay. We also comment on ultraviolet scenarios underlying the effective operator.

II. EFFECTIVE LONG-RANGE INTERACTIONS
We are interested in processes where right-and left-handed electrons are emitted along with a scalar φ considering as a first approach, only (V + A) and (V − A) currents. The effective Lagrangian can then be written as with the Fermi constant G F , the Cabbibo angle θ C , and the leptonic and hadronic currents j µ L,R =ēγ µ (1 ∓ γ 5 )ν and J µ L,R =ūγ µ (1 ∓ γ 5 )d, respectively. The proton mass m p is introduced in the exotic interactions as normalization to make the effective coupling constants φ RL and φ RR dimensionless, in analogy to the effective operator treatment of 0νββ decay [10,11]. In Eq. (1), we omit exotic operators with left-handed lepton currents; as in the standard long-range case, such contributions will be additionally suppressed by the small neutrino masses [11]. We instead focus on the process depicted in Fig. 1 where the SM (V − A) Fermi interaction, the first term in Eq. (1), meets one of the exotic operators. In this case, the momentum part in the numerator of the neutrino propagator contributes, rather than the mass. In Eq. (1) we consider the first generation electron and neutrino only. Generalizing to three flavours amounts to promoting the φ RX couplings to 3 × 3 matrices in generation space, ( φ RX ) αi (α = e, µ, τ , i = ν 1 , ν 2 , ν 3 ). The final decay rate will then be proportional to where U is the SM lepton mixing matrix.

III. DECAY RATE AND DISTRIBUTIONS
We base our calculation of the 0νββφ decay rate and kinematic distributions on Doi et al. [12]. The details of the calculation are given in the Appendix; here we outline the main features. Summing over all intermediate nuclear states N , the amplitude of 0 + I → 0 + F 0νββφ decay can be written as Here, X = L, R correspond to φ RL , φ RR , µ N = E N − E I + Q ββ /2 + m e with E I and E N the energies of the initial and intermediate nucleus, respectively. The energies of the two outgoing electrons and the Majoron are E 1,2 and E φ , respectively, and the available kinetic energy release Q ββ . The nucleon and lepton currents are defined as We consider that the internal neutrino propagates between the interaction points x and y with momentum q µ = (ω, q). From Eqs. (2) and (4) one can see explicitly the required antisymmetry of the amplitude under the exchange of the two electrons.
In Eq. (2), the Majoron energy E φ is added or subtracted depending on whether the electron labelled 1 or 2 is being emitted from the exotic operator. The Majoron makes a crucial difference, as E φ goes together with (E 1 − E 2 ) and not with the term proportional to the intermediate nuclei energy µ N as for an ordinary Majoron. A dependence on E φ will thus appear through the matrix element in addition to that through the phase space.
The differential decay rate for the 0 + → 0 + 0νββφ decay can then be written as [12] with the axial coupling g A of the nucleon and the radius R of the nucleus. The magnitudes of the electron spatial momenta are denoted p 1,2 and 0 ≤ θ ≤ π is the angle between the emitted electrons. In Eq. (5), the Majoron energy is determined as E φ = Q ββ + 2m e − Majoron decay (spectral index n = 1) are given for comparison.
decay rate can be found in the Appendix, where we show in detail the derivation of the differential decay rate. Therein we use the nuclear matrix elements listed in Table I and Coulomb-corrected relativistic electron wave functions.
From Eq. (5), the total decay rate and thus half life is calculated by performing the integral of a(E 1 , E 2 ) over all energies within the allowed phase space limits E 1 , E 2 ≥ 0 and E 1 + E 2 ≤ Q ββ + 2m e . In addition, we determine and discuss several distributions below. We will show results for the φ RL and φ RR versions of the effective operators where we consider only one of these to be present at a time. We assume the exotic φ Majoron to be massless in our calculations and comment on massive φ in the discussion below. For our numerical evaluation we focus on two isotopes: (i) 136 Xe, for which the KamLAND-Zen collaboration [13] currently provides the most stringent constraints; (ii) 82 Se used by NEMO-3 and the upcoming SuperNEMO experiments [14] that can measure the detailed electron topology.
For all experimental searches, the crucial distribution is with respect to the sum of the kinetic energies of the detected electrons. With the SM 2νββ decay as irreducible background to any exotic signal, it is important to calculate it precisely. In Fig. 2 (left), we compare the normalized total electron kinetic energy distribution of 0νββφ decay with that of 2νββ decay and ordinary 0νββJ Majoron decay (with spectral index n = 1) for the isotope 136 Xe. The distribution associated with φ RL is very similar to ordinary 0νββJ decay, while the introduction of a hadronic right-handed current in the φ RR term changes considerably the shape of the distribution. In both cases, the spectral index still corresponds to n = 1 with the characteristic onset near the kinematic endpoint. We emphasize that because of the different shape, a dedicated signal over background analysis is required to determine the experimental sensitivity on the effective parameters φ RL and φ RR precisely. NEMO-3 and SuperNEMO are able to measure the individual electron energies. In right-handed current scenarios without emission of a Majoron, the single energy distribution exhibits a distinctive valley-type shape. This occurs as the dominant term is for the corresponding RR term, as a result of the antisymmetry with respect to electron exchange. 1 In our case, depicted in Fig. 2 (right), part of the energy is being carried away by the Majoron, shifting the distribution towards lower electron energies and softening the characteristic valley-type distribution for φ RR . The distribution does not vanish for E 1 − m e = 1 2 Q ββ (as in the ordinary right-handed current case), but is still significantly different from that of ordinary Majoron emission.
The distribution with respect to both electron energies is depicted in Fig. 3

(top panel)
in the Appendix. It exhibits an even more pronounced difference between the φ RR mode and 2νββ. This may be used experimentally to improve the sensitivity through kinematic selection criteria, counteracting the effect of the less peaked total energy distribution, cf. One can use also angular correlations to distinguish between left-handed and righthanded currents [15,16] .70 (electrons are dominantly emitted collinearly) and k φ RR = −0.05 (electrons are emitted nearly isotropically) in our 0νββφ scenarios with φ RL and φ RR , respectively, for 82 Se. For comparison, the angular correlation factor for SM 2νββ decay is k 2νββ = −0.66 and k J = −0.80 for ordinary Majoron emission, i.e. the electrons are dominantly emitted back-to-back.
Finally, we estimate the sensitivity of existing and planned future double beta decay  searches on the effective coupling strength φ RL and φ RR of 0νββφ decay. We would like to emphasize again that due to the different total electron energy distribution, a dedicated signal over background analysis is required to determine the constraints precisely.
Experiments such as NEMO-3 and SuperNEMO can also improve their sensitivity due to the non-standard decay topology, especially for φ RR . As detailed in the Appendix, a requirement that any one electron has a kinetic energy of E i − m e > Q ββ /2 can for example reduce the 2νββ background by an order of magnitude. Here, we simply estimate the sensitivity by comparing our predictions for the 0νββφ decay half life T 1/2 = ln 2/Γ with the experimental constraints on ordinary (n = 1) Majoron emission. We use the most stringent limits for 82 Se by NEMO-3 [14] and for 136 Xe by KamLAND-Zen [13]. For future prospects, we estimate that experimental Majoron search sensitivities may reach T Se 1/2 ≈ 10 24 y (e.g. with the help of angular and energy selection cuts at SuperNEMO) and T Xe 1/2 ≈ 10 25 y. 2 The corresponding limits on φ RL and φ RR are shown in Table I, where only one effective operator is assumed to present at a time.

IV. DISCUSSION
Searches for Majorons or Majoron-like particles are a staple in double beta decay experiments. So far, they only cover the case where the neutrino involved couples via 2 The corresponding 0νββ decay sensitivities of the planned SuperNEMO [20] and nEXO experiments [21] may improve by O(100), but this requires an experimental approach that is essentially backgroundfree. This is not possible for Majoron emission with a continuous total electron energy spectrum. the SM (V − A) charged current interaction. This is clearly a well-motivated minimal choice but it is worthwhile to explore other scenarios. In this letter, we have discussed one such alternative where a Majoron-like particle φ is emitted from effective operators with (V + A) leptonic currents, cf. Fig 1 (center). The future sensitivities on the effective couplings φ RL and φ RR shown in Table I may be translated into effective operator scales Λ NP ≈ 1.3 TeV and 270 GeV, respectively, using 1/Λ 3 NP = φ RX G F cos θ C /( √ 2m p ). As noted before, we assume a massless φ in deriving these limits; they remain essentially unchanged for masses small compared to Q ββ , m φ 0.2 MeV and are of the same order for m φ 1 MeV, but will deteriorate as m φ → Q ββ (for a recent analysis in ordinary Majoron emission, see [9]).

An ultraviolet scenario generating the effective operators in Eq. (1) is suggested in
Left-Right symmetric models [22][23][24][25] where the SM W and ν are replaced by their right-  Fig. 1 (right). We can then identify leading to the estimate T Xe 1/2 where g R is the gauge coupling constant and θ R C the equivalent of the Cabibbo angle both associated with the SU (2) R of the Left-Right symmetric model. Alternatively, it is also possible to trigger the φ RL mode through the W R -W mixing θ. Its value is generically expected to be θ = κg R m 2 This is more stringent due to the better sensitivity on φ RL in Table I. Choosing the righthanded neutrino mass m N to be as low as 100 MeV is strictly speaking not allowed in the effective operator treatment which requires m N p F ≈ 100 MeV, but it may be more natural in a scenario where the mass of N is generated through the vacuum expectation value of φ, m N = y N φ . In fact, choosing m N to be smaller and abandoning the effective operator treatment may be more natural; the qualitative arguments should hold as above though a dedicated calculation of 0νββφ would be required. In addition, the contribution to 0νββφ via a heavy neutrino is expected to peak at m N ≈ p F with the above estimates applying to a good approximation [26].
Other ultraviolet completions do exist; to lowest dimension, the effective operator φ The Majoron may also couple derivatively, if originating as a Goldstone boson; this would increase the number of possible Lorentz-invariant combinations. Alternatively, if the exotic particle is a vector boson a µ , such as a dark photon, the fermion currents can couple to it via the vector field itself as well as its field strength tensor f µν . An even number of exotic neutral fermions χ may also be emitted but this would quickly increase the dimension of the corresponding effective operator. Instead, they may also originate from the internal neutrino via a dimension-6 operator of the form Λ −2 NP ννχχ [29]. Exploring such alternatives to the well-studied neutrinoless double beta decay is imperative in order to be able to draw reliable conclusions on the nature of neutrino mass generation. We here detail the computation of the amplitude and differential decay rate of the 0νββφ process. We follow the calculation of the standard long-range contributions presented in [12] and start from the effective Lagrangian with the SM charged current and the exotic 7-dimensional operators incorporating right-handed lepton currents and the Majoron φ, Here, G F is the Fermi constant, θ C is the Cabbibo angle and the leptonic and hadronic currents are defined as respectively.
To lowest order of perturbation, the amplitude for the process of 0 + I → 0 + F 0νββφ decay depicted in Fig. 1 (center) of the main text is The time-ordered product is expanded as with the chiral projectors defined as P L,R = 1 2 (1 ∓ γ 5 ). Using the neutrino propagator with momentum q and mass m ν , the highlighted term Ξ L µν (x, y) can be expressed as The amplitude needs to be antisymmetric under the exchange of the electrons e 1 and e 2 , and thus we generalize and E i is the energy of each electron.
We now perform the integral over the temporal variables. The integration over q 0 is straightforward by means of the residue theorem, with ω 2 = q 2 + m 2 ν . On the other hand, expanding the time-ordered product as and using the operator e iHt to extract the temporal dependence from the different wave functions, for example φ(y) = e iE φ y 0 φ(y), one can directly integrate over x 0 and y 0 obtaining the analogous expression to Eq. (C.2.19) in [12], . We anticipate the closure approximation and define the matrix element of the hadronic currents as In addition, the following properties under the exchange of position and electron energies were used in Eq. (A11), The integration over x 0 and y 0 in Eq. (A11) also provides the overall energy conservation We additionally assume that the Majoron φ is emitted predominantly in an S-wave configuration, φ(y) ≈ 1.
Considering the term between braces in Eq. (A11), one can write everything under the same exponential by interchanging x and y, It is furthermore useful to split the leptonic u L,R ρσ functions by separating out the part containing γ 5 as u L/R ρσ = 1 2 u ρσ ∓ u 5 ρσ . We then define These definitions become useful if one recalls that in the non-relativistic impulse approximation, the J L part of J LX ρσ acts on the n-th nucleon whereas the J X part acts on the m-th when performing the sum over all neutrons in the initial nucleus. The superscript ± in J ± ρσ thus indicates if the combination of currents is symmetric or antisymmetric under the interchange of m ↔ n. The same applies to F ± ρσ and F 5± ρσ . The closure approximation implies that the sum over all possible intermediate states is where H ωi are neutrino potentials defined as Now, the connection with the results of [12] can be done by contracting the leptonic and nuclear currents within the impulse approximation. The only change in our case is in the ω term, where the X and Y terms are functions of nuclear parameters and operators defined in Appendix C of [12]. The F α (5)± -terms are generated by the contraction of the hadronic and leptonic parts in Eqs. (A15)-(A17) factorizing out the dependence with the momentum q α from the leptonic part (see Eq. (C.2.25) in [12]). One trivially recovers the ω term in the expression (C.2.23) of [12] for E φ → 0.
Comparing Eq. (A20) with the results from [12], one can track the dependence with E φ in the decay rate down to Eq. (C.3.9) of [12]. The main change for 0 + I → 0 + F transitions is in the terms N 3 and N 4 where a contribution proportional to E φ appears explicitly, Here, α jk =Ã j (E 1 )Ã k (E 2 ) describe the Coulomb-corrected relativistic electron wave functions and ζ = 3αZ + (Q ββ + 2m e )R the correction of the electron P wave, with the fine structure constant α and the radius R and charge Z of the final state nucleus. The information about the electron wave functions is encoded iñ where we use an updated value from the same group [19].
with the Fermi factor where γ k = k 2 + (αZ) 2 , y = αZE/p and p = E 2 − m 2 e . In order to arrive at Eq. (A20) one should neglect the higher order terms E 2 12 , E 12 E φ and E 2 φ as they are suppressed with an extra denominator (ω+A i ) compared to Eq. (A19). The Z i terms are given in Eqs. (A25)-(A28) below and they contain the nuclear matrix elements and effective particle physics couplings. The Z i terms are the same as in [12], with the relevant couplings λ → φ RR and η → φ RL substituted. Note that the term with Z 1 in Eq. (C.3.9) from [12] related to the standard 0νββ decay disappears from Eq. (A21), as we are not considering the interaction L SM (x)L SM (y).
The above equations are valid when both φ RL and φ RR are present. For our numerical calculations, we use the Q ββ values and nuclear matrix elements M GT , χ F , etc. presented in Table II for 82 Se and 136 Xe. We use the following values for the remaining parameters: The differential decay rate for the 0 + → 0 + 0νββφ decay can then be written as [12] with Here, g A is the axial coupling of the nucleon and R is the radius of the nucleus. The magnitudes of the electron momenta are given by p i = E 2 i − m 2 e and 0 ≤ θ ≤ π is the angle between the emitted electrons. Throughout the above expressions, the Majoron energy is implicitly fixed by the electron energies as E φ = Q ββ + 2m e − E 1 − E 2 due to overall energy conservation.
The total decay rate Γ and the half life T 1/2 are then calculated as Γ = ln 2 The fully differential energy information is encoded in the normalized double energy This function, in terms of the kinetic energies normalized to the Q value, (E i − m e )/Q ββ , is plotted in the top row of Fig. 3 for the case of 0νββφ Majoron emission through φ RL (left) and φ RR (center) as well as for the SM 2νββ decay (right). The plots are for the isotope 82 Se but would be qualitatively similar for 136 Xe. As can be seen, the shapes depicted as contours are different between all three modes. Especially the φ RR exhibits an asymmetry in that one of the electrons takes the majority of the visible energy. If the individual electron energies can be measured, as e.g. in the NEMO-3 or SuperNEMO experiments, this can be exploited to enhance the signal over the 2νββ background. As an illustrating example, requiring that any one of the electrons in a signal event has a kinetic energy E i − m e > Q ββ /2 would reduce the 0νββφ-φ RR rate only by a factor of 2 but would suppress the 2νββ rate by a factor of 20. The distributions in Fig. 2 of the main text can be easily determined by appropriately integrating over dΓ dE 1 dE 2 . In addition to the energies, the angle between the electron momenta also contains useful information. The so-called angular correlation defined by is a function of the individual electron energies which can take values between −1 (the two electrons are dominantly emitted back-to-back) and +1 (the two electrons are dominantly emitted collinearly). For 82 Se it is plotted in the bottom row of Fig. 3 in the three modes of interest. As expected from angular momentum considerations, the electrons are dominantly emitted back-to-back in the SM 2νββ decay with (V − A) lepton currents, α < 0 for all energies. For φ RL , they are dominantly emitted collinearly, α > 0 for all energies. In the case of φ RR , the behaviour is complex due to the asymmetry of the amplitude under the exchange of electrons and nuclear recoil effects. The correlation α changes sign, with α > 0 when any one electron has a kinetic energy E i − m e > Q ββ /2 but α < 0 when both electrons each have a kinetic energy E i − m e < Q ββ /2. Note that Fig. 3 provides the full kinematical information in each mode; all measurable quantities can be constructed from these distributions.
As discussed in the main text, averaged over all energies, the electrons are actually emitted almost isotropically for φ RR . This is quantified by integrating Eq. (A31) over all energies analogous to Eq. (A34) to yield the angular distribution with the average angular correlation factor k.