Probing Beyond the Standard Model Physics with Double-beta Decays

Nuclear double-beta decays are a unique probe to search for new physics beyond the Standard Model. Still-unknown particles, non-standard interactions, or the violation of fundamental symmetries would affect the decay kinematic, creating detectable and characteristic experimental signatures. In particular, the energy distribution of the electrons emitted in the decay gives an insight into the decay mechanism and has been studied in several isotopes and experiments. No deviations from the prediction of the Standard Model have been reported yet. However, several new experiments are underway or in preparation and will soon increase the sensitivity of these beyond the Standard Model physics searches, exploring uncharted parts of the parameter space. This review brings together phenomenological and experimental aspects related to new-physics searches in double-beta decay experiments, focusing on the testable models, the most-sensitive detection techniques, and the discovery opportunities of this field.

Double-β decays are nuclear transitions in which an isotope changes its atomic number Z by two units, while maintaining constant its mass number A: This process can occur only if the conversion of two protons into two neutrons leads to a more bounded nuclear configuration, with a positive Q-value defined as Q ββ ≈ M (A, Z)c 2 − M (A, Z + 2)c 2 [1], where M refers to the mass of the isotopes.Different theory models predict different decay final states.However, in general, they all envision two electrons to conserve the electric charge, along with neutral particles.These decays are a second-order weak process.Thus, they can be observed only in isotopes for which the otherwise dominant first-order transitions -and in particular, the single-β decay -are strongly suppressed.This is the case for a limited number of even-even nuclei for which the single-β decay is energetically forbidden, as the attractive nuclear pairing interaction makes them more bound than their odd-odd neighbors but less than their even-even second-neighbors (see Figure 1).Alternatively, candidate isotopes for the observation of double-β transitions are also those for which the single-β decay is suppressed by the large mismatch between the total angular momentum of the initial and final nuclei [2].Double-β decays were postulated in 1935 [3] by Maria Goeppert-Maier, who pointed out how these transitions could proceed through the "simultaneous emission of two electrons and two (anti)neutrinos": Due to the pairing term in the semi-empirical mass formula, β − transitions of even/even nuclei to their odd/odd isobaric neighbor can be energetically forbidden, whereas, in a secondorder process, β − β − decay is allowed.
In this case, the decay would "appear as the simultaneous occurrence of two transitions, each of which does not fulfill the law of conservation of energy separately".This final state is the only one allowed by the Standard Model of particle physics and is typically referred to as "two-neutrino double-β" (2νββ) decay.
The discovery of double-β transitions traces back to the observation of the decay daughter isotope in a geochemical experiment with 130 Te [4].Only 40 years later, double-β decays were observed in real-time through a calorimetric measurement of the energy released by the electron emitted in the decay of 82 Se [5].The electron summed energy distribution was found to have a continuous shape well compatible with the expectations for Goeppert-Maier's 2νββ decay (see Figure 2).To date, double-β transitions compatible with the 2νββ final state have been observed in 9 isotopes with half-life values in the range of 10 18 − 10 24 yr [6], making it one of the rarest processes ever measured by humankind.
Theories beyond the Standard Model (BSM) predict a variety of new, additional final states involving the emission of electrically neutral exotic particles, which would affect the decay kinematics and alter the energy and angular distributions of the electrons emitted in the decay.Similarly, also the violation of fundamental symmetries would affect the decay kinematic.This is why double-β decays have been a unique window to new physics for several years.However, despite many past and recent experiments, no deviations from the standard Model's 2νββ-decay expectations have yet been observed.Double-β decay searches have been historically carried out in experiments primarily built to study a hypothesized final state with only two electrons (i.e.without neutrinos), in which the summed electron energy is expected to be precisely equal to Q ββ rather than a continuous distribution (see Figure 2).Such "neutrinoless double-β" (0νββ) decay is indeed predicted by our leading theories explaining the matter-antimatter asymmetry and origin of neutrino masses.Its search has been an indisputable priority of the particle-physics community for the last two decades.Currently, ton-scale double-β-decay experiments are under preparation with the ultimate goal of increasing the sensitivity to 0νββ-decay half-life values by two orders of magnitudes compared to current constraints, up to 10 28 yr.We refer the reader to Ref. [7] for a recent review on 0νββ decay.
Thanks to their planned ultra-low background level and huge target mass, the next-generation double-β decay experiments will provide high-precision measurements of the 2νββ-decay, enabling new opportunities to discover new final states.This motivates the timing of our work.Indeed, although several review articles have discussed double-β decay searches in recent years, they have typically focused on the 2νββ-decay final state.For instance, we would like to highlight Ref. [8] that gives an excellent review of 2νββ-decay's history and [9] for a recent summary of the field.Differently from past works, our review is the first one focused on the BSM searches, covering their theoretical and experimental aspects, current constraints, and ongoing endeavor to improve the experimental sensitivities.Future 0νββ decay experiments with leading sensitivity will measure only the summed energy of the two final state electrons -aiming to distinguish the 0νββ decay peak at Q ββ from the continuous 2νββ decay distribution -without having access to the single electron energy or electron angular correlation.Therefore, this review focuses on how BSM physics would affect the summed electron energy distribution and how this can be exploited to search for BSM double-β decays.
The manuscript is organized as follows.We first briefly discuss the Standard Model 2νββ decay in Sec.II for completeness.We then review in Sec.III the phenomenology of double-β transitions mediated by Beyond-Standard-Model processes.Section IV is dedicated to their experimental signatures exploited by experiments.It also contains an original contribution to the field, which is the first derivation of analytical formulas describing the experimental sensitivity as a function of key experimental parameters and systematic uncertainties.Finally, Sec.V summarises the latest experimental results and the prospects for next-generation experiments.

II. STANDARD-MODEL ALLOWED DOUBLE-β DECAY
Double-β transitions can be uniquely identified by the production of the daughter nucleus.With this principle, the double-β decay of 130 Te was first detected in 1950 with the first geochemical experiment.The detection of an excess of 130 Xe was proof of the double-β decay of the initial nucleus and allowed a first determination of its half-life [4].Even if this result was initially not considered seriously, it represents the first observation of 2νββ decay, as it became clear after 15-20 yr.
In fact, after both the 2νββ and 0νββ decays were proposed (1935)(1936)(1937)(1938)(1939), the first estimate of the half-life of the two processes strongly favored the 0νββ decay.This was predicted with a half-life of the order of 10 15 yr, whereas a half-life of about 10 21 yr was predicted for the Standard-Model allowed 2νββ decay.Only the discovery of parity violation in 1957 and the determination of the V-A nature of weak interaction made clear that the probability of 0νββ decay had to be much smaller than that of 2νββ decay.
In the following years, new geochemical experiments were performed, confirming the observation of the 2νββ decay of 130 Te [10,11] and observing for the first time the 2νββ decay of 82 Se [11] and 128 Te [12].
A key milestone in the history of double-β decays is the theoretical work of Doi, Kotani, and Takasugi in 1985 [13], who first computed the energy and angular distributions of the electrons emitted in double-β decays, providing a clear signature to distinguish the Standard-Model allowed 2νββ decay from BSM double-β decays.In the same work, Figure 2 appeared for the first time.
Following these theory developments, the first direct observation of a 2νββ decay was made in 1987 using a 82 Se time projection chamber, which could measure the summed energy distribution of the two electrons and extract the 82 Se 2νββ-decay half-life based on 36 observed events.This was a turning point in the history of 2νββ decays, and in the following ten years 2νββ decay was directly observed in seven isotopes:  96 Zr [26].In addition to the summed electrons' energy distribution, which was measured for all the above-mentioned isotopes, the NEMO-2 experiment also measured the single electron energy and the angular distributions of the electrons for 100 Mo, 116 Cd, 82 Se, and 96 Zr [18,20,22,24,26].In the same years, the 2νββ decay of 238 U was first observed in a radiochemical experiment [27].
By the end of the 20th century, the energy and halflife of 2νββ decay had been measured for 10 isotopes, considering direct, geochemical, and radiochemical experiments.In addition, in the last two decades, direct observation of the 2νββ decay of 130 Te and 136 Xe were reported in 2010 [28] and 2011 [29], respectively.Table I summarises the double-β decaying isotopes for which the 2νββ decay was directly observed.
The precision with which experiments were able to determine the half-life of 2νββ decays has increased over the years, with modern experiments reaching a percentlevel precision for most of the isotopes.The most precise measurement of the 2νββ decay half-life are summarised in Table I.
In addition to the Goeppert-Maier's 2νββ decay (Equation 2) or β − β − transition, depending on the relative numbers of protons and neutrons in the nucleus, three additional second-order transitions are allowed in the Standard Model: The energy released in the three processes listed above is smaller compared to the β − β − decay in Equation 2. Consequently, these processes have lower probabilities compared to the β − β − decay due to the smaller phase space, and experimentally they are much more challenging to observe.In the following, we will always refer to the β − β − process as double-β decay.
In nature, there are 35 isotopes that can undergo 2νββ decay, and 34 more that can undergo β + β + , ECEC, and ECβ + [39].In fact, not all of them fulfill the experimental requirements, e.g., high isotopic abundance, large Q ββ value, and compatibility with experimental technologies.To date, only the nine isotopes listed in Table I were used in direct search experiments.
The rate of 2νββ decay can be calculated following Fermi's golden rule for β decay.To a good approximation, the kinematic part (phase space of the leptons emitted in the decay) and the nuclear part (matrix element responsible for the transition probability between the two nuclear states) can be factorized as: where G 2ν is the phase-space factor and is obtained by integrating over the phase space of the four leptons, and M 2ν is the nuclear matrix element (NME) and deals with the nuclear structure of the transition.While the phasespace factor can be calculated exactly, the NME is much more difficult to evaluate and relies on nuclear structure models.
Experiments measure the distribution of the summed kinetic energy of the two electrons (K): To a first approximation, the shape of this distribution is determined only by the phase space.The contribution of the NME to it is small and primarily affects the absolute value of the transition probability.NEMO-3 [38] Table I.Direct observations of 2νββ decays.The year of the first observation is indicated for each isotope, together with the most precise half-life determination.The shown uncertainty is the sum in quadrature of the statistical and systematic uncertainties, when available.
The summed energy of the two electrons emitted in the 2νββ decay is continuously distributed between 0 and the endpoint at the Q ββ value due to the neutrinos escaping the detector and carrying away part of the energy 1 .This is shown in Figure 2, compared to the 0νββ decay, for which a δ function at Q ββ is expected because all the transition energy goes into the kinetic energy of the two electrons.
In conclusion, we have seen that even if the first observation of 2νββ decay is commonly traced back to 1950 and the first geochemical experiment with 130 Te, it took many years for the community to convince themselves that the production of the daughter nucleus observed in this first experiment was exactly the result of Goeppert-Maier's 2νββ decay.To date, precision measurements of the 2νββ decay of several isotopes and the agreement between the distribution of multiple observables (electrons summed energy, single electron energy, and angular distributions) with their theoretical prediction are striking evidence for the 2νββ decay and exclude a large part of the parameter space for many BSM theories.

