e1Corresponding author<EMAIL_ADDRESS>e2Present address: Department of
Physics & Astronomy, University College London, Gower Street, London, WC1E
6BT, UK

11institutetext: Univ Lyon, Université Lyon 1, CNRS/IN2P3, IP2I-Lyon, F-69622,
Villeurbanne, France 22institutetext: National Research Centre “Kurchatov
Institute”, Kurchatov Complex of Theoretical and Experimental Physics, 117218
Moscow, Russia 33institutetext: Dipartimento di Fisica, Sapienza Università di
Roma, P.le Aldo Moro 2, I-00185, Rome, Italy 44institutetext: INFN, Sezione di
Roma, P.le Aldo Moro 2, I-00185, Rome, Italy 55institutetext: Gran Sasso
Science Institute, L’Aquila I-67100, Italy 66institutetext: INFN, Laboratori
Nazionali del Gran Sasso, I-67100 Assergi (AQ), Italy 77institutetext:
Department of Physics, University of California, Berkeley, California 94720,
USA 88institutetext: Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay,
France 99institutetext: Nikolaev Institute of Inorganic Chemistry, 630090
Novosibirsk, Russia 1010institutetext: Dipartimento di Fisica, Università di
Milano-Bicocca, I-20126 Milano, Italy 1111institutetext: INFN, Sezione di
Milano-Bicocca, I-20126 Milano, Italy 1212institutetext: Institute for Nuclear
Research, 03028 Kyiv, Ukraine 1313institutetext: INFN, Sezione di Roma Tor
Vergata, Via della Ricerca Scientifica 1, I-00133 Rome, Italy
1414institutetext: Institute for Astroparticle Physics, Karlsruhe Institute of
Technology, 76021 Karlsruhe, Germany 1515institutetext: IRFU, CEA, Université
Paris-Saclay, F-91191 Gif-sur-Yvette, France 1616institutetext: Nuclear
Science Division, Lawrence Berkeley National Laboratory, Berkeley, California
94720, USA 1717institutetext: Key Laboratory of Nuclear Physics and Ion-beam
Application (MOE), Fudan University, Shanghai 200433, PR China
1818institutetext: Institute for Data Processing and Electronics, Karlsruhe
Institute of Technology, 76021 Karlsruhe, Germany 1919institutetext:
International Center for Advanced Training and Research in Physics (CIFRA),
409, Atomistilor Street, Bucharest-Magurele, 077125, Romania
2020institutetext: Finnish Institute for Educational Research, University of
Jyväskylä, P.O. Box 35, FI-40014 Jyväskylä, Finland 2121institutetext: Center
for Theoretical Physics, Sloane Physics Laboratory, Yale University, New
Haven, Connecticut 06520-8120, USA 2222institutetext: Massachusetts Institute
of Technology, Cambridge, MA 02139, USA 2323institutetext: Physik Department,
Technische Universität München, Garching D-85748, Germany 2424institutetext:
Department of Modern Physics, University of Science and Technology of China,
Hefei 230027, PR China 2525institutetext: Université Grenoble Alpes, CNRS,
Grenoble INP, LPSC/LSM-IN2P3, 73500 Modane, France 2626institutetext:
Laboratory of Nuclear Problems, JINR, 141980 Dubna, Moscow region, Russia
2727institutetext: Faculty of Mathematics, Physics and Informatics, Comenius
University in Bratislava, 842 48 Bratislava, Slovakia 2828institutetext:
Institute of Experimental and Applied Physics, Czech Technical University in
Prague, 128 00 Prague, Czech Republic2929institutetext: Université Grenoble
Alpes, CNRS, Grenoble INP, SIMAP, 38402 Saint Martin d’Héres, France
3030institutetext: John de Laeter Centre for Isotope Research, GPO Box U 1987,
Curtin University, Bentley, Western Australia, Australia

# Probing Beyond the Standard Model physics using the improved description of
100Mo $2\nu\beta\beta$ decay spectral shape with CUPID-Mo

C. AugierIPNL  A. S. BarabashITEP  F. BelliniSapienza,INFN-Roma  G.
BenatoGSSI,LNGS  M. BerettaUCB  L. BergéIJCLab  J. BillardIPNL  Yu. A.
BorovlevNIIC  L. CardaniINFN-Roma  N. CasaliINFN-Roma  A. CazesIPNL  E.
Celie1,GSSI,LNGS  M. ChapellierIJCLab  D. ChiesaMilano,INFN-Milano  I.
DafineiINFN-Roma  F. A. DanevichKINR, INFN-Roma2  M. De JesusIPNL  T.
Dixone2,IJCLab  L. DumoulinIJCLab  K. EitelKIT-IK  F. FerriCEA-IRFU  B. K.
FujikawaLBNLNSD  J. GasconIPNL  L. GironiMilano,INFN-Milano  A. GiulianiIJCLab
V. D. GrigorievaNIIC  M. GrosCEA-IRFU  D. L. HelisLNGS  H. Z. HuangFudan  R.
HuangUCB  L. ImbertIJCLab  A. JuillardIPNL  H. KhalifeCEA-IRFU  M.
KleifgesKIT-IPE  V. V. KobychevKINR  Yu. G. KolomenskyUCB,LBNLNSD  S. I.
KonovalovITEP  J. KotilaJyv,Jyv2,yale  P. LoaizaIJCLab  L. MaFudan  E. P.
MakarovNIIC  P. de MarcillacIJCLab  R. MariamIJCLab  L. MariniLNGS  S.
MarnierosIJCLab  X. F. NavickCEA-IRFU  C. NonesCEA-IRFU  E. B. NormanUCB  E.
OlivieriIJCLab  J. L. OuelletMIT  L. PagnaniniGSSI,LNGS  L. PattavinaLNGS,TUM
B. PaulCEA-IRFU  M. PavanMilano,INFN-Milano  H. PengUSTC  G. PessinaINFN-
Milano  S. PirroLNGS  D. V. PodaIJCLab  O. G. PolischukKINR,INFN-Roma  S.
PozziINFN-Milano  E. PrevitaliMilano,INFN-Milano  Th. RedonIJCLab  A. RojasLSM
S. RozovJINR  V. SanglardIPNL  J. A. ScarpaciIJCLab  B. SchmidtCEA-IRFU  Y.
ShenFudan  V. N. ShlegelNIIC  F. Šimkovic CUB,CTUP  V. SinghUCB  C. TomeiINFN-
Roma  V. I. TretyakKINR,LNGS  V. I. UmatovITEP  L. VagneronIPNL  M.
VelázquezUGA  B. WareCU  B. WelliverUCB  L. WinslowMIT  M. XueUSTC  E.
YakushevJINR  M. ZarytskyyKINR  A. S. ZolotarovaCEA-IRFU

(Received: date / Accepted: date)

###### Abstract

The current experiments searching for neutrinoless double-$\beta$
($0\nu\beta\beta$) decay also collect large statistics of Standard Model
allowed two-neutrino double-$\beta$ ($2\nu\beta\beta$) decay events. These can
be used to search for Beyond Standard Model (BSM) physics via $2\nu\beta\beta$
decay spectral distortions. 100Mo has a natural advantage due to its
relatively short half-life, allowing higher $2\nu\beta\beta$ decay statistics
at equal exposures compared to the other isotopes. We demonstrate the
potential of the dual read-out bolometric technique exploiting a 100Mo
exposure of 1.47 kg $\times$ y, acquired in the CUPID-Mo experiment at the
Modane underground laboratory (France). We set limits on $0\nu\beta\beta$
decays with the emission of one or more Majorons, on $2\nu\beta\beta$ decay
with Lorentz violation, and $2\nu\beta\beta$ decay with a sterile neutrino
emission. In this analysis, we investigate the systematic uncertainty induced
by modeling the $2\nu\beta\beta$ decay spectral shape parameterized through an
improved model, an effect never considered before. This work motivates
searches for BSM processes in the upcoming CUPID experiment, which will
collect the largest amount of $2\nu\beta\beta$ decay events among the next-
generation experiments.

††journal: Eur. Phys. J. C

##  1 Introduction

Neutrinoless double-$\beta$ ($0\nu\beta\beta$) decay is a hypothetical nuclear
decay not allowed by the Standard Model (SM). The potential observation of
$0\nu\beta\beta$ decay would demonstrate the violation of the $B-L$ symmetry
of the SM. Moreover, it would also prove the Majorana nature of neutrinos,
offering crucial insights into the fundamental symmetries governing particle
interactions. Beyond the primary focus on $0\nu\beta\beta$ decay, the pursuit
of large masses, extended data collection periods, and appropriate
$\beta\beta$ isotopes in $0\nu\beta\beta$ experiments lead to a high
collection of events from two-neutrino double-$\beta$ ($2\nu\beta\beta$) decay
Albert  _et al._ (2014); Agostini  _et al._ (2023a); Adams  _et al._ (2021);
Gando  _et al._ (2019); Barabash  _et al._ (2018); Azzolini  _et al._ (2023a);
Arnold  _et al._ (2019); Augier  _et al._ (2023a). $2\nu\beta\beta$ decay is a
second-order weak process occurring in the same $0\nu\beta\beta$ decay
sources. It is particularly interesting for its relevance in the search for
Beyond Standard Model (BSM) processes. Indeed, many theories predict the
existence of exotic double-$\beta$ decays, called this way because they are
characterized by a continuum energy distribution of the two emitted electrons
similar to the $2\nu\beta\beta$ one (Fig. 1). These include the emission of
new exotic particles such as scalar bosons known as “Majorons” or massive
sterile neutrinos. Additionally, there is the possibility of observing
violations of fundamental symmetries like the Lorentz invariance.

