Two-Neutrino Double Beta Decay with Sterile Neutrinos

Usually considered a background for experimental searches for the hypothetical neutrinoless double beta decay process, two-neutrino double beta decay nevertheless provides a complementary probe of physics beyond the Standard Model. In this paper we investigate how the presence of a sterile neutrino, coupled to the Standard Model either via a left-handed or right-handed current, affects the energy distribution and angular correlation of the outgoing electrons in two-neutrino double beta decay. We pay particular attention on the behaviour of the energy distribution at the kinematic endpoint and we estimate the current limits on the active-sterile mixing and effective right-handed coupling using current experimental data as a function of the sterile neutrino mass. We also investigate the sensitivities of future experiments. Our results complement the corresponding constraints on sterile neutrinos from single beta decay measurements in the 0.1 - 10 MeV mass range.

Abstract: Usually considered a background for experimental searches for the hypothetical neutrinoless double beta decay process, two-neutrino double beta decay nevertheless provides a complementary probe of physics beyond the Standard Model. In this paper we investigate how the presence of a sterile neutrino, coupled to the Standard Model either via a left-handed or right-handed current, affects the energy distribution and angular correlation of the outgoing electrons in two-neutrino double beta decay. We pay particular attention on the behaviour of the energy distribution at the kinematic endpoint and we estimate the current limits on the active-sterile mixing and effective right-handed coupling using current experimental data as a function of the sterile neutrino mass. We also investigate the sensitivities of future experiments. Our results complement the corresponding constraints on sterile neutrinos from single beta decay measurements in the 0.1 -10 MeV mass range.

Introduction
Sterile neutrinos are among the most sought-after candidates of exotic particles. The main motivation for their existence is the fact that the Standard Model (SM) does not contain right-handed (RH) counterparts of the left-handed neutrino states participating in electroweak interactions, in contrast to the quarks and charged leptons. Their absence is purely because they are required to be singlets under the weak SU (2) L and have zero weak hypercharge if they are to participate in a Yukawa interaction with the left-handed neutrino states and the SM Higgs. They are thus truly sterile with respect to the SM gauge group and in this paradigm can only manifest themselves through an admixture with the active neutrinos. Thus, sterile neutrinos can in fact be considered to be any exotic fermion that is uncharged under the SM gauge interactions; unless protected by some new symmetry they will mix with the active neutrinos as described above. Hence, sterile neutrinos are also often referred to as heavy neutral leptons.
An important feature of their mixing with the active neutrinos is the resulting impact on the light neutrino masses. Neutrino oscillations [1] imply that the active neutrinos have small but non-zero masses. By adding a SM-singlet, RH neutrino field ν R per generation to the SM, neutrinos can become massive. A so-called Dirac mass term can be generated through the Yukawa interaction with the SM Higgs, though the coupling required is tiny with effects unmeasurable experimentally. In any case, the sterile states will be allowed to acquire a so-called Majorana mass, unless protected by lepton number conserving symmetry, modifying the spectrum and nature of neutrinos considerably. This of course refers to the well-known type I seesaw mechanism [2][3][4][5][6]. The sterile neutrinos were initially considered to be very heavy (m N ∼ 10 14 GeV) in order to generate the correct light neutrino masses, but there is now a strong theoretical and experimental incentive to consider sterile neutrinos at accessible energies. Fig. 1 summarises the current constraints on the active-sterile mixing strength |V eN | 2 in the regime 1 eV < m N < 10 TeV derived from numerous experiments. The most stringent limits from fixed target and collider experiments can be found in the mass range 1 GeV < m N < 100 GeV, a region motivated by leptogenesis models.
Lighter sterile neutrino masses, while challenging to accommodate due to the constraints from astrophysics, are still of interest, especially around m N ∼ 10 keV where sterile neutrinos may act as warm dark matter. In the regime 10 eV < m N < 1 MeV, nuclear beta decays are currently the only laboratory-based experimental method able to probe sterile neutrinos. Neutrinoless double beta (0νββ) decay is an exception, setting stringent limits over the whole range in Fig. 1, but only if the sterile neutrinos are Majorana fermions -for sterile Dirac neutrinos or quasi-Dirac neutrinos with relative splittings ∆m N /m N 10 −4 [7], the constraints vanish or become weaker respectively. In addition, if the sterile neutrinos in question are wholly responsible for giving mass to the active neutrinos via the type I seesaw mechanism, the contributions from active and sterile neutrinos to 0νββ decay cancel each other, see Sec. 3.1. Thus, especially around m N ∼ 1 MeV, the current constraints are rather weak, of the order |V eN | 2 few ×10 −3 . As mentioned above, the constraints arise from searches for kinks in the electron energy spectrum and measurements of the f t value of various beta decay isotopes, see Sec. 3.2 for a brief review. This weakening of limits motivates the use of novel methods to constrain the activesterile mixing in this mass regime. In this work we assess the potential of 0νββ decay experiments being sensitive to kinks in the background two-neutrino double beta (2νββ) decay spectrum caused by the presence of sterile neutrinos in the final state with masses m N 1 MeV. This is fully analogous to the corresponding searches in single beta decays but 2νββ decaying isotopes typically have Q values of a few MeV and are thus expected probe sterile neutrinos in such a mass range. The 2νββ decay process is of course very rare so it may at first seem difficult to achieve high enough statistics. While 2νββ decay is indeed not expected to improve the limits considerably, the 2νββ decays spectrum will be measured to high precision in several isotopes as 0νββ decay is searched for in ongoing and future experiments. The relevant data to look for sterile neutrinos in 2νββ decays will be available, which, generally speaking, can be used to look for signs of new physics in its own right [8,9].
In addition to a truly sterile neutrino, i.e. one that inherits the SM charged-current Fermi interaction albeit suppressed by the active-sterile mixing, we also consider RH current interactions of the 'sterile' neutrino, e.g. arising in left-right symmetric models. Such interactions change the angular distribution of the electrons emitted in 2νββ decay [8]. We parametrise all interactions in terms of effective operators of the SM with a light sterile neutrino, suitable in 2νββ decays with characteristic energies of 10 MeV. This paper is organised as follows. In Sec. 2 we introduce the effective operators relevant for our discussions. In Sec. 3 we briefly review the current limits on the activesterile mixing squared |V eN | 2 with a focus on the mass regime m N ∼ 1 MeV. The calculation of the 2νββ decay spectrum with the emission of one sterile neutrino is described in Sec. 4. Sec. 5 introduces our statistical procedure and presents the estimated current limits and prospective future sensitivities from sterile neutrino searches in 2νββ decay as our results. We conclude in Sec. 6.