III. BEYOND-STANDARD-MODEL PHYSICS IN DOUBLE-β DECAY
This section reviews double-β decays in the presence of BSM physics.
Suppose new particles were involved in a double-β decay, or any new physics affected the phase space of the two electrons emitted in the 2νββ decay.In that case, the expected distribution of the summed electron energy would differ from that predicted for the SM 2νββ decay (equation 5, figure 2).This is the main feature used to 1 To be precise, the maximum energy is Q ββ minus the mass of the two emitted neutrinos.This is typically neglected as Q-values are at the MeV-energy scale, while neutrino masses are smaller than the eV-energy scale [40,41].
search for these BSM decays in the experimental data.Single electron energy distributions and electron angular distributions are also primarily determined by the phase space, and they are characteristic of the physics model and decay final state.Experiments that can also measure these distributions can strongly enhance their sensitivity in distinguishing the SM allowed 2νββ decay from BSM double-β decays.This review mainly discusses the first aspect, given that the calorimetric measurement approach pursued in the next generation of 0νββ decay experiments with leading sensitivities gives access only to the summed electron energy distribution.We classify the models into three groups.The first one contains those predicting the existence of new particles -either bosons or fermions -which are emitted in the decay, replacing one or both of the 2νββ's neutrinos.The second one includes theories in which fundamental symmetries such as Lorentz covariance or Pauli's exclusion principle are violated.The last group covers nonstandard interactions, like right-handed leptonic currents and strong neutrino self-interactions.
A. New particles

Bosons
In the early 80s, an attractive approach to the neutrino mass problem was considered in which the neutrinos are Majorana particles with small masses arising from the spontaneous breakdown of the global B-L symmetry.In these models, a massless Goldstone boson should exist, which was called the "Majoron".Several realizations of this idea were proposed, mainly differing by the weak isospin (I) properties of the Majoron and leading to different phenomenology.
The first model was proposed by Chikashige, Mohapatra, and Peccei [42,43].In this model, the Majoron arises from a Higgs singlet (I = 0) and gives rise to small neutrino masses via the "see-saw" mechanism.On the other hand, the Majoron coupling to neutrino is so small that it would be very hard to test it through laboratory experiments.Shortly after, Gelmini and Roncadelli proposed a model in which the Majoron arises from a Higgs triplet (I = 1) [44].In this second case, a stronger coupling to neutrinos would be possible.A third case was considered later in 1987, in relation to solar neutrino oscillations, where the Majoron arises from a Higgs doublet (I = 1/2), possibly leading to strong coupling to neutrinos [45][46][47].
If the Majoron Yukawa coupling to neutrinos were sufficiently strong (∼ 10 −5 −10 −3 ), this would have interesting consequences for particle physics, astrophysics, and cosmology.Already in 1981, starting from the Gelmini-Roncadelli Majoron model, a rich phenomenology was derived [48].One of the most interesting consequences considered at that time was the possibility of emitting such a Majoron in double-β decays, giving rise to a new final state in which two electrons and a Majoron (and no neutrinos) are present.The decay rate of this process was calculated, and constraints on the neutrino-Majoron coupling were set using existing limits on the decay rate of 0νββ decay [48].In fact, at that time, the possibility of distinguishing 2νββ, 0νββ, or double-β decay with the emission of a Majoron was not conceived, and the latter was regarded only as "an interesting possibility, which may confuse the analysis of double-β decay experiments" [49].Only the impressive theoretical work of Doi, Kotani, and Takasugi in 1985, in which the energy and angular distribution of the electrons emitted in the double-β decay with the emission of a Majoron were calculated for the first time [13], provided a clear experimental signature to distinguishing different doubleβ decay channels.The same authors revised and updated these calculations in a successive work [50].In the same years, a supersymmetric model was developed, which would lead to the emission of two Majorons in double-β decays [51].
It was clear since early works that the Majoron with nontrivial weak isospin properties (such as the triplet and doublet Majoron), which would have appreciable coupling to neutrinos, would also have a strong coupling to the other leptons and would necessarily be discovered in the upcoming lepton-lepton collisions machines [48].As the first LEP data on the invisible width of the Z boson became available in the early 90s, the number of active light neutrino species was limited to three, ruling out both triplet and doublet Majoron models [52,53].
On the other hand, in the same years, the first direct double-β decay experiments started to search for double-β decay with Majoron emission, and an excess of events below Q ββ was observed, which could be compatible with the energy distribution predicted by Doi Et al. for such decay [54][55][56].The only available Majoron model not contradicting LEP data was the singlet Majoron model.Nevertheless, in its original Chikashige-Mohapatra-Peccei formulation, the Majoron is so weakly coupled to neutrinos that it cannot produce detectable effects in double-β decay and, therefore, not explain the experimental data.
The sum of these events motivated in the following years the flourishing of a number of new models able to reconcile the results on the Z decay width with a neutrino-Majoron coupling strong enough to explain the event excess: new models with a singlet Majoron [57], models in which the Majoron carries a non-zero lepton number [58,59], models predicting the emission of two Majorons [60], and models in which, departing from the original conception of the Majoron as Goldstone boson, the Majoron arises as the component of a massive gauge boson [61] or a bulk field [62].
This ensemble of models can lead to two different final states, corresponding to the emission of one or two Majorons, which we will indicate with J: The rate of the double-β decay with the emission of one or two Majorons can be expressed as: where g Jα is the neutrino-Majoron coupling, M J(JJ) α the NME, and G J(JJ) α the phase-space factor.All three terms depend on the particular model, which we indicated with the subscript α.Systematic calculations of phase-space factors and NMEs for a number of Majoron models were performed in [63] for many isotopes.More recently, improved calculations of the summed electron energy distributions have been performed, leading to improved calculations of the phase-space factors [64].Improved calculations of the NMEs have also been performed lately [65].
If one or two Majorons are emitted in the double-β decay, they would escape any detector and carry away part of the decay energy.In analogy with the 2νββ decay, the summed electron energy is continuously distributed between 0 and Q ββ , and its exact shape is primarily determined by the phase space.In turn, the phase space depends on the Majoron model, particularly on the effective neutrino-Majoron interaction Lagrangian leading to the Majoron-emitting double-β decay.This can be parameterized to-a-first-approximation with a spectral index n: Figure 3 shows the summed electron energy distribution for different Majoron models compared to the SM 2νββ decay distribution.
The spectral index is commonly used to group models predicting the same experimental signature, i.e. the same summed electron energy distribution.These models are not distinguishable by the experiments.Table II shows a summary of all the Majoron models grouped by Table II.Different Majoron models which predict double-β decays with the emission of one or two Majorons.The third column indicates the model's spectral index (n), the fourth column indicates whether the Majoron is a Goldstone boson or not, and the last column indicates the leptonic charge (L) of the Majoron.Models with a leptonic charge different from zero preserve the lepton number symmetry.
the number of emitted Majorons in the second column, the spectral index in the third column, and the Majoron's properties in the last two columns.The fourth column indicates whether the Majoron is a Goldstone boson or not, whereas the last column shows the Majoron's leptonic charge.Models in which the Majoron carries a lepton number different from 0 preserve the lepton number symmetry.Without an independent test of Lepton number violation, these models are experimentally indistinguishable from the corresponding lepton number nonconserving processes.
All the models discussed so far focus on the emission of one or two Majorons originating from the intermediate neutrino exchanged in the process.Despite the differences among the models, i.e. the different effective neutrino-Majoron interaction Lagrangians, all of them assume the SM V-A structure of the charged currents involving leptons and quarks.We will refer to them as "classical" models.Recently a new scenario has been considered, in which a Majoron-like particle (φ) is emitted in the double-β decay (φββ decay).In this case, the interaction is described by an effective dimensionseven operator, with right-handed lepton current and right/left-handed quark current [66].The Feynman diagram of this process is shown in figure 4 in comparison to the Majoron emission in classical models.The coupling strength between the neutrino and the Majoron-like φ is RL if the effective operator contains left-handed quark current, and RR when the effective operator contains right-handed quark current.The two cases have been considered separately, with only one of the two operators being present at a time.
The energy distribution predicted for the φββ decay is also shown in figure 3. The distribution associated with RL is very similar to the classical Majoron emission models leading to n = 1.On the other hand, introducing a hadronic right-handed current in the RR term changes the shape of the distribution considerably.
In all the previous discussions, we always assumed the Majoron to be massless.However, many of the models already presented do not prevent the Majoron from being a light particle [60,66,67].This possibility became extremely popular because light Majorons could be a dark matter candidate [68,69].If the Majoron mass is below the Q ββ , double-β decay with the emission of a Majoron can still happen.In this case, the end-point of the energy distribution is shifted to Q ββ − m J , where m J is the Majoron mass.