Majorons are massless Nambu-Goldstone bosons resulting from the spontaneous
$B-L$ symmetry breaking in the low-energy regime Gelmini and Roncadelli
(1981); Chikashige  _et al._ (1980) and could play a role in the
$0\nu\beta\beta$ decay coupling to the Majorana neutrinos Georgi  _et al._
(1981). Despite many theories being disfavored by the accurate measurements of
the width of the Z boson decay Schael  _et al._ (2006), currently, different
models predict the emission of one or two Bamert  _et al._ (1995) Majorons in
the $0\nu\beta\beta$ decay, these are denoted as

$(A,Z)\rightarrow(A,Z+2)+2e^{-}+\chi_{0}\quad(\text{$\beta\beta\chi_{0}$}),$
(1)

or

$(A,Z)\rightarrow(A,Z+2)+2e^{-}+2\chi_{0}\quad(\text{$\beta\beta\chi_{0}\chi_{0}$}).$
(2)

In these models, the Majoron carries a non-zero lepton number Burgess and
Cline (1994) or it is a component of a massive gauge boson Carone (1993) or a
“bulk” field Mohapatra  _et al._ (2000). Recently, schemes with massive
Majoron-like particles have become popular since they could play the role of
dark matter Berezinsky and Valle (1993); Heeck (2017); Brune and Päs (2019);
Cepedello  _et al._ (2019). The signature of these decays can be distinguished
from the SM $2\nu\beta\beta$ one thanks to the different spectral index $n$,
which determines the position of the maximum intensity in the two emitted
electrons’ energy spectrum. In particular, the differential decay rate is
proportional to

$\frac{d\Gamma}{dT}\propto\left(Q_{\beta\beta}-T\right)^{n},$ (3)

where $T$ is the total kinetic energy of the two electrons emitted, and
$Q_{\beta\beta}$ is the Q-value. The value of the spectral index is $n=5$ for
the SM $2\nu\beta\beta$ decay, while for the Majoron-emitting modes, it can be
$n=1,2,3$, or 7 (left panel Fig. 1).

Some BSM theories involve a Lorentz invariance Violation (LV) and the CPT
(Charge-Parity-Time reversal) symmetry violation terms in the Lagrangian.
These theories have been developed in the Standard Model Extension (SME)
framework, in such a way that the SM gauge invariance is preserved Colladay
and Kostelecky (1997, 1998). The neutrino sector was extensively studied in
the SME framework Kostelecky and Mewes (2004a, b, 2012) and the majority of
these effects are experimentally investigated through neutrino oscillations
and time-of-flight measurements Kostelecky and Russell (2011). Nevertheless,
four operators equally change all neutrino energies and have no impact on
oscillations, which can be studied through weak decays Kostelecky and Tasson
(2009); Díaz  _et al._ (2013) and are the so-called countershaded operators,
labeled as “oscillation free” (of). The interaction of these operators with
neutrinos modifies their four-momentum in a way that

$q^{\alpha}=(\omega,\mathbf{q})\longrightarrow
q^{\alpha}=\left(\omega,\mathbf{q}+\mathbf{a}_{of}^{(3)}-\dot{a}_{of}^{(3)}\hat{\mathbf{q}}\right),$
(4)

where $\mathbf{a}_{of}^{(3)}$ encodes all the 3 directional of components
while $\dot{a}_{of}^{(3)}$ represents the anisotropic one Kostelecky and
Tasson (2009); Díaz  _et al._ (2013). This results in a perturbation of the
$2\nu\beta\beta$ decay rate, which can be written as the sum of two components

$\Gamma=\Gamma_{SM}+10\dot{a}_{of}^{(3)}\Gamma_{LV},$ (5)

where $\Gamma_{SM}$ is the SM $2\nu\beta\beta$ decay rate and $\Gamma_{LV}$ is
the perturbation induced by the Lorentz violation. The energy distribution of
the two electrons emitted in the LV $2\nu\beta\beta$ decay is shifted to
higher energies compared to the SM $2\nu\beta\beta$ decay (central panel Fig.
1).

Figure 1: Standard Model $2\nu\beta\beta$ decay (black) compared with
different exotic double-$\beta$ decay spectra. Left panel: comparison of
different Majoron-emitting modes with spectral indexes $n=1,2,3$, and 7.
Central panel: Lorentz-violating $2\nu\beta\beta$ decay represented in the
case of pure LV $2\nu\beta\beta$ decay (dashed line) and summed to the SM
$2\nu\beta\beta$ one with an arbitrary value of $|$$\dot{a}_{of}^{(3)}$ $|$
(solid line). Right panel: $2\nu\beta\beta$ decay with sterile neutrino
emission ($\nu N\beta\beta$) for $m_{N}$ = 0.5 MeV (blue) and 1.5 MeV (green)
in case of pure $\nu N\beta\beta$ decay (dotted line) and mixed with the SM
$2\nu\beta\beta$ one assuming $\sin^{2}\theta=0.1$ (solid line).

In the past years, some BSM theories have hypothesized the existence of a
sterile neutrino $N$ with a mass $m_{N}$ at accessible energies Chattopadhyay
and Patel (2017); Lee  _et al._ (2013); Shaposhnikov (2007); Deppisch and
Pilaftsis (2011); Kersten and Smirnov (2007); Pilaftsis (2005); Bolton  _et
al._ (2020). The sterile neutrino can also be considered a candidate for dark
matter Dasgupta and Kopp (2021). Experimentally, the parameter of interest is
the active-sterile mixing strength $\sin^{2}\theta$, which determines the
mixing angle between the electron neutrino flavor and the sterile one. The
existence of a sterile neutrino with a mass $m_{N}<Q_{\beta\beta}$ can induce
an effect on the $2\nu\beta\beta$ decay spectral shape Bolton  _et al._
(2021); Agostini  _et al._ (2021), implying the emission of one sterile
neutrino in the decay

$(A,Z)\rightarrow(A,Z+2)+2e^{-}+\bar{\nu}+N\quad(\nu N\beta\beta),$ (6)

while the emission of two sterile neutrinos is considered negligible. The
$0\nu\beta\beta$ experiments are sensitive to the sterile neutrino mass range
from 0.1 MeV up to 3 MeV (depending on $Q_{\beta\beta}$) due to kinematical
conditions. This region is particularly interesting due to the relative
weakness of the actual constraints coming from single $\beta$-decay
experiments, $\sin\theta\sim$ 10-3–10-2 Holzschuh  _et al._ (2000);
Schreckenbach  _et al._ (1983); Deutsch  _et al._ (1990); Bellini  _et al._
(2013); Derbin  _et al._ (2018). The effect on the total decay rate is given
by

$\Gamma=\cos^{4}\theta\Gamma_{SM}+2\cos^{2}\theta\sin^{2}\theta\Gamma_{\nu
N},$ (7)

where the first term $\Gamma_{SM}$ accounts for the SM $2\nu\beta\beta$ decay
and the second term $\Gamma_{\nu N}$ represents the $\nu N\beta\beta$ decay
Bolton  _et al._ (2021); Agostini  _et al._ (2021). The energy spectrum in the
case of pure $\nu N\beta\beta$ decay events is characterized by a shift in the
Q-value determined by the mass $m_{N}$ of the sterile state and a smaller
Phase Space Factor (PSF) compared to $2\nu\beta\beta$ decay (right panel Fig.
1).

In this framework, scintillating cryogenic calorimeters are one of the most
promising technologies for $0\nu\beta\beta$ decay searches Augier  _et al._
(2022); Azzolini  _et al._ (2022) and the study of $2\nu\beta\beta$ decay
spectral shape Augier  _et al._ (2023a); Azzolini  _et al._ (2023a). These
detectors offer outstanding capabilities in terms of energy resolution,
radiopurity, background rejection, detection efficiency, and mass scalability
Adams  _et al._ (2022); Augier  _et al._ (2022, 2023b); Azzolini  _et al._
(2022). CUPID-Mo exploits this technology to demonstrate the performances and
the experimental feasibility of a mid-scale experiment with Li2MoO4-based
bolometers. In this paper, we present the results of the search for exotic
double-$\beta$ decays of 100Mo exploiting the full data taking of the CUPID-Mo
experiment. This analysis relies on the background model which provides a
detailed description of the background sources releasing energy in the CUPID-
Mo detector Augier  _et al._ (2023b, a). Differently from the previous
analyses, CUPID-Mo is the first experiment to account for the systematic
uncertainties of the $2\nu\beta\beta$ decay spectral shape, parameterized
through the improved $2\nu\beta\beta$ decay description Šimkovic  _et al._
(2018); Niţescu  _et al._ (2021), in the search for exotic double-$\beta$
decays. The promising results obtained in this analysis motivate the interest
in investigating the CUPID potential in the search for new physics processes
involving distortion of the $2\nu\beta\beta$ decay spectral shape. CUPID is a
next-generation ton-scale experiment aiming to reach unprecedented
sensitivities on $0\nu\beta\beta$ decay of 100Mo using Li2100MoO4 cryogenic
calorimeters Armstrong  _et al._ (2019). After one year of data taking, the
100Mo exposure in CUPID will be about 150 times higher than the CUPID-Mo,
providing impressive statistics on the SM $2\nu\beta\beta$ decay. Finally,
with this work we analyze the main limits in searching for these BSM processes
with cryogenic calorimeters and propose possible solutions to overcome these
problems in future searches with CUPID.