Effective Interactions with Sterile Neutrinos
We consider the SM with the addition of a gauge singlet fermion N , i.e. the sterile neutrino. As we consider the second-order weak process of 2νββ decay, we restrict ourselves to the first generation of SM fermions. For processes with energies 100 GeV we can describe the relevant weak processes using the effective SM Fermi interaction. The sterile neutrino inherits the Fermi interaction, but is suppressed by the active-sterile mixing V eN . In addition, we allow the sterile neutrino to participate in exotic RH V + A interactions. The effective Lagrangian taking into account the above takes the form with the tree-level Fermi constant G F , the Cabbibo angle θ C , and the leptonic and hadronic respectively. The SM electroweak radiative corrections are encoded in δ SM . The active-sterile mixing is V eN and the XY encapsulate effects from integrating out new physics giving rise to V + A currents of the sterile neutrino. We neglect any further effective operators, such as exotic contributions to the SM Fermi interaction and RH currents with the active neutrino [8].
In Eq. (2.1), ν and N are 4-spinor fields of the light electron neutrino and the sterile neutrino. They are either defined to be Majorana fermions, ν = ν L +ν c L , N = N c R +N R (i.e. a Majorana spinor constructed from the left-handed Weyl spinor and its charge-conjugate) or Dirac fermions ν = ν L +ν R , N = N R +N L (a Dirac spinor constructed from two different Weyl fields). The calculation of 2νββ decay is not affected by this, i.e. it is insensitive to the Dirac versus Majorana character. If the neutrinos are Majorana the constraints from 0νββ decay must be considered.

Constraints on Sterile Neutrinos
In this section we review the constraints on the active-sterile mixing strength squared |V eN | 2 as a function of the sterile neutrino mass m N . We mainly concentrate on limits in the 0.1 MeV < m N < 3 MeV mass range. This is because 2νββ decay measurements are only sensitive to sterile neutrino masses below the Q value of the 2νββ decay process, which is of order Q ∼ 1 − 3 MeV for the isotopes of interest. The relevant constraints in this range come from the non-observation of 0νββ decay, single beta decay spectra, sterile neutrino decays and cosmological probes. We will see that the same constraints also apply broadly to the RH current couplings | LR | 2 and | RR | 2 .
As an overview we show in Fig. 1 the existing |V eN | 2 constraints over the mass range 1 eV < m N < 10 TeV; for further information on each labelled constraint see Sec. 4 of Ref. [10] and references therein. It is interesting to note the relative weakness of the upper limits from single beta decay experiments in the range 0.1 MeV < m N < 3 MeV. Mixing strengths are nonetheless excluded down to |V eN | 2 10 −7 −10 −6 and |V eN | 2 10 −14 −10 −11 by 0νββ decay and cosmological probes, respectively. It is crucial though to emphasise that the former constraints are model-dependent and can be avoided if neutrinos are Dirac fermions or if the sterile neutrinos are responsible for the light neutrino mass generation. Cosmological constraints rely on modelling of the early universe and can be avoided in extended scenarios where the sterile neutrinos have exotic interactions with a dark sector [11][12][13][14]. This therefore motivates looking at the sensitivities of current and future 2νββ decay measurements but we first look at the existing constraints within the region of interest in more detail.