Fermions
In many extensions of the Standard Model, new spin 1/2 particles, singlet under the Standard Model gauge group, are introduced in relation to the question of neutrino mass generation or dark matter.Currently, the most popular exotic fermion is the sterile neutrino N .Sterile neutrinos are neutral and right-handed Standard Model singlet fermions that interact with ordinary matter only through mixing with the active neutrinos.We refer the reader to Ref. [70] for a recent review of the theoretical and experimental motivation for sterile neutrinos, as well as their phenomenological consequence.In a variant of this scenario, the singlet fermion could be furnished with a Z 2 symmetry, so it can only be produced in pairs.In general, when a new exotic fermion is introduced in the theory, its mass and coupling to the Standard Model particles are free parameters of the model.It is left to laboratory experiments and astrophysical and cosmological observations to probe the vast allowed parameter space.
In 1980, Shrock examined the possibility of searching for sterile neutrinos in β decays [71].The admixture of one or more sterile neutrino states would create a discontinuity in the β decay spectrum similar to the discontinuity that a non-zero neutrino mass is expected to create at the endpoint.The position and amplitude of this kink would give information on the mass of the sterile neutrino and its mixing with the active neutrinos.Since then, several β decay experiments have searched for sterile neutrinos and, still today, they set the most stringent constraints in the mass range between ∼ 10 eV and ∼ 1 MeV [72][73][74][75][76].
In analogy with the case of single-β decays, sterile neutrinos with masses below few MeV could be produced in double-β decays.This possibility was recently discussed in [77,78].In [78], the production of exotic fermions in double-β decays was also extended to models in which the single production is forbidden by additional symmetries while the pair production is allowed.Such a model could be realized with the neutral singlet fermion χa potential dark matter candidate -being charged under a discrete Z 2 symmetry to make it stable.This new fermion could interact with neutrinos through an effective four-fermion scalar interaction of the form g χ ννχχ.This case is particularly interesting because such a particle cannot be produced in single-β decays, and double-β decays represent a unique discovery opportunity for laboratory experiments.
In general, models predicting the existence of light exotic fermions coupling with the Standard Model neutrinos can lead to two additional double-β decay final states, corresponding to the emission of one or two exotic fermions, which we will indicate with f : Sterile neutrinos can be produced via both 9a and 9b .Summed electron energy distributions of the doubleβ decay into one sterile neutrino with a mass of 600 keV and the double-β decay into two Z2-odd fermions with a mass of 300 keV in comparison to the SM 2νββ decay distribution for the case of 76 Ge isotope.An arbitrary normalization is used for illustrative purposes.Adapted from [78] decay channels.Provided that both decay channels are kinematically allowedi.e. the Q ββ value must be larger than the sterile neutrino mass for the 9a decay channel, and twice the sterile neutrino mass for the 9b decay channel -the total double-β decay rate would become an incoherent sum of three channels: where sin 2 θ represents the mixing angle between active and sterile neutrinos.The first term accounts for the SM 2νββ decay, the second one for the decay in which one of the two neutrinos is replaced by a sterile neutrino (equation 9a with f = N ), and the last one for the decay into two sterile neutrinos (equation 9b with f = N ).We shall notice here that this last term is strongly suppressed by a factor of sin 4 θ, making it negligible for experimental searches.
Each term of the sum can be expressed factorized as the product of the NME and the phase space factor: where a and b indicate the emitted particles, either the active or the sterile neutrinos.While the NME is the same for the three terms, the phase space factor is affected by the presence of one or two sterile neutrinos in the final state, and so is the summed electron energy distribution: The presence of a massive sterile neutrino in the final state affects the kinematics of the decay such that the endpoint of the summed electron energy distribution is shifted at Q ββ − m N , where m N is the sterile neutrino mass.The energy distribution for the double-β decay into one sterile neutrino with a mass of 600 keV is shown in figure 5 for the case of 76 Ge isotope.
Z 2 -odd fermions, which we refer to as χ fermions, can be produced in double-β decay via the 9b channel.To be kinematically allowed, the Q ββ value of the decay must be larger than twice the mass of the χ fermions.The rate of the double-β decay with the emission of two fermions χ can be expressed as: where g χ is the coupling between neutrinos and the χ fermions, m e the electron mass, and R the nuclear radius.The NME for this decay can be take to good approximation as the NME for the 0νββ decay, M 0ν .The presence of two massive χ fermions in the final state affects the kinematics of the decay as for the case of sterile neutrinos: the endpoint of the summed electron energy distribution is shifted by Q ββ − 2m χ , where m χ is the χfermion mass.This is shown for m χ = 300 keV in figure 5 for the case of 76 Ge isotope.

Lorentz violation
Lorentz invariance is one of the fundamental symmetries of the SM of particle physics.The breakdown of Lorentz and CPT symmetries at the Plank scale is an interesting feature of many theories of quantum gravity, such as string theory [79].Despite direct studies of physics at this ultrahigh energy scale are far from the reach of current accelerator-based experiments, some suppressed effects could arise at lower energies and be observable with the actual experimental technologies.
The general framework that characterizes Lorentz violation in the SM is the Standard Model Extension (SME) [80,81].This is an effective quantum field theory that includes all possible operators that can be constructed with the SM fields and that introduce Lorentz violation but preserve the SM gauge invariance.The development of the SME has led to experimental searches for Lorentz violation in all different sectors of physics, including matter, photon, neutrino, and gravity [82,83].A data table of the current constraints is compiled in [84] annually.
The behavior of neutrinos in the presence of Lorentz and CPT violation has been extensively studied using the SME framework, and its related coefficients classified [85][86][87].Most of these coefficients are currently constrained by neutrino oscillations experiments.However, four coefficients only affect the neutrino phase space and escape detection through the measurement of neutrino oscillations.The corresponding operators are an example of counter-shaded Lorentz-violating operators and they might have escaped detection to date even though their effect is large compared to the suppression given by the Plank scale [82].These coefficients, commonly referred to as oscillation-free (of) coefficients, can be studied in weak decays, such as single-β decay or double-β decay [88,89].
The interaction of neutrinos with the counter-shaded operator modifies their four-momentum: where a (3) of encodes the oscillation-free coefficients and å (3) of is the isotropic component å of ) 00 / √ 4π.Considering this modification, the decay rate of the 2νββ decay in the SME framework can be written as the sum of two terms: where Γ SM is the SM decay rate and δΓ LV is the perturbation term due to the introduction of Lorentz violation.Similarly to the previous cases, we can factorize each term as the product of the NME and the phase space factor (see equation 11).Lorentz violation does not affect the NME but appears as a kinematic effect modifying the phase space factor, thus the summed electron energy distribution: The modification of the phase space d(δG LV )/dK comes from the change of the differential element of the antineutrino momentum: The integration over all anti-neutrino orientation performed to obtain the summed electron energy distribution in the case of double-beta decays implies that only isotropic effects are observable.Hence, the spectrum only depends on å(3) of .The energy dependency in the phase space of the perturbation term can be approximated as δG LV ∼ (Q ββ − K) 4 .Using the same terminology introduced for the Majoron, the spectral index of this perturbation is n = 4.
On the other hand, the spectral index of the SM term is n = 5.Therefore, a non-zero value of the coefficient å (3) of , which implies a non-zero contribution of the perturbation term, produces a distortion of the spectrum of double-β decays compared to the SM expectation.The energy distribution of the Lorentz violating perturbation term is shown in figure 6 compared to the SM 2νββ decay distribution.