##  2 Experimental setup

CUPID-Mo is an experiment at Laboratoire Souterrain de Modane (LSM) in France
aiming to search for $0\nu\beta\beta$ decay of 100Mo. It was installed in the
EDELWEISS cryostat Armengaud  _et al._ (2017), optimized for low background
measurements. The detector acquired data at a stable temperature of $\sim$20
mK from March 2019 to June 2020, for a total 100Mo exposure of 1.47 kg
$\times$ y. The experiment is made up of 20 cryogenic calorimeters consisting
of Li2100MoO4 crystals produced from molybdenum enriched in 100Mo to $\sim$97%
and read-out by Neutron Transmutation Doped (NTD) Ge thermistors. Each crystal
has a cylindrical shape and an average weight of 210 g. Li2MoO4-based
bolometers are among the best detectors in this field due to the excellent
energy resolution (7.7 $\pm$ 0.4 keV FWHM at 3034 keV Augier  _et al._ (2022);
Poda and Giuliani (2017)) and the intrinsic radio-purity Augier  _et al._
(2023b). Moreover, these scintillating crystals permit the identification of
the $\alpha$ background from the $\beta/\gamma$ interactions Augier  _et al._
(2022). To detect the scintillation light and allow particle identification,
thin Ge wafers, which also operate as cryogenic calorimeters, are positioned
between Li2100MoO4 crystals along the tower. These Light Detectors (LDs) are
also read out with NTD Ge thermistors and coated with a $\sim$70 nm layer of
SiO2 to increase the light collection. In this setup, each Li2100MoO4 crystal
faces the top and the bottom with LDs, as an exception for the crystals on the
top floor which face only one LD. A cylindrical copper holder and PTFE
(polytetrafluoroethylene) pieces constitute the supporting structure for each
crystal and the adjacent LDs. In addition, a reflective foil (3M Vikuiti™)
surrounds each Li2100MoO4 crystal to increase the light collection. The
detector set-up consists of 20 bolometers arranged in five towers, four
modules for each tower. These are suspended using stainless steel springs to
reduce the vibrational noise. More details on the detector structure are in
Ref. Armengaud  _et al._ (2020). The cryogenic set-up is composed of five
copper screens corresponding to the different thermal stages (300 K, 100 K, 50
K, 1 K, and 10 mK). Two lead shields aim to screen the detector from the
environmental $\gamma$ radioactivity, an internal 14 cm thick Roman lead
shield installed at the 1 K stage and a 20 cm thick external lead shield
Armengaud  _et al._ (2013). In the same way, two 10 cm and 55 cm thick
polyethylene shields (internal and external, respectively) are installed to
shield against environmental neutrons Armengaud  _et al._ (2020). Finally, the
entire setup is surrounded by plastic scintillators acting as muon veto.

##  3 Experimental data

Experimental data are acquired as a continuum time stream and digitized with a
sampling frequency of 500 Hz both for Li2100MoO4 detectors and LDs. The data
are divided into twelve datasets, where each dataset corresponds to about 1
month of data taking and is sub-divided into a series of runs. Each run is
characterized by a period of $\sim$24 hours of stable data taking. At the
beginning and the end of each dataset, specific calibration runs are acquired
by deploying a 232Th/238U source in the vicinity of the detector. The
characteristic $\gamma$-lines of 232Th and 238U produce several peaks in
Li2100MoO4 crystals enabling the calibration of all the detectors. The light
detectors cannot be calibrated with the 232Th/238U sources since $\gamma$-rays
are not fully contained in the LD volume. For this reason, specific runs are
acquired by using an intense 60Co source, able to stimulate the production of
fluorescence 100Mo X-rays at $\sim$17 keV used to calibrate LDs.

Experimental data are acquired and processed with C++ based software packages
developed by previous bolometric experiments Di Domizio  _et al._ (2018);
Armengaud  _et al._ (2020). All the triggered events are acquired in a
3-seconds time window (1 second of pre-trigger). Each Li2100MoO4 detector is
associated with its adjacent LD(s). First, we apply a series of cuts aiming to
remove periods of detector instabilities and not-optimized data takings. We
reconstruct the pulse amplitudes using an optimal filter Gatti and Manfredi
(1986). This method allows the production of a new pulse with the maximum
signal-to-noise ratio by selectively weighting the frequency components of the
signal and suppressing those that are more affected by noise. Then, the
filtered waveforms are corrected in order to remove the amplitude’s dependence
on the initial detector temperature. For more details on the data processing
see Ref. Augier  _et al._ (2022). The detectors’ energy scale is determined
through calibration runs, which provide the calibration function parameters
calculated using the $\gamma$-ray peaks produced by the 232Th/238U sources. A
quadratic function with a zero intercept is used to calibrate the
$\beta/\gamma$ energy region. The $\alpha$ events appear in the background
spectrum in a high energy range that spans from 4 to 10 MeV. Due to the
different detector responses to $\gamma$-rays and $\alpha$-particles, the
$\alpha$ spectrum shows a mis-calibration of about 8%. In this case, we use
the 210Po $\alpha$-peak to re-calibrate the $\alpha$ region. Given the
granularity of the detector, the information about the timing of each event
with respect to the other detectors is extremely useful for studying the
topology of radioactive decays, as explained in the next section. In order to
do that, we tag each event occurring in a $\pm$10 ms time window with other
detectors. The variable describing the number of channels that triggered a
pulse in the same time window is called multiplicity.

##  4 Data selection

The data selection for the source spectra employed in the background model fit
(explained in Sec. 6) is based on particle identification and time-coincidence
criteria. The $\alpha$ identification relies on the scintillation light
detected by the adjacent LD(s) coupled to each Li2100MoO4 crystal. For each
event, we calculate an energy-independent variable that takes into account the
amount of scintillation light detected in both LD(s) to distinguish the
$\alpha$ events from the $\beta/\gamma$ ones. The cut on the normalized light
signals allows the identification of the $\beta/\gamma$ events with an
efficiency of $(99.4\pm 0.4)\%$ Augier  _et al._ (2022).

Following the approach used in other bolometric experiments Alduino  _et al._
(2017); Azzolini  _et al._ (2023a), we select the following spectra:

  * •

$\mathcal{M}_{1\alpha}$: $\alpha$ events with multiplicity = 1 and an energy
in the 3 to 10 MeV range.

  * •

$\mathcal{M}_{1\beta/\gamma}$: $\beta/\gamma$ events with multiplicity = 1 and
an energy in the 0.1 to 4 MeV range.

  * •

$\mathcal{M}_{2}$: the sum of energies released in two crystals coincident in
time (multiplicity = 2) within the 0.2–4 MeV range.

To select a clean data sample, we apply a series of quality cuts based on
pulse shape and light yield. The cut on the pulse shape exploits the principal
component analysis Augier  _et al._ (2022). In order to provide an almost
clean sample of signal-like events, we exploit $2\nu\beta\beta$ decay events
between 1–2 MeV. We also require that the acquired waveforms contain only one
pulse in the 3-second window and they satisfy the anti-coincidence with the
muon veto.

The efficiencies for the data quality cuts are evaluated using the
$\gamma$-peaks from $\mathcal{M}_{1\beta/\gamma}$ and $\mathcal{M}_{2}$
spectra. The efficiency of the muon veto and multiplicity selection is
assessed using the 210Po peak. Finally, we use random noise triggers to
evaluate the pile-up efficiency, which represents the probability of a pulse
being superimposed with another pulse within a 3-second window. Further
details on the data selection and efficiency evaluation can be found in Refs.
Augier  _et al._ (2023b, 2022). The resulting efficiencies from the cuts
applied in CUPID-Mo data are $\varepsilon_{1}=(88.9\pm 1.1)$% for
$\mathcal{M}_{1\beta/\gamma}$, $\varepsilon_{2}=(83.3\pm 2.5)$% for
$\mathcal{M}_{2}$, and $\varepsilon_{3}=(94.7\pm 1.0)$% for
$\mathcal{M}_{1\alpha}$.

##  5 Monte Carlo simulations

We reproduce the signature of the background sources by producing a series of
Monte Carlo simulations. For this purpose, we employ the version 10.04 of
Geant4 Agostinelli  _et al._ (2003). The geometry implemented in simulations
faithfully reproduces the experimental structure from small detector
components to cryostat vessels and radiation shields (see Ref. Augier  _et
al._ (2023b) for more details). We generate radioactive decays in the
experiment components using both Decay0 Ponkratenko  _et al._ (2000) and
Geant4 (G4RadioactiveDecay library). Particle propagation through the
experimental geometry employs the Livermore low-energy physics models liv
(2023).