Neutrinoless Double Beta Decay
If we consider the active and sterile neutrinos to be purely Dirac fermions, lepton number is conserved and 0νββ decay is forbidden. Searches for this decay will thus not provide constraints on the active-sterile mixing of Dirac neutrinos.
In the Majorana case, if n S sterile neutrinos are added to the SM with masses m N i and active-sterile mixing strengths V eN i (we assume for simplicity a single active state ν e ), the inverse of the half-life T 0ν 1/2 for the 0νββ decay process can be written using the interpolating formula (3.1) Here G 0ν is the phase space factor, g A is the axial vector coupling, M 0ν is the light neutrino exchange nuclear matrix element and p 2 is the average momentum transfer of the process [15,16]. By considering a single sterile neutrino with mass m N and neglecting the contribution from the active neutrinos, the constraint in Fig. 1 is derived using the current experimental bounds.
If the heavy states are related to the light state by a seesaw relation, then must be satisfied. Thus, if the sterile states are lighter than the 0νββ decay momentum transfer, m N i p 2 , the 0νββ decay rate vanishes and the corresponding constraint in Fig. 1 disappears. Sterile neutrinos have been discussed in the context of 0νββ decay in detail in Refs. [10,17,18].

Beta Decay
Electron neutrinos are produced in the beta decays of unstable isotopes via the LH chargedcurrent interaction. If the active-sterile mixing strength |V eN | 2 or RH couplings | LR | 2 , | RR | 2 are non-zero, sterile neutrinos can be produced if their masses are smaller than the Q value of the process. For a large enough m N the emission results in a distortion or 'kink' in the beta decay spectrum and associated Kurie plot.
The beta decay spectrum with respect to the kinetic energy of the emitted electron can be written for a single sterile neutrino with mixing as the incoherent sum where we neglect the light neutrino masses in the standard contribution. Due to unitarity, the contribution from the light neutrinos is reduced by the active-sterile mixing strength.
The sterile neutrino contribution gives rise to a kink in the spectrum of relative size |V eN | 2 and at electron energies E e = Q − m N . Alternatively, in the case the sterile neutrinos are produced by a RH current, the SM contribution is no longer reduced as a result of unitarity.
This weakens the upper limits on | LR | 2 and | RR | 2 compared to |V eN | 2 , though the effect is negligible for upper bounds below 10 −2 . Kink searches have been conducted for a variety of isotopes with different Q values, making them sensitive to a range of sterile neutrino masses. Shown in Fig. 1 [28] and 187 Re [29], assuming there to be a single sterile state. With smaller Q values, 3 H and 187 Re provide constraints over the range 1 eV < m N < 1 keV. It can be seen that 45 Ca, 64 Cu, 144 Ce-144 Pr and 20 F in the mass range of interest provide slightly weaker upper bounds (between 10 −3 and 10 −2 ) compared to 63 Ni and 35 S at lower masses.

Sterile Neutrino Decays
A sterile neutrino produced in the beta decay of a neutron-rich isotope in a reactor or a light element in the sun can decay before detection via the channels N → ννν and N → e + e − ν. The former channel is mediated by a neutral current and the latter via either a neutral or charged current. The latter also requires the sterile neutrino mass to be m N > 2m e . At tree-level (in the single-generation case) the total decay rate is given approximately by where the factor of 2 is present in the Majorana case. For RH currents, the factor |V eN | 2 is replaced by Reactor experiments with neutrino energies ∼ 10 MeV are sensitive to sterile neutrinos with masses in the range 1 MeV < m N < 10 MeV. Limits have been set by searches at the Rovno [30] and Bugey [31] reactors. Sterile neutrino decays were also searched for by the Borexino experiment [32] which was sensitive to heavy neutrinos with masses up to 14 MeV produced in the decays of solar 8 B nuclei. Borexino enforces the relatively stringent limit |V eN | 2 10 −6 − 10 −5 for m N ∼ 10 MeV.

