Search for Light Exotic Fermions in Double-Beta Decays

The Standard Model of Particle Physics predicts the double-$\beta$ decay of certain nuclei with the emission of two active neutrinos. In this letter, we argue that double-$\beta$ decay experiments could be used to probe models with light exotic fermions through the search for spectral distortions in the electron spectrum with respect to the Standard Model expectations. We consider two concrete examples: models with light sterile neutrinos, singly produced in the double-$\beta$ decay, and models with a light $Z_2$-odd fermion, pair produced due to a $Z_2$ symmetry. We estimate the discovery potential of a selection of double-$\beta$ decay experiment and find that future searches will test for the first time a new part of the parameter space of interest at the MeV-mass scale.

The Standard Model of Particle Physics predicts the double-β decay of certain nuclei with the emission of two active neutrinos. In this letter, we argue that double-β decay experiments could be used to probe models with light exotic fermions through the search for spectral distortions in the electron spectrum with respect to the Standard Model expectations. We consider two concrete examples: models with light sterile neutrinos, singly produced in the double-β decay, and models with a light Z2-odd fermion, pair produced due to a Z2 symmetry. We estimate the discovery potential of a selection of double-β decay experiment and find that future searches will test for the first time a new part of the parameter space of interest at the MeV-mass scale.

I. INTRODUCTION
Many models of new physics contain new spin 1/2 particles, singlets under the Standard Model gauge group, possibly related to the mechanism of neutrino mass generation, or to the dark matter of the Universe. An archetype of such fermions is the sterile neutrino, also called right-handed neutrino. In a variant of this scenario, the singlet fermion is furnished with a Z 2 symmetry so it can only be produced in pairs. Unfortunately, in either of these scenarios, the mass of the light exotic fermion and the coupling strength to the Standard Model particles are free parameters of the model, leaving a vast parameter space that must be probed in laboratory experiments, or with astrophysical and cosmological observations. From the phenomenological standpoint, sub-GeV sterile neutrinos are particularly attractive, since they could be produced in charged lepton or in hadron decays and possibly be discovered in a laboratory experiment. In fact, there is an intense experimental program aiming to detect signals of sterile neutrinos in the sub-GeV range, either in dedicated experiments, or as a by-product of an experiment initially designed for a different purpose. The current limits are fairly stringent below ∼ 100 keV and above ∼ 100 MeV (for an updated summary of the constraints, see e.g. [1]). However, the intermediate mass region remains relatively unexplored and some of the current constraints date back to the 1990s [2,3]. The scenario with Z 2 -odd singlet fermions has been even less studied, and, to the best of our knowledge, there is no laboratory constraint for the mass range ∼ 100 keV −100 MeV.
In this letter, we explore the possibility of searching for the production of light exotic fermions in double-β decays. These decays are nuclear transitions in which the atomic number Z increases by two units while the nucleon number A remains constant. The process is allowed by the Standard Model as long as two neutrinos and two electrons are emitted: (A, Z) → (A, Z + 2) + 2e + 2ν (2νββ decay). If light exotic fermions exist, they will also be produced in double-β decays, replacing one or both neutrinos in the final state. While neutrinos or exotic fermions are challenging to detect, the energy distribution of the electrons emitted in the process can be accurately measured. The shape of this distribution carries information about the full kinematic of the process and can be used to reconstruct which other particles have been emitted. The fermion mass range accessible with this kind of search extends from a few hundred keV to a few MeV, i.e., from the energy threshold of the experiments till the Q-values of the decay.
In the next couple of decades, several experiments with the capability of measuring the electron energy distribution in 2νββ decays will get online [4]. The primary physics goal of these experiments is the discovery of double-β decays with only the two electrons in the final state: (A, Z) → (A, Z + 2) + 2e. Observing this so-called neutrinoless double-β (0νββ) decay would have far-reaching implications. It would prove that the lepton number is not a global symmetry and that neutrinos have a Majorana mass component [5,6]. This physics goal already justifies major investments for experimental infrastructure. Nevertheless, there are several other valuable physics searches that can benefit from the low-background and large target mass of these experiments [7][8][9][10]. Our proposed search for light exotic fermions provides a valuable extension of their physics program.
We will discuss and estimate the discovery potential of double-β decay experiments under two illustrative scenarios that extend the Standard Model by adding light exotic fermions. In the first one, we assume the existence of a massive sterile neutrino N that mixes with the electron neutrino. This opens the possibility to observe a double-β decay final state in which a standard neutrino is replaced by a sterile neutrino (νN ββ decay): (A, Z) → (A, Z + 2) + 2e + ν + N , with a modified spectrum due to the sterile neutrino mass. This is indeed the same principle as the one used in kink searches with single-β decays [11][12][13][14][15], which currently provide the strongest laboratory bounds in our mass range of interest [2,3,16,17]. As we will see, the presence of light exotic fermions does not create a kink in the total double-β energy spectrum but still creates a continuous distortion that is detectable by 0νββ decay experiments.
The second scenario that we will consider is characterized by the presence of a new symmetry that forbids the production of a single exotic fermion in double-β decays. This is a typical scenario for models adding new exotic fermions χ as dark matter candidates, which are charged under a Z 2 symmetry to make them stable. This kind of models cannot be tested through single-β decays, but would result in a new final state in double-β decays (χχββ decay). Thus, double-β decay experiments can provide the first laboratory-based constraints on several models.
In the following, we will review the phenomenology of the aforementioned scenarios and present our calculations of how the electron energy distribution is affected by the model parameters. The experimental identification of distortions in the electron energy distribution requires a spectral fit whose accuracy is limited by statistical and systematic uncertainties. We will discuss how to set up a statistical analysis and evaluate the impact of systematic uncertainties. Finally, we will derive sensitivity projections for the most promising 0νββ experiments, assuming the systematic uncertainties and experimental performance achieved by running experiments.