We use the detector response model to reproduce the simulated spectra with the
measured data. We model the energy resolution with a Gaussian shape where the
mean is the energy deposited in the simulation and the standard deviation is
derived from the experimental data. Pulses falling below the energy threshold
of $<$40 keV are excluded, and the multiplicity is reproduced within a $\pm$10
ms time window. We parameterize the scintillation light energy measured by the
LD as a second-order polynomial of Li2100MoO4 detectors energy to reconstruct
light signals in data. The LD energy resolution is modeled as a Gaussian with
a standard deviation depending on the energy deposited in the corresponding
Li2100MoO4 crystal. This approach ensures faithful replication of the light
yield cuts in the simulations as observed in the experimental data. Data
selection efficiencies are also taken into account in the simulations by
generating a random number $r$ for each event uniformly distributed between 0
and 1.

##  6 Background model

In our previous work, we conducted a thorough background model analysis
utilizing a Markov Chain Monte Carlo (MCMC) Bayesian fit sampling the joint
posterior probability density function (p.d.f.) of the model parameters
utilizing the JAGS software Plummer  _et al._ (2003); Augier  _et al._
(2023b). In the list of background components utilized in the fit, we included
232Th and 238U contaminations in all detector and cryostat components. We
account for the breaks in the secular equilibrium of 232Th and 238U sources by
producing separate simulations of their sub-chains. Specifically, in the 232Th
chain, break points occur at 228Ra and 228Th, while in the 238U chain, break
points are 234U, 230Th, 226Ra, and 210Pb. Additionally, we consider other
contributions, such as 40K contamination in the springs and the outermost
cryogenic thermal shield, 60Co from cosmogenic activation in all copper
components, and ${40}$K, 87Rb and 90Sr/90Y in crystals. Decays are generated
within the bulk of components and on the surface for nearby elements,
following an exponential density profile $e^{-x/\lambda}$, where $\lambda$ is
a variable depth parameter (for most of the surface contaminants $\lambda$ =
10 nm, see Ref. Augier  _et al._ (2023b) for details). Initially, the
$2\nu\beta\beta$ decay of 100Mo was simulated in the background model under
the Single-State Dominance (SSD) hypothesis, utilizing exact Dirac electron
wave functions Kotila and Iachello (2012). Subsequently, we incorporated an
improved description Šimkovic  _et al._ (2018); Niţescu  _et al._ (2021) of
the $2\nu\beta\beta$ decay into the fit, with marginalization over the
theoretical uncertainty of the spectral shape Augier  _et al._ (2023a). The
fit also includes the $2\nu\beta\beta$ decay of 100Mo to 100Ru 0${}_{1}^{+}$
excited state, pile-up produced by $2\nu\beta\beta$ decay events occurring in
the same crystal, and random coincidences between two crystals.

Finally, a total of 67 sources are included in the fit. We modeled all the
priors as non-negative uniform p.d.fs. The only informative priors are set for
the $2\nu\beta\beta$ decay from 100Mo to 100Ru 0${}_{1}^{+}$ excited state
(with a half-life of $T_{1/2}=(6.7\pm 0.5)\times 10^{20}$ y Barabash (2020)),
the stainless steel springs contamination (from screening measurements), the
accidental coincidences, determined from the rate of single events, and the
$2\nu\beta\beta$ decay pile-up, estimated with the calibration runs (see Ref.
Augier  _et al._ (2023b) for details). Variable binning is employed for the
three spectra to ensure sufficient counts in each bin, thereby minimizing the
impact of statistical fluctuations. A minimum bin size of 15 keV is chosen for
$\mathcal{M}_{1\beta/\gamma}$ and $\mathcal{M}_{2}$, while a 20 keV minimum
bin size is used for $\mathcal{M}_{1\alpha}$. The minimum number of counts in
each bin is 50 for $\mathcal{M}_{1\beta/\gamma}$ and 30 for $\mathcal{M}_{2}$
and $\mathcal{M}_{1\alpha}$. Additionally, each peak is selected to be fully
contained within one bin to mitigate the systematic effects of the detector
response on the results.

##  7 BSM analysis

The search for exotic double-$\beta$ decays is performed by individually
incorporating the new physics spectra into the background model fit. Each MCMC
sample reconstructs the joint posterior p.d.f. of the normalization
parameters. The normalization parameter $N_{j}$ of the j-th BSM process is
directly related to the decay rate through the formula

$\Gamma_{j}=\dfrac{N^{j}_{MC}\cdot N_{j}}{T\cdot\varepsilon\cdot
N_{\beta\beta}}.$ (8)

Here, $N^{j}_{MC}$ denotes the total number of decays generated in the
simulation, $N_{\beta\beta}$ is the total number of 100Mo atoms, $T$ is the
experiment lifetime and $\varepsilon$ is the efficiency. In our case, the
efficiency has been included in the simulations and, therefore, already
accounted into the normalization parameter. The value of the product $T\cdot
N_{\beta\beta}$ in CUPID-Mo corresponds to $(99.7\pm 0.2)\times 10^{23}$ 100Mo
emitters $\times$ y. All new physics spectra in this analysis are simulated
using exact Dirac wave functions with finite nuclear size and electron
screening Kotila and Gráf (2023). We use the SSD approximation for all the BSM
spectra since the difference with their corresponding higher-state dominance
(HSD) approximated spectrum in the total number of events is negligible
compared to the systematic uncertainties of the model. We define as reference
fit the configuration of sources described above. We consider two possible
scenarios: the “standard” background model fit, as described in Ref. Augier
_et al._ (2023b), which uses the SSD hypothesis for 100Mo $2\nu\beta\beta$
decay, and the Improved Model fit (IM), detailed in Ref. Augier  _et al._
(2023a), which utilizes the improved description for the 100Mo
$2\nu\beta\beta$ decay spectral shape Šimkovic  _et al._ (2018); Niţescu  _et
al._ (2021). Previous analyses from many experiments in this field use either
the SSD or the HSD approximation to model the $2\nu\beta\beta$ decay spectral
shape Agostini  _et al._ (2022); Azzolini  _et al._ (2023b); Arnold  _et al._
(2015, 2019); Barabash  _et al._ (2018); Kharusi  _et al._ (2021); Gando  _et
al._ (2012). For the first time in the search for exotic double-$\beta$
decays, the uncertainties of the SM $2\nu\beta\beta$ decay spectral shape are
marginalized in the fit. Following the improved description, the
$2\nu\beta\beta$ decay rate or, more specifically, the PSF can be decomposed
as

$G^{2\nu}=G_{0}^{2\nu}+\xi_{31}^{2\nu}G_{2}^{2\nu}+\frac{1}{3}\left(\xi_{31}^{2\nu}\right)^{2}G_{22}^{2\nu}+\left[\frac{1}{3}\left(\xi_{31}^{2\nu}\right)^{2}+\xi_{51}^{2\nu}\right]G_{4}^{2\nu},$
(9)

where $G_{0},G_{2},G_{22}$ and $G_{4}$ are the PSFs for different terms in the
Taylor expansion of the lepton energies and $\xi_{31}$, $\xi_{51}$ are
parameters depending on the ratios of the Gamow-Teller Nuclear Matrix Elements
(NMEs). For 100Mo, the values of the PSFs are $G_{0}$ = 3.303 $\times$
$10^{-18}$ yr-1, $G_{2}$ = 1.509 $\times$ $10^{-18}$ yr-1, $G_{22}$ = 4.320
$\times$ $10^{-19}$ yr-1, and $G_{4}$ = 1.986 $\times$ $10^{-19}$ yr-1
Šimkovic  _et al._ (2018). In the IM fit, the four spectra corresponding to
$G_{0},G_{2},G_{22}$, and $G_{4}$ are simulated using exact Dirac wave
functions with finite nuclear size and electron screening Kotila and Gráf
(2023), and they are included in the fit separately. To accurately describe
the $2\nu\beta\beta$ decay spectral shape, the fit model is modified to
marginalize over the $\xi$ parameters. Since $\xi_{31}$ and $\xi_{51}$ are
strongly anti-correlated, a Gaussian prior is placed on the ratio
$\xi_{31}/\xi_{51}$ with a mean equal to the SSD prediction and a conservative
5% uncertainty Augier  _et al._ (2023a). This choice relies on the nuclear
structure calculations, where the value of $\xi_{31}/\xi_{51}$ can be reliably
obtained Šimkovic  _et al._ (2018); Coraggio  _et al._ (2022). Within the SSD
hypothesis, the value of the ratio is 0.367 Šimkovic  _et al._ (2018). For
more details on the IM fit, see Ref. Augier  _et al._ (2023a).

We perform a series of additional fits to assess the systematic uncertainties.
These are listed as follows:

  * •

The dominant contribution at low energies comes from 232Th and 238U
contaminations in cryostat components. However, the signatures of these
sources in close (10 mK) and far (300 K) components are almost degenerate,
leading to a possible mis-modeling. To estimate the uncertainty related to the
source location, we alternatively remove the 300 K and 10 mK sources from the
fit.

  * •