Cosmological and Astrophysical Constraints
The presence of sterile states with mixing strengths |V N | 2 (and/or the presence of RH currents) has wide-ranging consequences for early-universe observables. These include the abundances of light nuclei formed during Big Bang Nucleosynthesis (BBN), temperature anisotropies in the Cosmic Microwave Background (CMB) radiation and the large-scale clustering of galaxies [33]. Deviations from the standard smooth, isotropic background evolution (and perturbations around this background) impose severe constraints -the region between the grey lines labelled CMB+BAO+H 0 (an upper limit) and BBN (a lower limit) is excluded. These limits are highly sensitive however to the production and decay mechanism of the sterile state and can be relaxed in certain models.
The main constraint to consider in the 0.1 MeV < m N < 3 MeV mass range is the upper limit labelled CMB + BAO + H 0 . Via the active-sterile mixing or RH current, sterile states are populated in the early-universe and they decouple when the Hubble expansion overcomes the interaction rate with the SM particles. It is then possible for these states to decay at later times to produce non-thermally distributed active neutrinos, modifying the amount of extra radiation measured at recombination, ∆N eff , beyond the usual value including active neutrino oscillations, N eff 3.046. Useful probes include the CMB shift parameter R CMB , the first peak of the Baryon Acoustic Oscillations (BAO) and the Hubble parameter H(z) inferred from type Ia supernovae, BAO and Lyman-α data. These exclude values of m N and |V eN | 2 corresponding to lifetimes up to the present day, where the condition that N does not make up more than the observed matter density Ω sterile < Ω DM ≈ 0.12 h −2 also applies. This constraint can be evaded in exotic models [11][12][13][14], for example those that inject additional entropy and dilute the dark matter (DM) energy density.

Double Beta Decay Rate with a Sterile Neutrino
Considering one sterile neutrino N with mass m N < Q ββ few MeV and a SM chargedcurrent as in Eq. (2.1) with additional suppression by the active-sterile mixing strength V eN allows for the possibility that in 2νββ decay oneN is emitted (νN ββ) instead of ā ν e (we assume that N is long-lived and does not decay within the detector, thus being invisible). The final state is different from the standard 2νββ decay and thus there is no interference between νN ββ and 2νββ. There is also no anti-symmetrisation with respect to the two different neutrinos in νN ββ. Moreover, a RH lepton current can be also assumed to be associated with the emission of the sterile neutrino, which further affects the 2νββ observables, mainly the angular correlation of the outgoing electrons.
In order to write down expressions for the 2νββ and νN ββ decay rates, including the possibility of RH currents, let us start with the general expression [34] denote the energies of initial and final nuclei, electrons and antineutrinos, respectively. The magnitudes of the associated spatial momenta are p e i = |p e i | and pν i = |pν i | and m e and m ν denote the electron and neutrino masses. The phase space differentials are dΩ e 1 = d 3 p e 1 /(2π) 3 , etc. Note that in our calculations we neglect the mass of the light neutrino being emitted and we retain only the mass m N of the heavy neutrino.
After integrating over the phase space of the outgoing neutrinos, the resulting differential 2νββ decay rate can be generally written in terms of the energies 0 ≤ E e 1 , E e 2 ≤ Q+m e of the two outgoing electrons, with Q = E i − E f − 2m e , and the angle 0 ≤ θ ≤ π between the electron momenta p e 1 and p e 2 as [34] where with G β = G F cos θ C (G F is the Fermi constant and θ C is the Cabbibo angle).
The quantities A 2ν and B 2ν in Eq. (4.2), generally functions of the electron energies, include the integration over the neutrino phase space, where we have used Eν 1 due to energy conservation and kept the dependence on the neutrino masses, although in the SM case they can be safely neglected. In turn, the quantities A 2ν and B 2ν , generally functions of the electron and neutrino energies, are calculated below using the nuclear and leptonic matrix elements.
The rate corresponding to νN ββ decay then differs only by the non-negligible mass of the sterile neutrino entering the neutrino energy and, most importantly, the integration bounds. Consequently, the corresponding rate can be obtained from the above by a simple substitution ν 1 → N , ν 2 → ν and neglecting the mass m ν . As shown later in this section, in the standard case with only LH lepton currents the quantities A 2ν and B 2ν do not depend on neutrino masses; hence, the main effect of the sterile neutrino mass is the shrunk electron energy distribution given by the effectively smaller Q value, now given by Q = In our calculations we take the S 1/2 spherical wave approximation for the outgoing electrons, i.e.
Here,p e = p e /|p e | denotes the direction of the electron momentum, χ s is a two-component spinor and g −1 (E e ) and f +1 (E e ) stand for the radial electron wave functions depending on the electron energy E e . As commonly done, we approximate them with their values at the nucleus' surface, i.e. at distance R from the centre of the nucleus. The neutrinos, being neutral, can be simply described as plane waves in the long-wave approximation, Eν +mν χ s . (4.7)