Violation of Pauli exclusion principle
Neutrinos have many peculiarities among all the known particles.They are the only neutral leptons, which leaves the possibility for neutrinos to be Majorana particles.In addition, the smallness of their mass points to a different mass-generating mechanism compared to the standard coupling with the Higgs boson exhibited by all other particles.Therefore, neutrinos might have substantially different properties compared to the charged leptons.
Pauli's original formulation of the exclusion principle in 1925 postulated that two or more identical electrons cannot occupy the same quantum state within a quantum system simultaneously.This was successively extended to all fermions with half-integer spin and experimentally confirmed with extremely high precision.From the theoretical point of view, a local quantum field theory with violation of the Pauli principle has been discussed [90][91][92], and some difficulties highlighted [93,94].
The consequences of a change of neutrino statistics from fermionic to bosonic would have substantial cosmological and astrophysical consequences [95][96][97][98].However, even recent analyses of available cosmological data can set only weak bounds on neutrino statistics [99].With two identical anti-neutrinos in the final state, doubleβ decay is a unique process to test the violation of the Pauli's principle [98,100].Qualitative conclusions in [98] on 2νββ decay ruled out a pure bosonic neutrino, but not the possibility that neutrinos obey non-standard statistics, more general than Bose or Fermi ones [101].
If neutrinos obey a mixed statistic, the neutrino's state would be the combination of fermionic and bosonic states.As a consequence, the double-β decay amplitude would be given by the sum of two terms, corresponding respectively to the fermionic (anti-symmetric) and bosonic (symmetric) parts of the two anti-neutrino emissions: In the phase-space integration, the interference between the anti-symmetric and symmetric parts of the amplitude vanishes and the 2νββ decay rate becomes fermionic and pure bosonic neutrinos.In the decay rates for pure fermionic and pure bosonic neutrinos, both the kinematics terms and the NMEs are different.Defining the ratio the normalized differential decay rate can be written as The ratio r 0 determines the weight with which the bosonic component enters the total rate and the differential decay distribution.If r 0 is very small, a substantial modification of the energy distribution is expected only for sin 2 χ being very close to 1.In addition, the ratio r 0 needs to be calculated and depends on the values of the NMEs.Thus, it introduces an uncertainty due to the nuclear-structure calculations.
On the other hand, the normalized differential decay rate for pure fermionic dΓ f /Γ f and pure bosonic dΓ b /Γ b neutrinos does not depend on any nuclear model assumption and are shown in figure 7. The spectrum for bosonic neutrinos is softer, with the maximum shifted at lower energy by a factor of about 15%, compared to the pure fermionic spectrum.
The calculations performed in [100] predicts the ratio r 0 for 100 Mo and 76 Ge to be 0.076 and 0.0014, respectively.The small ratio predicted for 76 Ge limits the sensitivity of double-β decay experiments with 76 Ge to spectral distortions due to a partly bosonic neutrino.
Experimental searches for bosonic or partly bosonic neutrinos with double-β decay experiments could use not only the shape of the distributions but also the ratios between the rates of the transitions to the excited states and the ground states if the first were observed [100].C. Non-standard interaction

Right-handed leptonic currrents
The SM 2νββ decay is a second-order transition involving weak left-handed V-A currents with the strength given by the Fermi constant G F .Some BSM theories, such as Left-Right symmetric models with unbroken lepton number [103,104], predict the existence of V+A lepton currents, which can mediate double-β decays [102].The new physics effects can be modeled through effective charged current operators containing V+A lepton currents.The strength of these non-standard interactions is given by XR G F , where the small dimensionless coupling XR encapsulates the new physics effects.
Right-handed current interactions are independent of the Majorana or Dirac nature of neutrinos and do not necessarily violate the lepton number.If the neutrino is a Majorana particle, the operators associated with LR and RR violate the total lepton number by two units and give rise to extra contributions to the 0νββ decay [105].In this case, 0νββ decay searches set stringent limits of the order LR 3 × 10 −9 , RR 6 × 10 −7 [106].On the other hand, if there exists a sterile neutrino state ν R that combines with ν L to form a Dirac neutrino, the right-handed current interactions do not necessarily violate lepton number [104].The strong theoretical interest is therefore supported by the fact that their observation, along with the non-observation of lepton number violation, would indicate that neutrinos are Dirac fermions [104].
Direct experimental constraints on these operators are set by neutrons and different single-β decays but are rather feeble ( LR , RR 6 × 10 −2 ) [107].Searches at the LHC are also possible [108][109][110] but are generally model dependent and require some caveat on the use of the effective operator analysis at high energies.
In the presence of right-handed leptonic currents, the amplitude of the 2νββ decay would be calculated as a coherent sum of the three Feynman diagrams shown in figure 8: the SM second-order transition with two lefthanded interactions with the strength given by G 2 F (figure 8a), a transition involving one right-handed interaction with strength XR G 2 F (figure 8b), and a secondorder transition with two right-handed interactions with An arbitrary normalization is used for illustrative purposes.
Adapted from [102].8c).Nevertheless, to the lowest order in the exotic coupling XR , the decay rate can be expressed as an incoherent sum of only two terms: where the first term is the SM decay rate and the second term is the contribution of right-handed current to the decay rate, suppressed by the coupling XR .In fact, the interference of the SM term (diagram 8a) with the diagram 8b is helicity suppressed by the masses of the emitted electron and neutrino as m e m ν /Q 2 ββ because of the right-handed nature of the exotic current.Higher orders in the exotic coupling XR , coming from the last diagram 8c and its interference with the SM term, are also negligible.
The phase-space factor and the NME differ in the SM decay rate and the BSM contribution.Thus, the presence of right-handed currents in double-β decay changes the total decay rate and the shape of the energy spectrum.Nevertheless, given the uncertainties in the NME calculations, the change in the total decay rate is not expected to be measurable.Instead, experiments may be sensitive to the change in the spectral shape.Figure 9 shows the 2νββ decay distribution in the SM compared to the distribution arising from the presence of right-handed currents.The deviation includes a spectrum shift to smaller energy and a flatter profile near Q ββ .

Neutrino self-interaction
The Hubble tension indicates the discrepancy between CMB and local measurement of the Hubble constant.This tension has grown to about 4σ, and if confirmed, it would require new physics BSM or a new cosmological model [111,112].
Introducing a νSI, i.e. a four-neutrino contact interaction, could resolve the Hubble tension.Such a νSI interaction can be written as G S (νν)(νν), and it would inhibit neutrino free-streaming in the early Universe if its strength is much larger than the Fermi effective interaction predicted by the SM, G S ∼ 10 9 G F [113,114].This new strong interaction would indicate the presence of new physics at a scale 1/ √ G S ∼ 10 MeV -1 GeV.In general, these strong νSI interactions are difficult to probe in laboratory experiments due to the absence of electrons or quarks involved.With some assumption on the origin of the νSI operator, constraints can be obtained from different physics observations [115,116], while no modelindependent constraint is currently available.The study of νSI in single-β decays has been considered [117].More recently, the search for νSI in double-β decays has also been proposed [118].
In the presence of νSI, independently of the Dirac/Majorana nature of neutrinos, the two neutrinos in double-β decay can be emitted via the corresponding effective operator, resulting in a νSI-induced 2νββ (2νSIββ) decay.The Feynman diagram of this process is shown in figure 10.The final state of the 2νSIββ decay is identical to that of the SM 2νββ decay.The contribution from νSI to the decay rate can be written as: where m e denotes the electron mass and R the radius of the nucleus.For an exact contact interaction of four neutrinos and neglecting the final state lepton momenta, the phase-space factor for the 2νSIββ decay is related to the phase-space factor of the 2νββ decay as The NME of 2νSIββ is the same as of 0νββ.In this scenario, no difference is expected in the summed electron energy distribution of the 2νSIββ decay compared to the SM 2νββ decay.Therefore, only the experimental measurements of the 2νββ decay rate can be used to constrain the contribution of νSI.This approach was used in [118] to determine upper limits on the coupling G S from the measured 2νββ decay rates of several double-β decay isotopes.Limits in the range G S /G F (0.32 − 2.50) × 10 9 were obtained.The sensitivity on G S is limited by the uncertainty of the NME ratio |M 0ν |/|M 2ν |.Cosmological data favoured a strong interactive regime with G S = 3.83 × 10 9 G F .Even including the theoretical NME uncertainties, all the considered isotopes can fully exclude the strongly interacting cosmologically favored regime [118].However, one should note that this bound applies only under the assumption that two electron neutrinos are involved in the νSI.This might not be the case if only muon neutrinos and tau neutrinos participate in νSI.
Possible distortions of the electron energy distribution could arise from the νSI contribution if the νSI operator were generated by light mediators.In this scenario, the energy dependence of the coupling G S could cause observable spectral distortions.In [118], the simplest case of an s-channel scalar mediator with a mass just above the kinematic threshold (M = Q ββ + 0.1m e ) was discussed.The coupling G S acquires the following energy dependence where M is the mediator mass and s ≡ p 2 , with p being the momentum of the mediator (in the context of the 2νSIββ, this is of the order s Q 2 ββ ).The value of G S at zero momentum transferred (G 0 S ) is denoted as G 0 S = g 2 /M 2 , with g the coupling between the mediator and the neutrino.Using G S in equation 24, the differential decay rate of the 2νSIββ decay can be calculated.The corresponding summed electron energy distribution is shown in figure 11, for a mass of the mediator M = Q ββ + 0.1m e .The energy spectrum of the 2νSIββ decay is shifted at lower energy compared to the 2νββ decay spectrum.This shift can be understood qualitatively: with the summed energy of the two electrons increasing, the energy available for the neutrinos is smaller, leading to a smaller value of s and hence a smaller value of G S .