In the background model fit, the decay chain 90Sr/90Y is included. Both these
isotopes decay through pure $\beta$-decays, producing a featureless spectrum
correlating with the $2\nu\beta\beta$ decay and BSM components. In the
background model fit we measure a 90Sr/90Y activity of $179^{+36}_{-32}$
$\mu$Bq/kg. However, the presence of 90Sr/90Y is still uncertain since other
unexplained contributions at low energies can induce its convergence. This
systematic test involves removing this contribution from the fit.

  * •

In the reference fit, some contributions present a posterior p.d.f. converging
to a value compatible with zero or show an exponential shape flattened to
zero. The minimal model is a fit performed by removing these contributions.

  * •

As described above, we use a minimum bin size of 15 keV on
$\mathcal{M}_{1\beta/\gamma}$, nevertheless, we consider also different values
for the binning (1 keV, 2 keV, and 20 keV).

  * •

Two fits are performed varying the energy scale by $\pm$1 keV to account for a
possible energy bias.

  * •

The theoretical uncertainties on the $2\nu\beta\beta$ decay Bremsstrahlung
cross-section may affect the accuracy of MC simulations on the
$2\nu\beta\beta$ decay spectrum. In order to assess this uncertainty, we
perform the fit with alternatives $2\nu\beta\beta$ decay spectra obtained by
varying the Bremsstrahlung cross-section by $\pm$10% Pandola  _et al._ (2015).

The uncertainties on the efficiency (1.2%) and the 100Mo enrichment (0.2%) are
directly marginalized in the posterior p.d.fs. with Gaussian priors. In the
following sections, the way to extract the physical parameters for the
different BSM processes and their results are described.

##  8 Results

Even though the search for all the exotic double-$\beta$ decays is conducted
with the same fitting procedure, the way to extract a limit on the physical
parameters is significantly different. Majoron-emitting decays do not
interfere with the SM $2\nu\beta\beta$ decay which is just considered an
independent process. In this case, the SM $2\nu\beta\beta$ decay acts as a
background. Conversely, the existence of $\nu N\beta\beta$ decay tends to
suppress the SM $2\nu\beta\beta$ decay rate by a factor $\cos^{4}\theta$. The
Lorentz-violating $2\nu\beta\beta$ decay also interferes with the SM
$2\nu\beta\beta$ decay rate by introducing an additional term on the PSF. In
the last two scenarios, increasing the sample size of $2\nu\beta\beta$ decay
events will result in lower statistical uncertainties. This will make it more
sensitive to deviations. 100Mo is one of the best isotopes among the
$2\nu\beta\beta$ decay emitters to detect Lorentz-violating $2\nu\beta\beta$
decay and sterile neutrino emissions because of its relatively short half-
life.

###  8.1 Majoron-emitting decays

In this study, we performed the analysis of four different Majoron-emitting
modes, corresponding to $n=1,2,3$, and 7. The parameter of interest is the
decay rate $\Gamma_{0\nu M}$ of the different processes, which can be easily
converted into a limit on the neutrino-Majoron coupling $g_{ee}^{M}$ through
the formula

$\Gamma_{0\nu M}=G_{0\nu M}\left|\left\langle
g_{ee}^{M}\right\rangle\right|^{2m}\left|M_{0\nu M}\right|^{2},$ (10)

where $G_{0\nu M}$ and $\left|M_{0\nu M}\right|$ are PSFs and NMEs for the
Majoron-emitting modes, while $m$ is the number of Majorons emitted. The 100Mo
values of $G_{0\nu M}$ and $\left|M_{0\nu M}\right|$ for different Majoron-
emitting modes are summarized in Table 1.

Table 1: 100Mo values of the PSF $G$ and NME for different Majoron-emitting decays. The values of $G$ are taken from Ref. Kotila  _et al._ (2015). For the Majoron mode with $n$ = 1, the NME values are taken from Ref. Agostini  _et al._ (2023b). For the other modes, the NMEs have been calculated in the framework of the interacting boson model Kotila and Iachello (2021). decay mode |  $G$ [$\times 10^{-18}$ y-1] | NME
---|---|---
$\beta\beta\chi_{0}$ ($n$ = 1) | 598 | 3.84–6.59
$\beta\beta\chi_{0}$ ($n$ = 2) | - | -
$\beta\beta\chi_{0}$ ($n$ = 3) | 2.42 | 0.263
$\beta\beta\chi_{0}\chi_{0}$ ($n$ = 3) | 6.15 | 0.0019
$\beta\beta\chi_{0}\chi_{0}$ ($n$ = 7) | 50.8 | 0.0019

No evidence of signal was found for any of the decays mentioned above setting
90% Credible Interval (CI) limits on their half-lives by integrating the
posterior p.d.f. on the corresponding decay rates. We performed the
aforementioned systematic tests to determine the systematic uncertainty. For
$n=1,2$, and 3, the systematics with a greater impact on the results are +10%
Bremsstrahlung, –1 keV energy scale, and 90Sr/90Y. In the Majoron mode $n=7$,
where the signature exhibits a spectrum shifted at lower energies compared to
the SM $2\nu\beta\beta$ decay (see Fig. 1), the situation is inverted, with
dominant effects from –10% Bremsstrahlung, +1 keV energy scale, and source
location. These tests demonstrate that uncertainties in modeling the low-
energy part of the $\mathcal{M}_{1\beta/\gamma}$ spectrum limit the
sensitivity for exotic double-$\beta$ decay searches. In a conservative
approach, we quote as the final limit the least stringent including
systematics. The same systematic checks are performed in the IM fit. The
uncertainty on the $2\nu\beta\beta$ decay spectral shape, parameterized in the
improved description, significantly reduces the sensitivity for Majoron-
emitting decays, as depicted in Fig. 2. This effect is more pronounced in
Majoron modes with $n<5$, where the signal relies mostly on the high-energy
side of the $2\nu\beta\beta$ decay spectrum. We report the limits on half-
lives and neutrino-Majoron coupling constants obtained in this analysis for
SSD and IM in Table 2. Only results with the SSD fit can be directly compared
with the other experiments, where the same assumption on the $2\nu\beta\beta$
decay spectral shape was made. In the SSD assumption, the obtained limits are
less stringent with factors of 1.8 for $n$ = 1, 1.7 for $n$ = 2, 2 for $n$ =
3, and 5.4 for $n$ = 7 compared to NEMO-3 Arnold  _et al._ (2015, 2019).
However, the 100Mo exposure available in CUPID-Mo ($\sim$1.5 kg $\times$ y) is
22 times less than the NEMO-3 exposure ($\sim$34 kg $\times$ y), demonstrating
the high sensitivity of the dual readout bolometric technique for these
searches.

Figure 2: Upper limits at 90% CI on the decay rate of different Majoron-emitting decays with the SSD approximation for $2\nu\beta\beta$ decay (green) and the improved description implemented in the fit (orange). The solid color corresponds to the limit obtained from the reference while the shaded one indicates the systematic uncertainty. The percentage difference between the results using the SSD and the IM fits (with the systematics) is highlighted in the orange column. Table 2: Lower limits on the half-life $T_{1/2}$ at 90% CI and upper limits on the neutrino-Majoron coupling constant for different Majoron-emitting modes are shown both for SSD and IM fits, including systematics.  Decay mode | $T_{1/2}$ | $g_{ee}^{M}$
---|---|---
limit SSD [yr] | limit IM [yr] | limit SSD | limit IM
$\beta\beta\chi_{0}$ ($n$ = 1) |  $>$ 2.4 $\times$ 1022 |  $>$ 1.6 $\times$ 1022 |  $<$ (4.0–6.9) $\times$ $10^{-5}$ |  $<$ (5.0–8.5) $\times$ $10^{-5}$
$\beta\beta\chi_{0}$ ($n$ = 2) |  $>$ 5.8 $\times$ 1021 |  $>$ 2.7 $\times$ 1021 | - | -
$\beta\beta\chi_{0}$ ($n$ = 3) |  $>$ 2.2 $\times$ 1021 |  $>$ 0.5 $\times$ 1021 |  $<$ 0.053 |  $<$ 0.112
$\beta\beta\chi_{0}\chi_{0}$ ($n$ = 3) |  $>$ 2.2 $\times$ 1021 |  $>$ 0.5 $\times$ 1021 |  $<$ 2.1 |  $<$ 3.1
$\beta\beta\chi_{0}\chi_{0}$ ($n$ = 7) |  $>$ 2.2 $\times$ 1020 |  $>$ 2.0 $\times$ 1020 |  $<$ 2.2 |  $<$ 2.3

###  8.2 Lorentz-violating $2\nu\beta\beta$ decay

The violation of Lorentz and CPT symmetries introduces a perturbation in the
SM-$2\nu\beta\beta$ spectrum. By adding the LV spectrum in the background
model, the measured decay rate of the LV perturbation $\Gamma_{LV}^{m}$ and
the SM $2\nu\beta\beta$ decay $\Gamma_{SM}^{m}$ can be evaluated from the
normalization parameters. Since the countershaded operator
($\dot{a}_{of}^{(3)}$) can assume negative values, under-fluctuations for the
LV component are allowed in the fit. The ratio of the decay rates is directly
proportional to $\dot{a}_{of}^{(3)}$ through