Purely Left-Handed Currents
The standard contribution to 2νββ decay given by the first term in the Lagrangian in Eq. (2.1) has been studied in great detail [35,36]. Sticking to the formalism outlined above, the decay rate is described by the functions and where we define Fermi and Gamow-Teller nuclear matrix elements (4.10) The electron mass m e in the above expression is inserted conventionally to make the nuclear matrix elements dimensionless. The lepton energies enter in Eq. (4.10) through the terms which satisfy −Q/2 ≤ ε K,L ≤ Q/2. In case of 2νββ decay with energetically forbidden transitions to the intermediate states, E n − E i > −m e , the quantity E n − (E i + E f )/2 = Q/2 + m e + (E n − E i ) is always larger than Q/2. The above expressions may be further simplified using several well-motivated approximations.
Isospin Invariance: Neglecting the isospin non-conservation in the nucleus, the double Fermi nuclear matrix elements vanish, i.e. M K F = M L F = 0. Therefore, Eqs. (4.12) and (4.13) then respectively acquire the approximate form and

.13)
Nuclear matrix element dependence on lepton energies: If we neglect the dependence of nuclear matrix elements on ε K,L , the nuclear and leptonic parts can be separated and we get 14) with the Gamow-Teller nuclear matrix element now defined as A better approximation is obtained by Taylor expansion of the nuclear matrix elements in the small parameters K,L [36]. Keeping terms up to the fourth power in K,L gives and Here, the nuclear matrix elements introduced are defined as This is the approximation we employ in our later numerical analyses.

Contribution with a Right-Handed Current
The non-standard contribution to 2νββ decay involving the RH currents proportional to the XR coupling, as appearing in the Lagrangian in Eq. (2.1), was calculated in Ref. [8].
The corresponding functions A 2ν and B 2ν entering Eq. (4.2) read Here, the dependence on the electron radial wave functions has been made explicit. Likewise, the terms proportional top 1 ·p 2 = cos θ combine to give In Eqs. (4.22) and (4.23), the terms proportional to m ν m N are small, as one of the emitted neutrinos is still assumed to be the light with m ν 0.1 eV. As in the SM case, for the purpose of numerical computations we approximate the above expressions with their Taylor expansions up to the fourth power in the small parameters K,L .