IV. STATISTICAL ANALYSIS AND SENSITIVITY
All the double-β decays introduced in section III result in the same event topology.The two electrons emitted in the decay are detected within a small detector region.The additional particles produced in the double-β decay process, either the two anti-neutrinos or one or more exotic particles, escape the detector carrying away part of the decay energy.Double-β decay experiments typically measure multiple observables for each event.These include for instance, the energy deposited within the detector, spatial and timing information, and variables related to the type of particles involved in the event.Signal and background events feature specific values of these observables that can be used to separate them.As previously mentioned, this review focuses on the energy deposited within the detector -the summed energy of the two electronswhich will be measured by all future 0νββ decay experiments.However, we must recall that any additional information, like the energy distribution of the single electrons or their angular correlation, can strongly enhance the sensitivity to BSM physics as it was shown in past searches [30,119,120].
Before discussing the statistical analysis to search for BSM decays, we shall divide them into three classes of models, which require slightly different statistical treatments.
a. Class I models.The first class contains all those BSM decays that are independent of the SM 2νββ decay.This condition is satisfied if the predicted phase space and NME for the BSM decay are different from the SM ones and if no correlation is expected between the BSM decay and the SM 2νββ decay.This is the case of the decays involving Majorons (III A 1) or Z 2 -odd massive fermions (second model in III A 2), but also the case of 2νββ decay realized via non-standard interaction (III C), in the assumption that any interference term between the SM 2νββ decay and the non-standard decay can be neglected.Despite the differences among the models, the resulting BSM decay would manifest as an additional contribution to the total energy spectrum, given by a continuous distribution in the same energy range of the SM 2νββ decay.Being the BSM decay subdominant, a small distortion of the observed energy spectrum compared to the SM expectation would be observed.This is shown in the first left panel of figure 12, for the Jββ decay (n=1) and a non-realistically large coupling g J , as an example.In the search for such a BSM decay, the parameter of interest is the number of BSM decay events, underlying the respective distribution.On the other hand, the dominant SM 2νββ decay acts as a background.
b. Class II models.The second class contains those BSM decays which are not independent of the SM 2νββ decay.This is the situation in which the BSM model not only leads to an additional BSM decay but also affects the total decay rate of the SM 2νββ decay.This is the case of the decays involving sterile neutrinos or bosonic neutrinos.The correlation between the SM 2νββ decay and the BSM one is given by the respective mixing, and the existence of the BSM decay also implies a decrease of the SM 2νββ decay by a factor of cos 4 θ and cos 4 χ, respectively.This is shown in the central panel of figure 12, for a sterile neutrino with mass m N = 500 keV and a large mixing of sin 2 θ = 0.15, as an example.Again, a small distortion of the observed energy spectrum compared to the SM-only expectation would be observed.In the search for such an BSM decay, the parameter of interest is proportional to the ratio between the number of events underlying the BSM decay distribution and the number of events underlying the 2νββ decay distribution.Therefore, in this case, the dominant SM 2νββ decay signal helps to constrain the BSM physics.c.Class III models.Some BSM theories do not predict the existence of additional decay, as in class I and class II models, but alter the prediction for the 2νββ decay.We classify these decays as class III.This is the case, for example, of the violation of Lorentz symmetry.The Lorentz violating term in the neutrino momentum affects the phase space of the two electrons emitted in the 2νββ decay and, therefore, the predicted energy distribution.This is shown in the right panel of figure 12, for non-realistically large values of the Lorentz violating coefficient å(3) of .The search for Lorentz violation in 2νββ decays corresponds to a search for deviations of the ob-served two-electron energy distribution compared to the SM expectation.In practice, the total decay rate for the Lorentz violating 2νββ decay can be approximated, to the first order, as a sum of the SM 2νββ decay rate and a perturbation term, whose size is regulated by the coefficient å(3) of , which is very small.While the phase space differs for the two terms of the sum, the nuclear part, the NME, is the same, introducing a correlation between the two terms.From the statistical point of view, class III models can be treated as class II models.The parameter of interest is proportional to the ratio between the number of events underlying the distribution of the LV perturbation and the number of events underlying the SM 2νββ decay distribution.Again, the dominant SM 2νββ decay signal helps to constrain the BSM physics.

A. Statistical analysis
In the search for spectral distortions due to the contribution of an BSM decay in the energy spectrum, the energy region of interest extends from the detector threshold to Q ββ .In this window, most of the observed events are attributed to the 2νββ decay N 2ν with additional contributions due to different background processes N others .The sum of the 2νββ decay and other background constitutes the so-called background model, which is known with a certain accuracy by the experiments.
To search for an BSM decay, a spectral fit of the energy spectrum is performed, 2 adding the BSM decay to the background model.Typically, the information from the background model is used to construct a likelihood function, which is then used in a frequentist or Bayesian approach to constrain the parameter of interest.The definition of the parameter of interest depends on the model to be constrained, as it was pointed out in the previous section.
The most important experimental parameters determining the sensitivity of the experiment are the exposure E, the background rate R bkg , 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 with a subdominant contribution due to other sources R bkg = R 2ν + R others .The systematic uncertainties related to the energy reconstruction can largely differ between experiments and are specific to the detector technology.
A precise evaluation of the sensitivity of an experiment requires considering experiment-specific information.However, it can be approximated by considering 2 Some experiments utilize a multi-variate approach, fitting multiple observable at the same time to better separate signal and background.Here we focus on the simplest approach of onedimensional fit, but the results can be easily generalized.
a counting analysis and Poisson statistics in the region of interest defined above and with a known background expectation given by R bkg .In this derivation, we neglect the systematic uncertainties.Their impact will be discussed later in this section.The dependence of the sensitivity on exposure E and background rate R bkg is different from class I models and class II and II models.
In the following, we will discuss the two cases separately.

B. Experimental sensitivity for class I models
Let's first consider a BSM decay which we classified as a class I model.This is, for example, the case of the existence of an exotic particle, e.g. the Majoron or the Z 2 -odd exotic fermion, with a coupling to neutrino given by g X . 3This is the parameter of interest we want to determine the sensitivity.
The precision with which a subdominant contribution -the number of BSM decays, N X -can be constrained is proportional to the fluctuations of the background in the analysis window, N bkg : The number of background events in the analysis window can be expressed as a function of the background rate and the exposure: and analogously the number of BSM decays: where we have expressed the BSM decay rate as a function of the parameter of interest g X , through the phasespace factor G and the NME M of the decay.Using equations 26 and 27 into equation 25, we obtain: In the last step, we have explicated the contribution of the 2νββ decay to the total background rate (R bkg = R 2ν + R others ).The sensitivity scales with the square root of the exposure, but it is limited by the background, to which the 2νββ decay contributes.In addition, equation 28 shows that uncertainties in the phase space and NMEs can limit the sensitivity.