$\dfrac{\Gamma_{LV}^{m}}{\Gamma_{SM}^{m}}=\dot{a}_{of}^{(3)}\cdot
10\cdot\dfrac{\delta G_{LV}}{G_{SM}},$ (11)

where the NMEs cancel out, $\delta G_{LV}$ is the PSF of the Lorentz
perturbation, and $G_{SM}$ is the SM $2\nu\beta\beta$ decay PSF. Finally, the
countershaded operator can be calculated from

$\dot{a}_{of}^{(3)}=C\cdot\dfrac{\Gamma_{LV}^{m}}{\Gamma_{SM}^{m}},$ (12)

where $C$ is a constant value. In this case, the collection of high
$2\nu\beta\beta$ decay statistics plays an important role in constraining
possible Lorentz-violating effects. In both the SSD and the IM fit, the
posterior p.d.f. converges within a range compatible with zero, setting a
double-sided limit at 90% CI on the negative and positive values of the
$\dot{a}_{of}^{(3)}$ parameter. We set the limits by taking into account the
strong anti-correlation of the LV $2\nu\beta\beta$ decay component with the SM
one by calculating the $\dfrac{\Gamma_{LV}^{m}}{\Gamma_{SM}^{m}}$ ratio for
each MCMC sampling. In the SSD fit, the value of the factor $C$ is calculated
from the values of the PSFs in the SSD assumption, $C=(299.0$ $\pm$ $0.2)$
$\times$ $10^{-6}$ GeV-1 Kotila and Gráf (2023), where the error comes mostly
from the uncertainty on the 100Mo $Q_{\beta\beta}$ and it is calculated as in
Ref. Nitescu  _et al._ (2020). In the IM fit, the PSF of SM $2\nu\beta\beta$
decay is explicitly calculated from the $\xi$ parameters.

Figure 3: Double-side limits at 90% CI on the Lorentz-violating countershaded
operator with the SSD approximation for $2\nu\beta\beta$ decay (green) and the
improved description implemented in the fit (orange). The solid color
corresponds to the limit obtained from the reference, while the shaded one
shows the effect of systematics.

The systematic checks listed above are performed on both the SSD and the IM
fits, showing dominant effects for Bremsstrahlung, energy scale, and 90Sr/90Y.
A summary of the results is represented in Fig. 3. The final limits including
systematic in the SSD fit are

$-8.1\times 10^{-6}<\dot{a}_{of}^{(3)}<2.2\times 10^{-6},$

while in the IM fit

$-6.5\times 10^{-6}<\dot{a}_{of}^{(3)}<2.5\times 10^{-5}.$

The strong anti-correlation between the LV spectrum and the $\xi$ parameters
leads to a large broadening of the $\dot{a}_{of}^{(3)}$ posterior, thus a
weaker limit. The most stringent limit on the LV $2\nu\beta\beta$ decay of
100Mo has been set by NEMO-3 Arnold  _et al._ (2019), corresponding to
$(-4.2<\dot{a}_{of}^{(3)}<3.5)\times 10^{-7}$, assuming the single state
dominance for the $2\nu\beta\beta$ decay. In our case, the SSD fit prefers
negative values of the $\dot{a}_{of}^{(3)}$ parameter, producing a stronger
limit on the positive values than the negative ones. Contrarily, in the IM fit
the $\dot{a}_{of}^{(3)}$ converges at positive values.

###  8.3 Sterile neutrino emissions

We investigated a mass range for sterile neutrino masses $m_{N}$ from 0.5 MeV
to 1.5 MeV, with 0.1 MeV steps. The active-sterile mixing strength can be
calculated for different sterile neutrino masses as Bolton  _et al._ (2021);
Agostini  _et al._ (2021)

$\sin^{2}\theta=\frac{G_{SM}}{2G_{\nu N}}\cdot\dfrac{\Gamma_{\nu
N}^{m}}{\Gamma_{SM}^{m}},$ (13)

where $G_{\nu N}$ is the PSF of the sterile neutrino spectrum, calculated for
each sterile neutrino mass considered Kotila and Gráf (2023). The same
considerations on $G_{SM}$ done in the previous section are valid also in this
case. No signal evidence is found for any of the considered masses, setting a
90% CI limit on $\sin^{2}\theta$. We identified as dominant systematic effects
–10% Bremsstrahlung, +1 keV energy scale, and minimal model. The same
systematics are performed in the improved model fit. All the obtained limits
are summarized in Table 3 and Fig. 4. For sterile neutrino masses lower than
1.2 MeV, the $\nu N\beta\beta$ decay component starts correlating with the
$2\nu\beta\beta$ one, and this effect is amplified when the $2\nu\beta\beta$
decay shape is described with the IM. In particular, for low values of $m_{N}$
the shape of $\nu N\beta\beta$ decay spectrum becomes more similar to the SM
component. In the framework of $0\nu\beta\beta$ experiments, GERDA set limits
on $\sin^{2}\theta$, spanning a range for sterile neutrino masses from 0.1 to
0.9 MeV, as reported in Fig. 4. Nevertheless, existing bounds on the active-
sterile mixing strength from $\beta$-decay and solar neutrino experiments
Holzschuh  _et al._ (2000); Schreckenbach  _et al._ (1983); Deutsch  _et al._
(1990); Bellini  _et al._ (2013); Derbin  _et al._ (2018) have already
excluded this region of the parameter space, setting limits on $\sin\theta$ in
the range 10-3–10-2 Holzschuh  _et al._ (2000); Schreckenbach  _et al._
(1983); Deutsch  _et al._ (1990); Bellini  _et al._ (2013); Derbin  _et al._
(2018).

Table 3: Upper limits on the active-sterile mixing strength $\sin^{2}\theta$ for different values of sterile neutrino masses $m_{N}$ Kotila and Gráf (2023). $m_{N}$ [MeV] | PSF [$\times 10^{-18}$] | Limit SSD | Limit IM | diff.
---|---|---|---|---
0.5 | 2.19 |  $<$ 0.033 |  $<$ 0.108 | +229%
0.6 | 1.87 |  $<$ 0.033 |  $<$ 0.085 | +156%
0.7 | 1.57 |  $<$ 0.035 |  $<$ 0.071 | +101%
0.8 | 1.28 |  $<$ 0.039 |  $<$ 0.065 | +64%
0.9 | 1.03 |  $<$ 0.045 |  $<$ 0.061 | +36%
1.0 | 0.91 |  $<$ 0.047 |  $<$ 0.055 | +17%
1.1 | 0.81 |  $<$ 0.049 |  $<$ 0.052 | +5.5%
1.2 | 0.71 |  $<$ 0.051 |  $<$ 0.051 | +0.2%
1.3 | 0.62 |  $<$ 0.053 |  $<$ 0.050 |  $-$4.9%
1.4 | 0.47 |  $<$ 0.063 |  $<$ 0.060 |  $-$5.0%
1.5 | 0.34 |  $<$ 0.074 |  $<$ 0.074 |  $-$0.2%
Figure 4: Limits on the active-sterile mixing strength $\sin^{2}\theta$ as a
function of the sterile neutrino mass $m_{N}$. The blue region represents the
region excluded by the GERDA including systematics, covering the $m_{N}$ range
0.1–0.9 MeV Agostini  _et al._ (2022). The green region is excluded by CUPID-
Mo assuming the improved model to describe the $2\nu\beta\beta$ decay spectral
shape and including systematics. The green dashed line shows the limit
obtained from the reference IM fit. The red lines represent the limits
obtained using the SSD assumption with and without systematics (solid and
dashed lines, respectively). Figure 5: Exotic double-$\beta$ decay spectra
compared to the experimental data with a number of counts corresponding to the
90% CI limit obtained in the SSD fit. The LV $2\nu\beta\beta$ decay spectrum
represents only the limit on the positive side. The black dots correspond to
the experimental data, the gray spectrum is the background reconstructed by
the background model while the blue spectrum represents the SM
$2\nu\beta\beta$ decay.

##  9 Conclusions and outlook