Decay Distributions and Total Rate
The kinematics of the electrons emitted in the decay is captured by the fully differential decay rate expressed in Eq. (4.2) depending on the (in principle) observable electron energies E e 1 , E e 2 and the angle θ between the electron momenta. All the information is contained  Table 1: Nuclear matrix elements calculated within the pn-QRPA with partial isospin restoration [36] assuming the effective axial coupling g A = 1.0.
by the quantities A 2ν and B 2ν presented above both for the standard LH (Eqs. Since quenching of the axial coupling g A is expected in the nucleus [37], we take g A = 1 instead of the usual value g nucleon A = 1.269 for a free neutron. Further, we use the 2νββ decay nuclear matrix elements from Ref. [36], as shown in Tab. 1.
With all the above ingredients we can now calculate the the total decay rate as well as various decay distributions potentially observable in 2νββ decay experiments.
Total electron energy and single electron energy: The 2νββ decay experiments measure primarily the distribution with respect to the total kinetic energy of the outgoing electrons, i.e. dΓ 2ν /dE K with E K = E e 1 + E e 2 − 2m e − m ν 1 − m ν 2 . Here, m ν 1,2 denote the masses of the emitted neutrinos, which can be safely neglected in the SM case, but we consider also a contribution involving a heavy sterile neutrino, in which case one of the masses becomes non-negligible and we denote it m N . We also neglect the recoil of the final state isotope which would change the endpoint by ∼ Q 2 /M 0.1 keV, with the Q 3 MeV and the mass of the nucleus M ≈ 76 − 136 GeV. Some experiments capture the energies and tracks of individual electrons, thus allowing for study of the single electron energy distribution dΓ 2ν /dE e 1 (the symmetry of the process ensures the distribution with respect to the second electron is identical) and the double differential distribution dΓ 2ν /(dE e 1 dE e 2 ). These distributions are calculated from Eq. (4.2) as where in the latter and E max We neglect the light neutrino masses m ν 1 = m ν 2 = 0 in the SM case and retain only the heavy neutrino mass in the sterile contribution, m ν 1 = m N , m ν 2 = 0. Given the fact that most experiments provide only the 2νββ decay distribution in dependence on the total kinetic energy of the electrons, in the following analysis we focus primarily on this observable.
Unlike in the case of single beta decay, the electron energy spectrum of 2νββ decay with the emission of a sterile neutrino does not manifest the typical square-root behaviour towards its endpoint. Hence, no kinks are expected to appear in the total spectrum including both the SM and the sterile neutrino contributions. As shown in the following section, the total electron energy distribution of 2νββ decay is rather smooth, behaving approximately as near the kinematic endpoint E K Q − m N .
Angular correlation factor and total decay rate: The integration over the electron energies leads to the equation describing the angular distribution of the decay. Here, Γ 2ν denotes the total 2νββ decay rate and K 2ν = Λ 2ν /Γ 2ν stands for the angular correlation factor, which are given by As the inclusion of RH current leads to the opposite sign of the angular correlation of the emitted electrons [8], it can be also used to distinguish the corresponding contributions, as analysed in the following section.

Constraints on Sterile Neutrino Parameters
We will now use the differential 2νββ decay rates derived in Sec. 4 to exclude regions of the sterile neutrino parameter space -namely, the sterile neutrino mass m N and mixing with the electron neutrino |V eN | 2 . To do this we will first outline a simple frequentist limit setting method. We will then use the non-observation of deviations from the SM 2νββ decay spectrum by 0νββ decay search experiments such as GERDA-II, CUPID-0, NEMO-3 and KamLAND-Zen to put upper limits on |V eN | 2 as a function of m N . We will also estimate upper limits from the forecasted sensitivities of future 0νββ decay experiments such as LEGEND, SuperNEMO, CUPID and DARWIN. Finally, we will compare these upper limits to existing constraints in the 0.

Statistical Procedure
To obtain upper limits on the mixing |V eN | 2 we follow the standard frequentist approach of Refs. [1,38]. Firstly, we define the total differential 2νββ decay rate as the incoherent sum of the sterile neutrino and SM rates for a given sterile mass m N and total kinetic energy explicitly writing the dependence on active-sterile mixing |V eN | 2 . The total differential rate depends on the sterile neutrino parameters ξ ≡ (m N , |V eN | 2 ) and E K .
In Fig. 2, the total differential decay rate in Eq. (5.1) is compared to the sterile neutrino contribution |V eN | 2 ·dΓ 2ν N /dE K (where both are normalised to the total SM decay rate Γ 2ν SM ) for the isotopes 100 Mo and 136 Xe. The respective Q values of the isotopes are indicated by the vertical dotted lines and the values m N = 1.0 MeV and |V eN | 2 = 0.5 are chosen. In the panel below we show the corresponding percentage deviation of the total differential rate from the SM rate, It can be seen that the magnitude of dΓ 2ν /dE K decreases with respect to the dΓ 2ν SM /dE K as the total kinetic energy increases, eventually plateauing at around −10%. This is because the sterile neutrino contribution |V eN | 2 dΓ 2ν N /dE K falls as E K increases above ∼ 1.0 MeV. Eventually its contribution is negligible, but there remains a suppression from the (1 − |V eN | 2 ) factor multiplying the SM contribution, which is particularly sizeable for the choice |V eN | 2 = 0.5. It is apparent from Eq. (5.2) that the deviation tends to a factor of −|V eN | 2 . The characteristic signature of the sterile neutrino is a relative increase of the differential rate for E K Q − m N .
Any experiment measuring the 2νββ decay spectrum will count a number of events N events distributed over a number of bins N bins in the total kinetic energy E K . In the presence of a sterile neutrino, the expected fraction of events ∆N (i) exp per bin will be the integral of dΓ 2ν /dE K over the width of the bin from the total kinetic energy E i to E i+1 , where the normalisation factor N is i.e. the total area enclosed by dΓ 2ν /dE K between kinetic energies E min and E max . The total number of expected events per bin will then be where we have also split the expected number of events into the number of signal and background events as The probability of the experiment observing N where the second equality holds via Wilks' theorem if there are a large number of events per bin [39]. From this we can construct the test-statistic q ξ = −2 log L(D|ξ) − log L(D|ξ) , (5.9) whereξ are the values of the sterile neutrino parameters that minimise the log-likelihood function. The quantity q ξ is expected to follow a χ 2 distribution with one degree of freedom. We assume that the experiment does not observe a spectrum deviating significantly from the SM prediction. We therefore set the number of observed events in Eq. (5.8) to N (i) N , 0). In reality, however, the experiment could be repeated many times and record a different value of N (i) obs each iteration. This fluctuation can be imitated by running a series of toy Monte Carlo simulations of the experiment. For every toy Monte Carlo there is a value of q ξ , with the relevant test-statistic becoming the median of these values. A representative data set is commonly used as a good approximation of the MC method in the large sample limit [40]. This is the so-called Asimov data set D A for which the observed number of events per bin N (i) obs equals the number of background events N (i) bkg [41]. Theξ that minimises the log-likelihood to −2 log L(D A |ξ) = 0 is then simplŷ ξ = (m N , 0) which matches our initial approach.
The magnitude of the test-statistic q ξ = −2 log L(D A |ξ) translates to a degree of compatibility between the Asimov data set and the sterile neutrino hypothesis with parameters ξ = (m N , |V eN | 2 ). For example, if both parameters are allowed to vary, combinations of the parameters giving q ξ 4.61 are excluded at 90% confidence level (CL). Rather than performing a two-dimensional scan of the parameters, we instead fix m N for values over the range ∼ 0.1 − 3 MeV and find the value of |V eN | 2 for which q ξ = 2.71, corresponding to the 90% CL upper limit on the mixing.
Finally we note that we have not yet included the effect of systematic uncertainties. Systematics altering the total number of observed events without leading to distortions in the spectrum can be accounted for by introducing the nuisance parameter η where σ η is a small associated uncertainty. The remaining systematic uncertainties are included in the quantity σ which will be used derive constraints in the next subsection. The same procedure can be applied to place upper limits on the RH current couplings | LR | 2 and | RR | 2 . As seen in the previous sections, the RH current modifies the total kinetic energy distribution to where the SM contribution is no longer reduced by the sterile neutrino mixing. The RH current also modifies the angular distribution to Eq. (4.27) with the total rate Γ 2ν and the angular correlation factor K 2ν given in terms of SM and RH current contributions as (5.14) Assuming | XR | 2 1, K 2ν can be Taylor expanded as where the SM contribution and RH current contributions, respectively, are The SM values are K 2ν SM = −0.627 for 100 Mo and K 2ν SM = −0.631 for 82 Se (the isotopes of experiments that are sensitive to the angular correlation factor, NEMO-3 and SuperNEMO, respectively). The α(m N ) factors are plotted for 82 Se (red) and 100 Mo (blue) in Fig. 3, which also indicates the values at m N = 0. The factor α(m N ) is positive, indicating a change of the angular distribution away from the back-to-back configuration of electrons in the SM V − A case. It is maximal for m N = 0 and is suppressed to zero as m N approaches the Q value.
Using the measured total kinetic energy distributions from all 2νββ decay experiments, the ξ = (m N , | XR | 2 ) parameter space can be constrained in the same was as (m N , |V eN | 2 ) described above, i.e. using the test-statistic in Eq. (5.12). In addition, the experiments NEMO-3 and SuperNEMO will measure a certain number of events N  (2-5) × 10 4 10 6 -10 7 5 (0.5, 0.5)   KamLAND-Zen ( 136 Xe, blue). We also show a combined constraint (black dashed) found by summing the log-likelihoods of the experiments (each minimised with respect to a separate nuisance parameter η). It can be seen that the upper limits worsen for smaller and larger values of the sterile mass in the range 0.1 MeV < m N < 3 MeV, with the most stringent upper bound being found at m N similar to the peak energy of the associated spectrum. The constraints are compared to pre-existing constraints (shaded areas) from single beta decay experiments and sterile neutrino decays. While NEMO-3 and KamLAND-Zen provide the best individual constraints (|V eN | 2 0.02), they are not as competitive as previous limits. However, it is promising that 2νββ decay is more sensitive for sterile masses 0.3 MeV < m N < 0.7 MeV where existing constraints are less stringent. Fig. 4 (right) shows the corresponding sensitivities estimated for the next generation of 0νββ decay experiments. The forecasted range of exposures given by the collaborations are often one or two orders of magnitude larger than those of the current generation. We estimate the total number of events N events seen in future by multiplying the current values by the ratio of future to current exposures. Energy resolutions are taken from the references in Tab. 2 and we assume an optimistic value of σ η ∼ σ f ∼ 0.5% for the systematic uncertainties. We compute the 90% CL sensitivity for both the higher and lower forecasted number of events in Tab  sensitivity (black dashed) using the largest predicted exposure of each experiment. For a given experiment the upper bounds exhibit the same improvement for sterile masses close to the maximum of the total differential decay rate. The most stringent upper limits come from CUPID and DARWIN, |V eN | 2 2.5 × 10 −3 , which would exclude the currently unconstrained region in the 0.3 MeV < m N < 0.7 MeV range.

A selection of current and next generation
Likewise, we estimate the current limits and future sensitivity on the RH couplings | LR | 2 and | RR | 2 from measuring the 2νββ decay energy distribution and angular correlation. In Fig. 5 we plot the upper limits at 90% CL on | LR | 2 and | RR | 2 as a function of the sterile neutrino mass m N . The blue solid line is the combined constraint from current 2νββ decay experiments using the total kinetic energy distribution, while the red solid line is the upper limit derived from the angular distribution measurement of NEMO-3 ( 100 Mo). The blue dashed line is the combined sensitivity from future 2νββ decay experiments, while the red dashed band indicates the sensitivity range from the angular distribution measurement of SuperNEMO ( 82 Se). The latter does not improve over the current limit as SuperNEMO is not expected to have a significantly increased exposure compared to NEMO-3, see Tab. 2. We therefore also indicate the sensitivity of a hypothetical 82 Se angular measurement with an exposure of 10 7 events (red dot dashed).
Due to the different total kinetic energy distribution for the RH current in Eq. (5.13) (no suppression of the SM rate), the combined constraints on | LR | 2 and | RR | 2 (dashed lines) are slightly weaker than the equivalent constraints on |V eN | 2 . The constraints from the NEMO-3 angular distribution are generally better, tending to a constant upper bound | XR | 2 10 −3 for m N 0.2 MeV. This roughly agrees with the result XR < 2.7 × 10 −2 in the massless case found in Ref. [8].

Conclusions
Measuring the kinematic endpoint in single beta decay is arguably the cleanest means to determine the absolute neutrino masses in a model-independent fashion. For the light active neutrinos in the SM, the most promising isotope for this is tritium ( 3 H) and its beta decay is currently measured in the KATRIN experiment [51] as well as the future Project 8 [52] and CRESDA [53] efforts. The same method can be applied to search for sterile neutrinos, not only in Tritium but in a host of beta decay isotopes where masses smaller than the respective Q value of the decay can be probed. The limits on the active-sterile mixing strength |V eN | 2 from such searches are summarised in Fig. 6. They are comparatively weak, of the order |V eN | 2 × 10 −2 − 2 × 10 −3 , in the sterile neutrino mass range 0.1 MeV < m N < 1 MeV.
In this work, we have analysed the prospects to search for sterile neutrinos using the same principle in 2νββ decay. If one of the two neutrinos emitted in the process is a heavier, sterile neutrino it will likewise affect the distribution with respect to the kinetic energy of the two electrons observed in the decay: the kinematic endpoint is shifted to lower values depending on the sterile neutrino mass and the active-sterile mixing will reduce the usual SM contribution. This is expected to be challenging because of the very long 2νββ decay half lives and small rates compared to single beta decay. Nevertheless, future searches for the lepton number violating 0νββ decay will push the envelope in terms of exposure and allow measuring 2νββ decay with up to 10 7 events. This data can then be used to probe exotic physics with 2νββ decay in its own right. Apart from sterile neutrino searches, other examples include exotic neutrino self interactions [9] and RH leptonic currents [8]. We have extended the latter analysis here to consider a RH V + A current for a sterile neutrino rather than the SM electron neutrino. As in Ref. [8], this gives rise to an anomalous angular distribution of the electrons in 2νββ decay.
To summarise the sensitivity we compare in Fig. 6 the current limits on |V eN | 2 from existing 2νββ decay (solid blue) to constraints from single beta decays and sterile neutrino decays over a wider range of masses, 100 eV < m N < 10 MeV. The blue curve uses the combined constraints from measurements of 2νββ decay electron energies. The red curve shows the current constraint on the effective RH coupling | XR | 2 using the NEMO-3 angular distribution measurement. The dashed curves indicate the corresponding future sensitivities. At lower masses both the current and future upper limits on |V eN | 2 cannot compete with existing constraints from 64 Cu and 144 Ce− 144 Pr beta decays. At higher masses they are also less stringent than constraints from Borexino, Bugey and Rovno. It is the 0.3 MeV < m N < 0.7 MeV range where 2νββ decay can provide competitive constraints  | XR | 2 (2νββ Angular) Figure 6: Current upper limits (solid blue) and future sensitivities (dashed blue) on the mixing strength |V eN | 2 between the electron and sterile neutrino as a function of the sterile mass m N . Likewise, the red curves give the current limit and future sensitivity on the RH coupling | XR | 2 using a measurement of the angular distribution in 2νββ decay. The shaded regions are excluded by existing searches in single beta decay and sterile decays in reactor and solar neutrino oscillation experiments.
in the future, though we expect that similar improvements from 20 F and 144 Ce− 144 Pr beta decays are also possible. The constraints on the RH coupling | XR | 2 using an angular distribution measurement in 2νββ decay is most sensitive for light sterile neutrino masses m N 0.1 MeV as the effect is phase space suppressed otherwise. We note, though, that the limits from single beta decays and the other processes shown strictly speaking apply to |V eN | 2 only and need to be re-evaluated for a heavy neutrino coupling through a RH current.
Our analysis demonstrates that 2νββ decay can be used to search for sterile neutrinos with masses lighter than m N ∼ 1 MeV. While current searches are not competitive with limits from single beta decays, future searches will have a much more increased statistics where effects of new physics can be tested. While sterile neutrinos in this mass range are also heavily constrained from astrophysical measurements and cosmological considerations, it is important to improve our understanding using all available data.