C. Experimental sensitivity for class II and III models
A slightly different result is obtained for class II and III models.Let's consider, for instance, the double-β decay into sterile neutrinos. 4In this case, the parameter of interest is the mixing angle sin 2 θ.Given that sin 2 θ also modifies the 2νββ decay rate, it is proportional to the ratio between the number of decays into sterile neutrino N νN and the number of 2νββ decay events through the respective phase-space factors.
The statistical uncertainty on this quantity can be computed through standard error propagation Because of the same arguments previously used to define the uncertainty on N X , the uncertainty on N νN will be In addition, we can write N νN in terms of N 2ν and sin 2 θ.Also, the correlation coefficient ρ N νN ,N2ν is proportional to the mixing angle, because of the relation 29.Putting everything together, we can rewrite equation 30 as: where we reintroduced the dependence on the phasespace factors in the last passage.The sensitivity scales with the square root of the exposure, it is limited by the background, but in this case, a high 2νββ decay rate is advantageous due to the dependence σ sin 2 θ ∝ 1/R 2ν .The sensitivity also depends on the ratio between the phase space factors G 2ν /G νN .
. Sensitivity of a double-β decay experiment with 76 Ge to sterile neutrinos with a mass of 500 keV, as a function of the exposure, background level, and systematic uncertainties [78].The results of a full frequentist analysis (markers) are compared to the expectation (lines) from equations 35 and 34.

D. Impact of the systematic uncertainties on the sensitivity
Given an experiment with exposure E and using an isotope with half-life T 2ν 1/2 , we expect the sensitivity to double-β decay into Majorons or Z 2 -odd fermions (or decays resulting from non-standard interaction) to scale with 1/(T 2ν 1/2 • E) while the sensitivity to the sterile neutrino mixing angle (or Lorentz violation and bosonic neutrinos) with T 2ν 1/2 /E.This means that experiments using an isotope with a long 2νββ decay half-life will be favored in the first kind of search as they will have a lower background rate.However, they will be disfavoured in the second kind of search where the sensitivity is linearly proportional to the 2νββ decay half-life.
In our derivation so far, we have neglected the impact of systematic uncertainties on the sensitivity.In a more general way, the sensitivity can be approximated as It is determined by both the statistical uncertainty σ stat , which was derived in equations 28 and 34 as a function of the most important experimental parameters, and the systematic uncertainty σ sys .
As long as the statistical uncertainty is dominant, the sensitivity improves by increasing the exposure approximately as f ∝ 1/E.The sensitivity saturates when the statistical uncertainty becomes comparable with the systematic one.This is illustrated in figure 13.The figure shows the sensitivity of a double-β decay experiment using 76 Ge to sterile neutrinos with a mass of 500 keV, as a function of the exposure, the background level, and the systematic uncertainties [78].The sensitivity computed using a full frequentist analysis is shown by the markers.This is compared to the expectation given by equation 35, where the statistical uncertainty is given by equation 34.

V. DOUBLE-β DECAY EXPERIMENTS AND CONSTRAINS
The 80-years-long history of double-β decay experiments has seen a variety of technologies and concepts being tested and developed over the years with the standing goal of reducing backgrounds and increasing the isotope mass.A pivotal time for the field was around the turn of the century when the unexpected discovery of neutrino masses [121] raised the question of whether that mass could be due to the peculiar mechanism conceived by Majorana [122][123][124], a hypothesis that can be proven by observing 0νββ decay [125].This boosted the interest for double-β decay experiments, setting in motion a process eventually culminating in the consolidation of five main detection technologies: high-purity germanium semiconductor detectors, cryogenic calorimeters, time projection chambers, large liquid scintillators, and tracking calorimeters.The following sections review the most recent experiments related to these technologies and the community's plan for the next-generation projects.We conclude with a summary of the state-ofart constraints on the search for BSM double-β decay in Sec.V G and prospects in Sec.VI.

A. HPGe semiconductor detectors
High-Purity Ge (HPGe) detectors have been a leading technology for double-β experiments since the very first 0νββ decay searches [126].HPGe detectors are semiconductor devices in which electron-hole charge carriers produced by ionization processes are collected by an electric field applied throughout an ultra-pure Ge crystal isotopically enriched in 76 Ge up to 92%.Crystals are grown through the Czochralski method [127] and thus are intrinsically ultra radio-pure.The typical detector size is 1-3 kg, requiring the simultaneous operation of multiple detectors to reach a large target mass.HPGe detectors have superior energy resolution, the best of any doubleβ decay experiment, while also providing information on the event topology.Double-β decays are fully contained within the active detector region, and no volumefiducialization is required to eliminate background leading to very high detection efficiency.
The most sensitive double-β decay searches based on HPGe detectors have been conducted by the Germanium Detector Array (Gerda) experiment [128] and the Majorana Demonstrator [129].Gerda was located at the Laboratori Nazionali del Gran Sasso (LNGS) in cen-tral Italy and operated about 40 kg of HPGe detectors directly immersed in a LAr volume instrumented to detect its scintillation light.The Majorana Demonstrator was located in the Sanford Underground Research Facility (SURF) in South Dakota and operated about 30 kg of HPGe detectors in two vacuum cryostats.The results of the Gerda and Majorana Demonstrator experiments have demonstrated the feasibility of building a ton-scale 76 Ge-based 0νββ decay experiment with ultralow background and superior energy resolution.With the Gerda and Majorana Demonstrator experiments now completed, the next generation experiment will be realized in the framework of the LEGEND project, following two stages named LEGEND-200 and LEGEND-1000 [130].LEGEND-200 has just started approaching physics data taking with 200 kg of HPGe detectors in the upgraded Gerda infrastructure.LEGEND-1000 is currently under preparation and is expected to come online towards the end of the decade.
The Gerda experiment has performed several searches for BSM double-β decays of 76 Ge.Its most sensitive search for Majorons-mediated decays led to half-life constraints of 6.4 • 10 23 yr, 2.9 • 10 23 yr, 1.2 • 10 23 yr, and 1.0 • 10 23 yr (at 90% C.L.), respectively for the decays with spectral index n = 1, 2, 3, and 7 [131].In the same work, limits on Lorentz violation and the decay into light exotic fermions have also been derived.The Lorentz violating isotropic coefficient å(3) of has been constrained to (−2.7 < å(3) of < 6.2)•10 −6 GeV at 90% C.L..In the search for double-β decays into sterile neutrinos, the most stringent limit was obtained for masses between 500 and 600 keV.For these masses, the 90% C.L. interval obtained on the mixing between sterile and active neutrinos is sin 2 θ < 0.013.In the search for double-β decays into Z 2 -odd fermions, limits on the decay half-life have been derived in the range (0.18 − 2.5) • 10 23 yr, at 90% C.L..The best limit of 2.5 • 10 23 yr was obtained for a mass of 400 keV and corresponds to a limit on the coupling constant g χ < (0.6 − 1.4) • 10 −3 MeV −2 , where the range is due to the NME uncertainties.

B. Cryogenic calorimeters
Cryogenic calorimeters, also referred to as bolometers, have been employed for 0νββ decay and dark-matter searches since the 80s [132].Bolometers are crystals coupled with a superconductive thermal sensors.Energy deposition within the crystal produce phonons, which warm up the detectors, and additionally scintillation light for specific materials.To reach the desire sensitive, the thermal sensors typically rely on the resistivity change of superconductive materials at the transition edge, and thus need to be operated at mK temperatures.
The cryogenic calorimeters are extremely versatile tools, as the crystal material can be chosen in order to study a variety of double-β decay isotopes.Further advantages of this detection technique are an excellent en-ergy resolution and very high detection efficiency.The limited crystals dimensions (typically between 0.2 and 0.8 kg) require the operation of thousands of cryogenic calorimeters to reach target masses of hundreds of kilograms or larger.This poses challenging requirements to the cryogenic infrastructure that needs to keep a ton of material stably at mK temperatures for years.
Currently, the largest bolometric experiment is the Cryogenic Underground Observatory for Rare Events (CUORE) at LNGS, which operates about 750 kg of TeO 2 crystals with natural isotopic composition (corresponding to 206 kg of 130 Te) in a large cryogen-free cryostat [133].The CUORE experiment successfully demonstrated the feasibility of a ton-scale bolometric experiment [134], leading to CUPID, the next generation bolometric experiment CUORE Upgrade with Particle Identification enabled by the usage of scintillating bolometers.As part of the R&D towards CUPID, two independent experiments have produced constraints on BSM physics.The first one is CUPID-0, at LNGS, which utilized ZnSe crystals enriched in 82 Se [135].The second one is CUPID-Mo at LSM in France, which utilized Li 2 MoO 4 crystals enriched in 100 Mo [136].