II. DOUBLE-β β β DECAY INTO STERILE NEUTRINOS
One of the simplest extensions of the Standard Model consists in adding to the particle content one spin 1/2 particle, singlet under the gauge group, N . The gauge symmetry allows a Yukawa coupling of N to the Standard Model Higgs doublet and lepton doublets, which leads upon electroweak symmetry breaking to a Dirac mass term which couples N to ν, m D νN . For this reason, N is commonly denominated right-handed neutrino or sterile neutrino. Furthermore, the gauge symmetry also allows a Majorana mass term for N , M N c N which we assume to be in the range 0.1-2 MeV. The two parameters of the model, m D and M N are usually recast in terms of the misalignment angle between the interaction and mass eigenstates, sin θ m D /M , and the mass of the heaviest eigenstate m N M .
If kinematically possible, sterile neutrinos could be produced in any decay process involving Standard Model neutrinos, due to the active-sterile mixing angle. The new decay channels lead to distortions in the energy distribution for the visible particles compared to the Standard Model expectations, and which then constitute a test for the sterile neutrino scenario. In this work we concentrate on the double-β decay. In the presence of sterile neutrinos, and provided their mass is below the Q-value of the decay, the decay channels νN ββ and N N ββ become kinematically possible. The differential energy spectrum Γ is given by the incoherent superposition of three channels, and can be expressed as: where T = (E e1 + E e2 − 2m e )/m e is the sum of the kinetic energies of the two electrons (normalized to the electron mass). This variable is kinematically restricted to be 0 ≤ T ≤ T 0 , T 0 − x N and T 0 − 2x N , for 2νββ, νN ββ and N N ββ respectively, with x N = m N /m e and T 0 = Q ββ /m e . Q ββ is the end-point energy assuming massless neutrinos and depends on the particular nucleus. The 2νββ decay has been studied in several works [18][19][20][21][22][23][24][25], always assuming vanishing neutrino masses. In this work we extend this analysis, leaving the masses of the invisible fermions as free parameters. To this end, we follow and generalize [23], and express the differential rate for a generic decay (A, Z) → (A, Z + 2) + 2e + a + b, with a and b being either ν or N , as: with M eff 2ν is the dimensionless nuclear matrix element (NME). The phase space factor is given by where with E I the energy of the initial nucleus, E N a suitable excitation energy in the intermediate nucleus and Here, E X and p X = E 2 X − m 2 X denote the energy and the modulus of the 3-momentum of the particle X = e 1 , e 2 , a, b with mass m X , and subject to the condition which for m a , m b = 0 shifts the end-point of the spectrum to lower values. Finally, the factor f (0) 11 originates from the Coulomb interaction of the electrons with the daughter nucleus, which we parameterize using the Fermi function [20] f (0) 11 In the limit a, b = ν and m ν = 0 we recover the results of [23]. The complete analysis of the spectrum requires a numerical evaluation of Eq. (2). However, we also identified accurate analytical approximations for the dominant 2νββ and νN ββ channels (N N ββ is negligible for small sin θ), which assumes the Primakoff-Rosen approximation [26] for the Fermi function and neglects the lepton energies in Then, neglecting also the active neutrino mass but keeping the sterile neutrino one, we obtain: where we have introduced the form factors The left panel of Fig. 1 displays the phase space resulting from our numerical calculation applied to 76 Ge (Q ββ 2039 keV), showing how the presence of a final massive invisible fermion modifies the energy spectrum, shifting the end-point as well as the peak to lower values as m N increases. We can also see that the analytical expression reproduces the full result to a good approximation, up to a global factor of 1.67. This factor is due to the approximation we considered, but it does not affect the analysis, which is sensitive to the G νN /G νν ratio. We stress that, contrary to the single-β decay spectrum, the end-point of the νN ββ spectrum is smooth and therefore the total spectrum does not manifest a kink at Q ββ − m N . Nevertheless, the spectrum differs from the standard one, allowing 0νββ experiments to probe this hypothesis by measuring the 2νββ spectrum very accurately.
In general, the production of a pair of exotic particles is strongly suppressed compared to the production of a single one and can be neglected. However, there are scenarios in which the single production is forbidden and only the double production can take place. These scenarios cannot be tested by single-β decay experiments, whereas 0νββ decay experiments are still sensitive to the double production and have the unique opportunity to explore these channels.
In this section we consider a variant of the previous model, in which the symmetry group is extended by a discrete Z 2 symmetry, possibly related to the dark matter sector, which is exact or mildly broken in the electroweak vacuum. We assume that all Standard Model particles are even under the Z 2 symmetry, while the neutral singlet fermion is odd. We will denote the Z 2 -odd neutral singlet fermion as χ to differentiate it from the sterile neutrino, since the Z 2 symmetry forbids the mass term m D νN characteristic of the latter. Correspondingly, the Z 2 symmetry forbids the single-β decay (A, Z) → (A, Z + 1) + e + χ, as well as the double-β decay (A, Z) → (A, Z + 2) + 2e + ν + χ. However, the double-β decay (A, Z) → (A, Z + 2) + 2e + 2χ is possible (χχββ decay). The search for spectral distortions in the double-β decay spectrum can therefore probe scenarios with light Z 2 -odd singlet fermions. The total differential decay rate receives in this scenario two contributions: where dΓ νν /dT is the Standard Model contribution, defined in Eq. (2), and dΓ χχ /dT is the exotic contribution.
Here, T and T 0 are defined as in Eq. (1), with x χ = m χ /m e . The exotic contribution is very model dependent. For definiteness, we will consider the following effective interaction between the active neutrinos and χ: with a constant coupling g χ . We use the results of [10], that considered a similar four-fermion scalar interaction as in Eq. (14) but for neutrino self-interactions, to relate the decay rate for χχββ to that of N N ββ. We obtain where M 0ν is the NME of the 0νββ process, R is the nuclear radius and the phase factor G N N is given in Eq. (3) (replacing m N by m χ ). The translation of experimental measurements into constraints on g χ will hence need to rely on NME calculations.
If the Z 2 stabilizing symmetry is exact, χ could play a role in cosmology. For example, if it thermalizes with the primeval plasma, it may significantly alter the successful predictions of the standard Big Bang nucleosynthesis scenario [27]. These strong conclusions could be modified, for instance, if the Z 2 symmetry is only approximate, so that χ is stable within the detectors, but is not cosmologically long-lived. In this paper we will adopt a phenomenological approach, and we will focus on the search of exotic light fermions produced in pairs in double-β decays, without addressing their possible implications on cosmology that depend on additional parameters (such as the lifetime).