We presented the results of the searches for exotic double-$\beta$ decays with
CUPID-Mo. The analysis exploited the precise spectral shape reconstruction
provided by the background model. No signal evidence has been found for any of
the BSM processes investigated, setting a 90% CI limit on the corresponding
new physics parameter. For the first time, the theoretical uncertainties of
the $2\nu\beta\beta$ decay spectral shape, parameterized through the improved
model description Šimkovic  _et al._ (2018), have been taken into account in
this type of search. This work demonstrated that uncertainties on the
$2\nu\beta\beta$ decay shape induce a significant reduction in sensitivity for
all the processes investigated, requiring better theoretical constraints,
higher statistics, and precise background reconstruction at low energies in
the next-generation experiments. A parallel analysis was performed using a
fixed spectral shape for the $2\nu\beta\beta$ decay (single-state dominance)
to compare the results with those of other experiments. The limits at 90% CI
on the experimental data are represented in Fig. 5. CUPID-Mo set stringent
constraints on the neutrino-Majoron coupling and the Lorenz-violating
countershaded operator, despite the relatively small exposure ($\sim$1.5 kg
$\times$ y), the limits are only a factor 2–10 less stringent than NEMO-3 ones
($\sim$34 kg $\times$ y) Arnold  _et al._ (2015, 2019). The search for $\nu
N\beta\beta$ decay has been performed for the first time using cryogenic
calorimeters. Exploiting the high Q-value of 100Mo, CUPID-Mo data allowed
constraining the active-sterile mixing strength for higher values of $m_{N}$
compared to GERDA. Nevertheless, that region of the parameter space was
already excluded by $\beta$-decay and solar neutrino experiments Holzschuh
_et al._ (2000); Schreckenbach  _et al._ (1983); Deutsch  _et al._ (1990);
Bellini  _et al._ (2013); Derbin  _et al._ (2018). The results of CUPID-Mo
demonstrate the potential of the bolometric technique for exotic
double-$\beta$ decay searches. These promising results motivate the interest
in investigating these BSM processes in the next-generation CUPID experiment.
This study extensively analyzed the main limitations in sensitivity. In
particular, the theoretical uncertainties on the $2\nu\beta\beta$ decay
spectral shape, the uncertainty on the presence of pure $\beta$-emitters in
crystals, and the small statistics are the primary limiting factors. In the
future, improvements in NME calculations can further constrain the $\xi$
parameters of the $2\nu\beta\beta$ decay improved description. For the next-
generation CUPID experiment, the possibility of measuring the 90Sr and 40K
concentration in Li2MoO4 crystals with a sensitivity $\lesssim 10^{-20}$ g/g
will be extremely helpful in overcoming the problem of pure $\beta$-decays.

##  10 Acknowledgments

This work has been performed in the framework of the CUPID-1
(ANR-21-CE31-0014) and LUMINEU programs, funded by the Agence Nationale de la
Recherche (ANR, France). We acknowledge also the support of the P2IO LabEx
(ANR-10-LABX0038) in the framework ”Investissements d’Avenir”
(ANR-11-IDEX-0003-01 – Project ”BSM-nu”) managed by ANR, France. The help of
the technical staff of the Laboratoire Souterrain de Modane and of the other
participant laboratories is gratefully acknowledged. We thank the mechanical
workshops of LAL (now IJCLab) for the detector holders fabrication and
CEA/SPEC for their valuable contribution in the detector conception. F.A.
Danevich, V.V. Kobychev, V.I. Tretyak and M.M. Zarytskyy were supported in
part by the National Research Foundation of Ukraine Grant No. 2020.02/0011.
O.G. Polischuk was supported in part by the project “Investigations of rare
nuclear processes” of the program of the National Academy of Sciences of
Ukraine “Laboratory of young scientists”. J. Kotila is supported by Academy of
Finland (Grant Nos. 3314733, 320062, 345869). F.Š. acknowledges support from
the Slovak Research and Development Agency under Contract No. APVV-22-0413 and
by the Czech Science Foundation (GAČR), project No. 24-10180S. Additionally
the work is supported by the Istituto Nazionale di Fisica Nucleare (INFN) and
by the EU Horizon2020 research and innovation program under the Marie
Sklodowska-Curie Grant Agreement No. 754496. This work is also based on
support by the US Department of Energy (DOE) Office of Science under Contract
Nos. DE-AC02-05CH11231, and by the DOE Office of Science, Office of Nuclear
Physics under Contract Nos. DE-FG02-08ER41551, DE-SC0011091; by the France-
Berkeley Fund, the MISTI-France fund and by the Chateau-briand Fellowship of
the Office for Science & Technology of the Embassy of France in the United
States. This research used resources of the National Energy Research
Scientific Computing Center (NERSC) and the IN2P3 Computing Centre. This work
makes use of the DIANA data analysis software and the background model based
on JAGS, developed by the CUORICINO, CUORE, LUCIFER, and CUPID-0
Collaborations. Russian and Ukrainian scientists have given and give crucial
contributions to CUPID-Mo. For this reason, the CUPID-Mo collaboration is
particularly sensitive to the current situation in Ukraine. The position of
the collaboration leadership on this matter, approved by majority, is
expressed at https://cupid-mo.mit.edu/collaboration#statement. The majority of
the work described here was completed before February 24, 2022.