C. Time projection chambers
The first direct observation of 2νββ decay was made using a Time Projection Chamber (TPC) [5].Since then, this technology has been at the forefront of 0νββ decay searches because of the combination of mass scalability and optimal background discrimination capabilities enabled by the 3D reconstruction of the event topology, position, and energy.TPCs are particularly well-suited to search for 0νββ decay of 136 Xe: Xe is a noble element that can be used directly in TPCs as a liquid or gas.On the other hand, TPC detectors have a limited energy resolution and require a multi-variate analysis to constrain background features close to the Q ββ , which are often not resolved (e.g. 214Bi γ line at 2447.7 keV, just below the 136 Xe Q ββ ).
The most sensitive liquid Xe TPC among the current generation of double-β decay experiments was EXO-200 at WIPP near Carlsbad New Mexico.EXO-200 was a single-phase liquid-Xe TPC, filled with 161 kg of 136 Xe [139].EXO-200 demonstrated the capabilities of a monolithic liquid-Xe TPC, which includes relatively good energy resolution, near maximal signal detection efficiency, and solid topological discrimination of backgrounds [140].Building on it, the next generation tonscale 136 Xe-based 0νββ decay experiment nEXO is currently being proposed [141].
High-pressure gaseous Xe TPCs for 2νββ-decay have been developed in the context of the NEXT project, which combines electroluminescence to reach good energy resolution (< 1% FWHM at Q ββ ) and charged-particle tracking for the active suppression of background [142].A concept for a ton-scale phase of NEXT is currently under investigation with the goal of operating a full ton of 136 Xe in the form of enriched Xe gas [143], and possibly deploying technologies for the identification of the 136 Xe decay daughter isotope.

D. Large liquid scintillators
Large liquid scintillator detectors have historically been the most mass-scalable technology used by 0νββ experiments.The scintillator can be doped with different isotopes of interest for these searches, including 136 Xe and 130 Te.Decays occurring within the detector create scintillation photons, which travel straight to the outer surface of the scintillator volume where they are read out by photo-multipliers tubes.The event position and energy release are reconstructed using the time-of-flight of the photons and their number.Multivariate analysis can also provide some information on the initial spatial extension of the energy deposition, giving an extra handle to tag background-like events.
The most important liquid scintillator detector in the field is operated by the KamLAND-Zen experiment in the Kamioka Mine in Japan [146].It consists of a nylon balloon placed at the center of the detector volume and filled with a liquid scintillator in which 136 Xe has been dissolved.After a successful first phase (KamLAND-Zen 400) with up to 340 kg of 136 Xe, the second phase (KamLAND-Zen 800) is currently running with about 680 kg of 136 Xe [147].With the first 1.6 years of data, KamLAND-Zen 800 produced a world-leading 0νββ decay half-life limit [147].The KamLAND-Zen collaboration is already preparing for the ton-scale phase (KamLAND2-Zen) in which about 1 ton of 136 Xe will be deployed [148].
The technology of large liquid scintillators is largely employed in neutrino experiments.Among future large liquid scintillator experiments, we shall mention the SNO experiment.This is a multi-purpose neutrino experiment located at SNOLAB in Canada [149].An upgrade of the experiment is planned after the completion of the main goals, which consists in loading the 780 tons of organic liquid scintillator with a double-beta decaying isotope to search for 0νββ decay [149].

E. Tracking calorimeters
Tracking calorimeters are the only technology capable of measuring with high accuracy the kinematic of electrons emitted in double-β decays, such as the singleelectron energy and the electron angular distribution.Measuring these quantities would give precious inputs to pin down the actual channel mediating the 0νββ decay [151].It would also strongly enhance the sensitivity to search for the BSM double-β decays discussed in this review.
The tracking capability is obtained by decoupling the double-β decay isotope from the detector.The target isotope is placed on a thin foil, immersed in a magnetic field, and surrounded by tracking and calorimetric layers.This configuration enables the measurement of the electron momentum through its bending in the magnetic field, and the measure of its energy when it enters the calorimeters.Unfortunately, it reduces the detection efficiency.The requirement of using very thin foils to minimize energy losses makes it extremely challenging to scale up the isotope mass.
The NEMO-3 experiment utilized this technology to search for 0νββ decay of several isotopes at the LSM in France.Masses from a few grams to a few kilograms of the isotopes of interest were deployed in separate sectors of the detector (6.99 g of 48 Ca [38], 0.932 kg of 82 Se [152], 9.4 g of 96 Zr [153], 6.914 kg of 100 Mo [154], 410 g of 116 Cd [119], and 36.6 g of 150 Nd [30]).A next-generation tracking calorimeter detector is the SuperNEMO Demonstrator, which is based on the technology demonstrated by NEMO-3 [155].In its first phase, the SuperNEMO Demonstrator will deploy one module with 7 kg of 82 Se.A future full-scale experiment is foreseen, consisting of multiple modules aiming for a total 82 Se mass of 100 kg.
The NEMO-3 experiment searched for Majoroninvolving double-β decays in several isotopes:  150 Nd.Limits on the half-life of the decays corresponding to spectral indexes n = 2, 3, and 7 have been derived only with 100 Mo, because of the lower statistics data sets and the higher background achieved with the other isotopes.The corresponding limits at 90% C.L. are: 9.9

F. Other technologies
We should mention the Aurora experiment at LNGS among the experiments using technologies other than those described in the previous subsections.It utilized more than 1 kg of radio-pure cadmium tungstate ( 116 CdWO 4 ) scintillating crystals enriched in the 116 Cd isotope [156].Even if this technology is not competitive in terms of 0νββ decay sensitivities and there is no concrete plan to scale it to a ton-scale 0νββ decay experiments, the Aurora experiment has produced competitive constraints in the search for BSM double-β decays of 116 Cd [35].They set limits on the half-life of Majorons emitting decays with the emission of Majorons at the order of 10 21 yr, while the most stringent limit was obtained for the n = 1 spectral index: T 1/2 > 8.2 • 10 21 yr (at 90% C.L.).In the same work, they searched for Lorentz violation and obtained a limit on the isotropic coefficient of å (3) of < 4.0 • 10 −6 GeV (at 90% C.L.).

G. Most sensitive constraints
In this section we summarise the most sensitive constraints reported by all experiments mentioned in the previous sections, grouping them based on the new physics searched.
a. Double-β decay with the emission of Majorons A summary of the latest results obtained by different double-β decay experiments is presented in table III.The most stringent limits, regardless of the isotope and the experiment, are obtained for the model corresponding to a spectral index n = 1.In fact, the energy distribution predicted for this decay differs the most from the SM 2νββ decay compared to other spectral indexes, as shown in figure 3.Among different experiments, the best limit on the half-life of the Jββ decay (n = 1 mode) Majorons-involving decays.The lower limits on the half-life are converted into upper limits on the neutrino-Majoron coupling constant using equation 7 with the axial vector coupling constant gA = 1.27 and the phase space factors from [64].The NME calculations for the spectral index n = 1 are taken from [7] and references therein, and for n = 3 and n = 7 from [65].
is obtained by the EXO-200 experiment with 136 Xe: 4.3 • 10 24 yr at 90% C.L.. EXO-200 also obtained the best limits on the half-life of the Jββ/JJββ decays (n = 3 modes): 6.3 • 10 23 yr.For the Jββ decay (n = 2 mode), the KamLAND-Zen experiment set the most competitive limit on the half-life at 1.0 • 10 24 yr, while the GERDA experiment set the most competitive limit on the half-life of the JJββ decay (n = 7 mode): 1.0 • 10 23 yr.Recently, EXO-200 has searched for BSM double-β decays in which a Majoron-like particle is emitted via non-standard righthanded currents were also investigated [144].Limits were derived on such a decay for both right-handed and lefthanded hadronic currents: 3.7 of comes from NEMO-3 and is at the order of 10 −7 , a factor 10 better than all the other experiments.This result is attributed to the much larger statistics of 2νββ decay events achieved by the NEMO-3 experiment (∼ 1.9 • 10 5 events in the analysis range).However, part of this difference comes from the statistical treatments as only the CUPID-0 and GERDA experiments used the approach highlighted in this review for class III models (see section IV), to which the Lorentz violating 2νββ decay belongs.Aurora, NEMO-3, and EXO-200 treated the perturbation introduced by Lorentz violation as an independent component in the fit, neglecting any correlation with the SM 2νββ decay distribution in the result.Nevertheless, the impact is hard to quantify and goes beyond the scope of this work.When comparing different results, it should also be considered that the limits on å(3) of depend on the calculated phase space ratio between the SM 2νββ decay and the LV perturbation.The different experiments reviewed in table IV used different calculations.In recent work, improved phase space calculations were performed, in which the Fermi functions are built with exact electron wave functions obtained by numerically solving a Dirac equation in a realistic Coulomb-type potential, including finite nuclear size and screening effects [158].Differences up to 30% for heavier nuclei were found between these improved calculations and the previous calculations using approximated analytical Fermi functions.As pointed out in [158], for example, there is a relevant difference between the newly calculated phase space ratio and the one used by the CUPID-0 collaboration in [137].c. 2νββ decay with bosonic neutrinos The experimental search for an admixture of bosonic and fermionic neutrinos through the search for distortions of the 2νββ decay spectrum has been performed only by the NEMO-3 experiment with 100 Mo [120].They obtained an upper limit on the bosonic neutrino contribution sin 2 χ < 0.27 at 90% C.L.. Also in this case, the statistical treatment does not follow the receipt given in IV for class II models as it neglects the correlation between the SM 2νββ decay (fermionic neutrinos) and 2νββ decay with bosonic neutrinos distributions.As already introduced in III B 2, searches with other isotopes than 100 Mo might be disfavoured by the small predicted ratio r 0 with which the fermionic and bosonic contributions are weighted in the total 2νββ decay rate.
d. Double-β decay into sterile neutrinos and Z 2odd fermions The experimental search for light exotic fermions, i.e. sterile neutrinos and Z 2 -odd fermions, has been performed only by the GERDA experiment with 76 Ge [131], which set a limit on the mixing between active and sterile neutrinos sin 2 θ < 0.013 for a sterile neutrino mass of m N = 500 keV.The limits get worse for lower and higher masses (sin 2 θ < 0.15 for m N = 100 keV, sin 2 θ < 0.050 for m N = 900 keV).GERDA has also set the first direct experimental constraints on the emission of two Z 2 -odd fermions, constraining the half-life of the corresponding decay to 2.5 • 10 23 yr for a mass m χ = 400 keV, which translates into a constraint on the coupling g χ of (0.6 − 1.4) • 10 −3 MeV −2 .Again, limits get worse for smaller and larger masses.