IV. DISCOVERY POTENTIAL OF FUTURE 0νββ ββ ββ DECAY EXPERIMENTS
A large experimental program has been mounted to search for the 0νββ decay of different isotopes and using different detection techniques [4]. The current-generation experiments set the ground for the identification of the  best techniques for further investment. The next-generation experiments are currently being proposed, designed, or in construction. Among the experiments proposed for the next decade, we will focus on LEGEND [28], nEXO [29] and CUPID [30]. We selected these experiments because their efficiencies and uncertainties in the search for massive fermions can be inferred from the 2νββ analysis published by their predecessors, i.e., GERDA [31], EXO-200 [32] and CUPID-Mo [33]. The parameters assumed for each experiment are listed in Tab. I. While the target isotope and the backgrounds vary across these experiments, the analysis to search for light exotic fermions is always conceptually the same. The energy window of interest goes from the detector threshold to the Q-value of the decay. In this window, the majority of the events is due to 2νββ decays (> 95% in the currentgeneration experiments). The other events can be due to a multitude of processes, for instance, natural radioactivity and cosmic rays. The most important parameters affecting the sensitivity of an experiment to light exotic fermions are the exposure E, the background rate R B and the systematic uncertainties due to the energy reconstruction σ sys . The exposure is given by the product of the number of observed nuclei and the observation time. The background rate is primarily given by the 2νββ decay rate -i.e. the inverse of the half-life -with a subdominant contribution due to the other sources, i.e. R B = R 2νββ + R other . The systematic uncertainties related to the event energy reconstruction can largely differ between experiments, but in general, their impact can be parameterized through an energy-dependent shape factor.
In the search for new physics, the upper limit on the rate of a new signal R U.L. scales approximately as the sum in quadrature of the statistical and systematic uncertainties. When the statistical uncertainty is dominant, the upper limit on the number of signal counts will be proportional to the fluctuations on the number of events in the analysis window, which is dominated by the background contributions. By comparing the upper limit on the number of signal counts (R U.L. · E) with the Gaussian fluctuations of the number of background counts ( √ R B · E), we obtain that: R U.L. ∝ R B /E. When the statistical uncertainty becomes comparable with the systematic one, the sensitivity starts to saturate and, eventually, further increasing the exposure does not improve the sensitivity.
To accurately quantify the sensitivity of the experiments, we have implemented a comprehensive frequentist analysis framework. Distortions of the double-β decay energy distribution due to light exotic fermions are searched through a binned maximum-likelihood fit based on a profile-likelihood test statistic [37]. Each process possibly contributing to the count rate in the energy window of interest should be added to the fit through a probability distribution function (PDF) that describes its expected event energy distribution. For this work, we use a PDF for the soughtafter signal and one for the dominant 2νββ decay. Both these PDFs are based on the calculations described in the previous sections. We use a third uniform probability distribution to account for other generic sub-dominant background sources. The actual shape of this third PDF affects only marginally our results. The parameters of the fit are typically the scaling factors of each PDF, i.e. the number of events attributed to each process 1 . The probability distributions of the test statistic are computed from large ensembles of pseudo-data generated under different hypotheses on the signal rate. This approach provides the right coverage by construction, including when the signal rate is close to the physical border.
Systematic uncertainties can bias the event energy reconstruction. Many detector-specific sources of bias should be considered. However, their overall impact can be parameterized through a shape factor with the form f (E) = 1 + a · E + b · E 2 + c/E where a, b and c are parameters that are considered to be known with limited accuracy (i.e. σ a , σ b and σ c ). Current-generation experiments are typically able to control the energy reconstruction bias at the per-cent level. To incorporate this systematic uncertainty in our analysis, we randomize a, b and c during the generation of the pseudo-data and run the fit using not-deformed PDFs [38]. The parameters are sampled from normal distributions centred at 0 and with a sigma of 10 −3 keV −1 , 10 −6 keV −2 and 10 −3 keV, respectively. We tested that this parameterization covers the maximal distortions estimated by the experiments. The result of this procedure is a broadening of the test statistic distribution and a reduction of the power of the statistical analysis. Figure 2 shows our projected sensitivity for a generic 76 Ge experiment and a sterile neutrino with a mass of 500 keV. The sensitivity for different background rates and exposure values is given in terms of the median 90% confidence level (C.L.) upper limit on sin 2 θ [37]. As mentioned above, the upper limit is proportional to the sum in quadrature of the statistical and systematic uncertainties. The computed sensitivities follow accurately the expected analytical expression: The parameters k 1 and k 2 cannot be computed without a full analysis. They depend on the isotope specific parameters (e.g. molar mass, phase space factor, Q-value and NME), the assumed mass of the sterile neutrino and the systematic uncertainties. They also capture the capability of the fit to differentiate between two spectral shapes. The vertical lines in the figure mark the exposure that has been collected by GERDA, as well as the exposure planned for the two phases of LEGEND. LEGEND has the potential to improve the sensitivity by an order of magnitude compared to GERDA, assuming that the systematic uncertainties can be kept below the statistical uncertainty. A significant improvement is expected even in the conservative scenario in which the systematic uncertainties are not reduced with respect to the current level.
After tuning the parameters k 1 and k 2 , Eq. (16) describes accurately the evolution of the sensitivity for all the considered experiments. It is also valid for the search of the Z 2 -odd fermions. In general, current-generation experiments are dominated by the statistical uncertainty while future experiments will need to fight against the systematic uncertainties. We will report projections based on the conservative scenario in which the systematic uncertainties will not be improved as well as the optimistic scenario in which they will be sub-dominant. These two projections define the ballpark for the sensitivity of future experiment and the results for intermediate scenarios can be interpolated from these two cases.
Our sensitivity projections for the search of sterile neutrinos and of the Z 2 -odd fermions are shown in Fig. 3. In both cases, the sensitivity evolution as a function of the fermion mass has a parabolic shape. Its minimum depends on the Q-value of the decay and corresponds to the fermion mass for which the experiment is most sensitive. The experiments quickly loose sensitivity towards vanishing masses. This is because the smaller is the fermion mass, the smaller is the spectral distortion. A similar loss of sensitivity occurs at larger masses, where the fermion mass approaches The experimental constraints are displayed through a band covering a range of sensitivities that goes from the most optimistic scenario (systematic uncertainty smaller than statistical uncertainty), till a conservative scenario (systematic uncertainty at the level of past analyses). Existing sterile neutrino constraints from single-β decay experiments [2,3,16,17] and solar neutrinos [39] are shown in the background. The right panel shows the median upper limit on the coupling constant between Z2-odd fermion and the neutrinos assuming the effective interaction in Eq. (14). The spread of the bands account for the systematic uncertainties as in the case of the sterile neutrinos, but additionally it covers also the full range of possible NME values found in literature [40][41][42][43][44][45][46][47][48][49][50]. the maximum energy available in the decay and the phase space shrinks 2 . The exposure, isotope properties, and experimental parameters define the offset of the curves. These projections are done assuming that each experiment can extend its analysis window down to arbitrarily small energies. If this is not the case, the upper edge of the curve will be lowered by the value of the energy threshold and the offset will slightly increase because of the reduction in detection efficiency. The sensitivity on the coupling constant g χ is computed from the decay rate of χχββ using Eq. (15). In addition to the systematic uncertainty considered for the sterile neutrino search, our constraints on g χ include also the uncertainties due to the NME calculations.