## References

  * Albert  _et al._ (2014) J. B. Albert _et al._ (EXO-200), Phys. Rev. C 89, 015502 (2014), arXiv:1306.6106 [nucl-ex] .
  * Agostini  _et al._ (2023a) M. Agostini _et al._ (GERDA, (GERDA Collaboration)*), Phys. Rev. Lett. 131, 142501 (2023a), arXiv:2308.09795 [nucl-ex] .
  * Adams  _et al._ (2021) D. Q. Adams _et al._ (CUORE), Phys. Rev. Lett. 126, 171801 (2021), arXiv:2012.11749 [nucl-ex] .
  * Gando  _et al._ (2019) A. Gando _et al._ (KamLAND-Zen), Phys. Rev. Lett. 122, 192501 (2019), arXiv:1901.03871 [hep-ex] .
  * Barabash  _et al._ (2018) A. S. Barabash _et al._ , Phys. Rev. D 98, 092007 (2018), arXiv:1811.06398 [nucl-ex] .
  * Azzolini  _et al._ (2023a) O. Azzolini _et al._ (CUPID), Phys. Rev. Lett. 131, 222501 (2023a), arXiv:2306.14654 [nucl-ex] .
  * Arnold  _et al._ (2019) R. Arnold _et al._ (NEMO-3), Eur. Phys. J. C 79, 440 (2019), arXiv:1903.08084 [nucl-ex] .
  * Augier  _et al._ (2023a) C. Augier _et al._ (CUPID-Mo), Phys. Rev. Lett. 131, 162501 (2023a), arXiv:2307.14086 [nucl-ex] .
  * Gelmini and Roncadelli (1981) G. B. Gelmini and M. Roncadelli, Phys. Lett. B 99, 411 (1981).
  * Chikashige  _et al._ (1980) Y. Chikashige, R. N. Mohapatra, and R. D. Peccei, Phys. Rev. Lett. 45, 1926 (1980).
  * Georgi  _et al._ (1981) H. M. Georgi, S. L. Glashow, and S. Nussinov, Nucl. Phys. B 193, 297 (1981).
  * Schael  _et al._ (2006) S. Schael _et al._ (ALEPH, DELPHI, L3, OPAL, SLD, LEP Electroweak Working Group, SLD Electroweak Group, SLD Heavy Flavour Group), Phys. Rept. 427, 257 (2006), arXiv:hep-ex/0509008 .
  * Bamert  _et al._ (1995) P. Bamert, C. P. Burgess, and R. N. Mohapatra, Nucl. Phys. B 449, 25 (1995), arXiv:hep-ph/9412365 .
  * Burgess and Cline (1994) C. P. Burgess and J. M. Cline, Phys. Rev. D 49, 5925 (1994), arXiv:hep-ph/9307316 .
  * Carone (1993) C. D. Carone, Phys. Lett. B 308, 85 (1993), arXiv:hep-ph/9302290 .
  * Mohapatra  _et al._ (2000) R. N. Mohapatra, A. Perez-Lorenzana, and C. A. de S. Pires, Phys. Lett. B 491, 143 (2000), arXiv:hep-ph/0008158 .
  * Berezinsky and Valle (1993) V. Berezinsky and J. W. F. Valle, Phys. Lett. B 318, 360 (1993), arXiv:hep-ph/9309214 .
  * Heeck (2017) J. Heeck, PoS NuFact2017, 138 (2017).
  * Brune and Päs (2019) T. Brune and H. Päs, Phys. Rev. D 99, 096005 (2019), arXiv:1808.08158 [hep-ph] .
  * Cepedello  _et al._ (2019) R. Cepedello, F. F. Deppisch, L. González, C. Hati, and M. Hirsch, Phys. Rev. Lett. 122, 181801 (2019), arXiv:1811.00031 [hep-ph] .
  * Colladay and Kostelecky (1997) D. Colladay and V. A. Kostelecky, Phys. Rev. D 55, 6760 (1997), arXiv:hep-ph/9703464 .
  * Colladay and Kostelecky (1998) D. Colladay and V. A. Kostelecky, Phys. Rev. D 58, 116002 (1998), arXiv:hep-ph/9809521 .
  * Kostelecky and Mewes (2004a) V. A. Kostelecky and M. Mewes, Phys. Rev. D 70, 031902 (2004a), arXiv:hep-ph/0308300 .
  * Kostelecky and Mewes (2004b) V. A. Kostelecky and M. Mewes, Phys. Rev. D 69, 016005 (2004b), arXiv:hep-ph/0309025 .
  * Kostelecky and Mewes (2012) A. Kostelecky and M. Mewes, Phys. Rev. D 85, 096005 (2012), arXiv:1112.6395 [hep-ph] .
  * Kostelecky and Russell (2011) V. A. Kostelecky and N. Russell, Rev. Mod. Phys. 83, 11 (2011), arXiv:0801.0287 [hep-ph] .
  * Kostelecky and Tasson (2009) V. A. Kostelecky and J. Tasson, Phys. Rev. Lett. 102, 010402 (2009), arXiv:0810.1459 [gr-qc] .
  * Díaz  _et al._ (2013) J. S. Díaz, A. Kostelecký, and R. Lehnert, Phys. Rev. D 88, 071902 (2013), arXiv:1305.4636 [hep-ph] .
  * Chattopadhyay and Patel (2017) P. Chattopadhyay and K. M. Patel, Nucl. Phys. B 921, 487 (2017), arXiv:1703.09541 [hep-ph] .
  * Lee  _et al._ (2013) C.-H. Lee, P. S. Bhupal Dev, and R. N. Mohapatra, Phys. Rev. D 88, 093010 (2013), arXiv:1309.0774 [hep-ph] .
  * Shaposhnikov (2007) M. Shaposhnikov, Nucl. Phys. B 763, 49 (2007), arXiv:hep-ph/0605047 .
  * Deppisch and Pilaftsis (2011) F. F. Deppisch and A. Pilaftsis, Phys. Rev. D 83, 076007 (2011), arXiv:1012.1834 [hep-ph] .
  * Kersten and Smirnov (2007) J. Kersten and A. Y. Smirnov, Phys. Rev. D 76, 073005 (2007), arXiv:0705.3221 [hep-ph] .
  * Pilaftsis (2005) A. Pilaftsis, Phys. Rev. Lett. 95, 081602 (2005), arXiv:hep-ph/0408103 .
  * Bolton  _et al._ (2020) P. D. Bolton, F. F. Deppisch, and P. S. Bhupal Dev, JHEP 03, 170 (2020), arXiv:1912.03058 [hep-ph] .
  * Dasgupta and Kopp (2021) B. Dasgupta and J. Kopp, Phys. Rept. 928, 1 (2021), arXiv:2106.05913 [hep-ph] .
  * Bolton  _et al._ (2021) P. D. Bolton, F. F. Deppisch, L. Gráf, and F. Šimkovic, Phys. Rev. D 103, 055019 (2021), arXiv:2011.13387 [hep-ph] .
  * Agostini  _et al._ (2021) M. Agostini, E. Bossio, A. Ibarra, and X. Marcano, Phys. Lett. B 815, 136127 (2021), arXiv:2012.09281 [hep-ph] .
  * Holzschuh  _et al._ (2000) E. Holzschuh, L. Palermo, H. Stussi, and P. Wenk, Phys. Lett. B 482, 1 (2000).
  * Schreckenbach  _et al._ (1983) K. Schreckenbach, G. Colvin, and F. Von Feilitzsch, Phys. Lett. B 129, 265 (1983).
  * Deutsch  _et al._ (1990) J. Deutsch, M. Lebrun, and R. Prieels, Nucl. Phys. A 518, 149 (1990).
  * Bellini  _et al._ (2013) G. Bellini _et al._ (Borexino), Phys. Rev. D 88, 072010 (2013), arXiv:1311.5347 [hep-ex] .
  * Derbin  _et al._ (2018) A. V. Derbin _et al._ , JETP Letters 108, 499 (2018).
  * Augier  _et al._ (2022) C. Augier _et al._ , Eur. Phys. J. C 82, 1033 (2022), arXiv:2202.08716 [nucl-ex] .
  * Azzolini  _et al._ (2022) O. Azzolini _et al._ (CUPID), Phys. Rev. Lett. 129, 111801 (2022), arXiv:2206.05130 [hep-ex] .
  * Adams  _et al._ (2022) D. Q. Adams _et al._ (CUORE), Nature 604, 53 (2022), arXiv:2104.06906 [nucl-ex] .
  * Augier  _et al._ (2023b) C. Augier _et al._ (CUPID-Mo), Eur. Phys. J. C 83, 675 (2023b), arXiv:2305.01402 [hep-ex] .
  * Šimkovic  _et al._ (2018) F. Šimkovic, R. Dvornický, D. Stefánik, and A. Faessler, Phys. Rev. C 97, 034315 (2018), arXiv:1804.04227 [nucl-th] .
  * Niţescu  _et al._ (2021) O. Niţescu, R. Dvornický, S. Stoica, and F. Šimkovic, Universe 7, 147 (2021).
  * Armstrong  _et al._ (2019) W. R. Armstrong _et al._ (CUPID), (2019), arXiv:1907.09376 [physics.ins-det] .
  * Armengaud  _et al._ (2017) E. Armengaud _et al._ , Eur. Phys. J. C 77, 785 (2017), arXiv:1704.01758 [physics.ins-det] .
  * Poda and Giuliani (2017) D. Poda and A. Giuliani, Int. J. Mod. Phys. A 32, 1743012 (2017), arXiv:1711.01075 [physics.ins-det] .
  * Armengaud  _et al._ (2020) E. Armengaud _et al._ , Eur. Phys. J. C 80, 44 (2020), arXiv:1909.02994 [physics.ins-det] .
  * Armengaud  _et al._ (2013) E. Armengaud _et al._ (EDELWEISS), Astropart. Phys. 47, 1 (2013), arXiv:1305.3628 [physics.ins-det] .
  * Di Domizio  _et al._ (2018) S. Di Domizio, A. Branca, A. Caminata, L. Canonica, S. Copello, A. Giachero, E. Guardincerri, L. Marini, M. Pallavicini, and M. Vignati, JINST 13, P12003 (2018), arXiv:1807.11446 [physics.ins-det] .
  * Gatti and Manfredi (1986) E. Gatti and P. F. Manfredi, Riv. Nuovo Cimento 9, 1 (1986).
  * Alduino  _et al._ (2017) C. Alduino _et al._ (CUORE), Eur. Phys. J. C 77, 13 (2017), arXiv:1609.01666 [nucl-ex] .
  * Agostinelli  _et al._ (2003) S. Agostinelli _et al._ , Nucl. Instr. Meth. A 506, 250 (2003).
  * Ponkratenko  _et al._ (2000) O. A. Ponkratenko, V. I. Tretyak, and Y. G. Zdesenko, Phys. Atom. Nucl. 63, 1282 (2000), arXiv:nucl-ex/0104018 .
  * liv (2023) Geant4 Physics Reference Manual (2023).
  * Plummer  _et al._ (2003) M. Plummer _et al._ , in  _Proceedings of the 3rd international workshop on distributed statistical computing_ , Vol. 124 (Vienna, Austria, 2003) pp. 1–10.
  * Kotila and Iachello (2012) J. Kotila and F. Iachello, Phys. Rev. C 85, 034316 (2012), arXiv:1209.5722 [nucl-th] .
  * Barabash (2020) A. Barabash, Universe 6, 159 (2020), arXiv:2009.14451 [nucl-ex] .
  * Kotila and Gráf (2023) J. Kotila and L. Gráf, Private communication (2023).
  * Agostini  _et al._ (2022) M. Agostini _et al._ (GERDA), JCAP 12, 012 (2022), arXiv:2209.01671 [nucl-ex] .
  * Azzolini  _et al._ (2023b) O. Azzolini _et al._ (CUPID-0), Phys. Rev. D 107, 032006 (2023b), arXiv:2209.09490 [hep-ex] .
  * Arnold  _et al._ (2015) R. Arnold _et al._ (NEMO-3), Phys. Rev. D 92, 072011 (2015), arXiv:1506.05825 [hep-ex] .
  * Kharusi  _et al._ (2021) S. A. Kharusi _et al._ , Phys. Rev. D 104, 112002 (2021), arXiv:2109.01327 [hep-ex] .
  * Gando  _et al._ (2012) A. Gando _et al._ (KamLAND-Zen), Phys. Rev. C 86, 021601 (2012), arXiv:1205.6372 [hep-ex] .
  * Coraggio  _et al._ (2022) L. Coraggio, N. Itaco, G. De Gregorio, A. Gargano, R. Mancino, and F. Nowacki, Phys. Rev. C 105, 034312 (2022), arXiv:2203.01013 [nucl-th] .
  * Pandola  _et al._ (2015) L. Pandola, C. Andenna, and B. Caccia, Nucl. Instrum. Meth. B 350, 41 (2015), arXiv:1410.2002 [physics.med-ph] .
  * Kotila  _et al._ (2015) J. Kotila, J. Barea, and F. Iachello, Phys. Rev. C 91, 064310 (2015), [Erratum: Phys.Rev.C 92, 029903 (2015)], arXiv:1509.05154 [nucl-th] .
  * Agostini  _et al._ (2023b) M. Agostini, G. Benato, J. A. Detwiler, J. Menéndez, and F. Vissani, Rev. Mod. Phys. 95, 025002 (2023b), arXiv:2202.01787 [hep-ex] .
  * Kotila and Iachello (2021) J. Kotila and F. Iachello, Phys. Rev. C 103, 044302 (2021), arXiv:2104.02327 [nucl-th] .
  * Nitescu  _et al._ (2020) O. Nitescu, S. Ghinescu, and S. Stoica, J. Phys. G 47, 055112 (2020), arXiv:2001.04859 [nucl-th] .