VI. OUTLOOK AND PROSPECTS
Hunting for the extremely rare 0νββ decay, existing experiments collected up to millions of 2νββ decay events.This statistics is expected to rapidly increase, as future experiments, with much larger target mass, will start taking data.We showed in this review that 2νββ decays could be used as powerful probes of new physics.
Since the first experimental searches for double-β decay with the emission of one or two Majorons, we have seen remarkable progress both in the theoretical description of the decays and in the experimental technologies.Improved and more precise calculations of the phase space factors and NMEs are available today, which are essential to convert the experimental constraints on the half-life of the decays into a coupling between the exotic particle, i.e. the Majoron, and neutrinos.On the other hand, the experiments reached incredible precision in the study of 2νββ decay with large statistics data samples and drastic reduction of the background compared to their predecessors, pushing the bounds on the half-life of these decays up to 10 24 yr.Limits on the neutrino-Majoron coupling g J for the Majoron model leading to n = 1 are also available from astrophysics.Supernova observations allow excluding the region of the parameter space 4 • 10 −7 < g J < 2 • 10 −5 by studying the role of Majorons in the Supernova explosion [159,160].Current double-β decay experiments completely excluded the re-gion above the lower Supernova bound.The combined results bring the upper bound down to g J < 4 • 10 −7 , far from the sensitivity of any future double-β decay experiment.Nevertheless, one should remember that Supernova bounds are model-dependent and rely upon additional assumptions.To date, double-β decays provide the best direct constraints on the neutrino-Majoron coupling.
Light exotic fermions can also be searched in doubleβ decays.Depending on the Q ββ of the double-β decay isotope, one or two exotic fermions with a mass between a few hundred keV and a few MeV can be emitted in double-β decay.In the search for sterile neutrino in this mass range, current double-β decay experiments provide bounds that are still weaker than the existing single-β decay bounds, as predicted in [77,78] and confirmed by GERDA results [131].Still, future experiments could reach unexplored regions of the parameter space, down to sin 2 θ ∼ 10 −3 − 10 −4 , for masses between 100 keV and 2000 keV.To date, no experiment exists or is planned with the capability of testing this part of the parameter space (100 keV < m N < 2000 keV). 5Therefore, future double-β decay experiments will provide the best direct constraints on the active-sterile neutrino mixing in the aforementioned mass range.In addition, double-β decay experiments offer a unique opportunity to test all those models in which only the double production of light exotic fermions is allowed, leading to the best direct constraints on the coupling between these exotic fermions and neutrinos of the order of g χ ∼ 10 −4 MeV −2 .
Studying the 2νββ decay spectrum can provide a sensitive test of the so-called counter-shaded Lorentz violation.Current experiments constrained the Lorentzviolating coefficient å(3) of at the level of |å of | < 2.0 • 10 −8 GeV, which is already more competitive of the constraints from double-β decays.Recently, the KATRIN experiment performed a similar analysis using a small set of available data, setting a limit at |å .This limit is ex-pected to further improve up to a sensitivity of 10 −9 GeV or more with the full KATRIN exposure [162].Future double-β decay experiments will be able to improve their current limits (in the best case scenario by a factor of 1/E), nevertheless, hardly reaching single-β decay experiments sensitivity.
A purely bosonic neutrino would substantially change the total double-β decay rate, therefore, the measured 2νββ decay half-life values.Several precision measurements of the 2νββ decay half-life of different isotopes completely ruled out thy hypothesis of a purely bosonic neutrino.On the other hand, experimental data does not completely exclude the hypothesis of a mixed statistic with a partly bosonic neutrino.The only experimental upper limit on the admixture of the bosonic component was set by the NEMO-3 experiment with 100 Mo.The sensitivity to spectral distortions depends on the ratio r 0 between the rates for purely bosonic and purely fermionic neutrinos, which involves phase space factors and NMEs of the two decays.Consequently, some isotopes are favored (i.e. 100 Mo) for future searches of bosonic neutrino admixture compared to others (i.e. 76Ge), for which the very small value of r 0 neutralizes possible effects induced by a partly bosonic neutrino.
Finally, 2νββ decays can be used as a probe of hidden non-standard interaction of neutrinos, like strong neutrino self-interactions and the presence of right-handed currents in weak interactions.Although it was shown that already current experiments could be competitive in constraining such non-standard operators [102,118], no experimental searches have been performed yet.

MFigure 1 .
Figure 1.Mass parabolas of nuclear isobars with even A.Due to the pairing term in the semi-empirical mass formula, β − transitions of even/even nuclei to their odd/odd isobaric neighbor can be energetically forbidden, whereas, in a secondorder process, β − β − decay is allowed.

Figure 2 .
Figure 2. Summed electron energy distribution for the SM 2νββ decay and the lepton number non-conserving 0νββ decay.An infinite energy resolution is assumed, and an arbitrary normalization is used for illustrative purposes.

Figure 3 .
Figure 3. Summed electron energy distribution for different Majoron models (the spectral index corresponding to each model is indicated) compared to the SM 2νββ decay distribution.The decays with the emission of a non-standard Majoron are also shown: they can be triggered by an effective seven-dimension operator, containing right-handed ( φ RR ) and left-handed ( φ RL ) hadronic current.The latter two were adapted from[66].An arbitrary normalization is used for illustrative purposes.

Figure 4 .
Figure 4. (Top) Feynman diagram for the double-β decay with the emission of a Majoron in classical models.(Bottom)Feynman diagram for the emission of a Majoron-like particle φ through an effective dimension-seven operator containing right-handed currents in double-β decay (adapted from[66]).
Figure 5. Summed electron energy distributions of the doubleβ decay into one sterile neutrino with a mass of 600 keV and the double-β decay into two Z2-odd fermions with a mass of 300 keV in comparison to the SM 2νββ decay distribution for the case of76 Ge isotope.An arbitrary normalization is used for illustrative purposes.Adapted from[78]

Figure 6 .
Figure 6.Summed electron energy distribution of the 2νββ decay in the SM and the perturbation term introduced by Lorentz violation (LV).An arbitrary normalization is used for illustrative purposes.

Figure 7 .
Figure 7. Summed electron energy distribution of the 2νββ decay for pure bosonic neutrinos compared to the case of pure fermionic neutrinos (SM 2νββ decay).An arbitrary normalization is used for illustrative purposes.

Figure 8 .
Figure 8. Feynman diagrams of the double-β decay (a) with two left-handed currents, i.e. the SM 2νββ decay, (b) with one right-handed current, and (c) with two right-handed currents.Adapted from [102].

Figure 9 .
Figure 9. Summed electron energy distribution of the 2νββ decay in the presence of right-handed lepton currents compared to the SM 2νββ decay (left-handed lepton currents).An arbitrary normalization is used for illustrative purposes.Adapted from[102].

Figure 11 .
Figure 11.Summed electron energy distribution of 2νSIββ decay where the νSI operator is generated by an s-channel mediator with a mass of M = Q ββ + 0.1me, compared to the SM 2νββ decay.An arbitrary normalization is used for illustrative purposes.Adapted from[118].

Figure 12 .
Figure 12.The three classes of BSM double-β decay searches.

Table III .
Comparison of the results obtained by different double-β decay experiments with different isotopes in the search for

Table IV .
• 10 24 yr and 4.1 • 10 24 yr, respectively.Summary of the results obtained by different double-β decay experiments in the search for Lorentz violation.