V. RESULTS AND OUTLOOK
The main results of our analysis are displayed in Fig. 3. The left panel shows the projected upper limits on the active-sterile mixing, parametrized as sin 2 θ, as a function of the sterile neutrino mass. We also show in the plot the current limits from single-β decay experiments [2,3,16,17] and from solar neutrinos [39]. As can be seen from the plot, the sensitivity of current double-β decay experiments is weaker than the existing limits, but only by a factor of a few. The larger exposure of future double-β decay experiments encourages a dedicated search for these exotic decay channels.
Indeed, our sensitivity study demonstrates the potential of future double-β decay experiments to improve the current limits on the active-sterile mixing angle in this mass range. Fig. 3 shows our projected sensitivities for LEGEND-1000, nEXO and CUPID, assuming no improvement in the systematic uncertainties with respect to the current-generation experiments (upper side of the bands), or assuming that the systematic uncertainties will be reduced below the statistical ones (lower side of the bands). In the most conservative scenario, the sensitivity of future searches will be comparable to the current limits. However, even with a modest reduction of the systematic uncertainties, next-generation experiments will explore uncharted regions of the sterile neutrino parameter space, even reaching sin 2 θ ∼ 10 −3 − 10 −4 for m N ∼ (100 − 1600) keV.
Double-β decay experiments have also the capability of probing models in which only the double production of light exotic fermions is allowed. This previously overlooked opportunity can lead to the first constraints on this kind of models. This is shown in the right panel of Fig. 3 for the Z 2 -odd singlet χ introduced in Sec. III. The sensitivities of current experiments for the effective coupling g χ lie between 10 −2 and 10 −3 MeV −2 , which could be improved up to 10 −4 MeV −2 with a favorable NME value and negligible systematic uncertainties.
Assuming that this effective interaction is originated at some scale Λ with O(1) dimensionless Wilson coefficient, this would imply that Λ 100 MeV. This is comparable to the typical momentum transfer in double-β decay experiments, and therefore the effective field theory approach employed in this work to recast the limits on the decay rate into limits on the coupling strength should be taken with a grain of salt. Yet, the search for distortions in the double-β decay spectrum due to the emission of two light exotic particles (fermions or scalars) is well motivated theoretically and deserves further investigation.
To conclude, in this letter we have proposed and explored the possibility of searching for the production of light exotic fermions using double-β decay experiments. We have derived the double-β decay spectrum when one or two light exotic fermions are emitted in the final state. We have also studied how these new channels can be searched for through a fit, and we have estimated the projected sensitivity for various future experiments considering the impact of both statistical and systematic uncertainties. We have found that double-β decay experiments constitute an promising avenue (and for some scenarios unique) to search for new physics. We therefore encourage the experimental collaborations to perform dedicated searches for this class of scenarios, properly including all the systematic and statistical uncertainties.
Note added. During the final stages of this project, we received a preprint by Bolton et al. [51], also exploring the possibility of searching for sterile neutrinos with double-β decay experiments. While they focus exclusively on sterile neutrinos (left and right-handed models) we considered generic single and double production of light exotic fermions. Their procedure to derive the sensitivity projections and their treatment of the systematic uncertainties differ from ours. Nevertheless, in the aspects where our analysis overlap, the results are qualitatively similar.