CUPID-Mo collaboration

††thanks: Now at: Northwestern University, Evanston, IL 60208, USA

# New measurement of double beta decays of 100Mo to excited states of 100Ru
with the CUPID-Mo experiment

C. Augier  Univ Lyon, Université Lyon 1, CNRS/IN2P3, IP2I-Lyon, F-69622,
Villeurbanne, France  A. S. Barabash  National Research Centre Kurchatov
Institute, Institute of Theoretical and Experimental Physics, 117218 Moscow,
Russia  F. Bellini  Dipartimento di Fisica, Sapienza Università di Roma, P.le
Aldo Moro 2, I-00185, Rome, Italy  INFN, Sezione di Roma, P.le Aldo Moro 2,
I-00185, Rome, Italy  G. Benato  INFN, Laboratori Nazionali del Gran Sasso,
I-67100 Assergi (AQ), Italy  M. Beretta  University of California, Berkeley,
California 94720, USA  L. Bergé  Université Paris-Saclay, CNRS/IN2P3, IJCLab,
91405 Orsay, France  J. Billard  Univ Lyon, Université Lyon 1, CNRS/IN2P3,
IP2I-Lyon, F-69622, Villeurbanne, France  Yu. A. Borovlev  Nikolaev Institute
of Inorganic Chemistry, 630090 Novosibirsk, Russia  L. Cardani  INFN, Sezione
di Roma, P.le Aldo Moro 2, I-00185, Rome, Italy  N. Casali  INFN, Sezione di
Roma, P.le Aldo Moro 2, I-00185, Rome, Italy  A. Cazes  Univ Lyon, Université
Lyon 1, CNRS/IN2P3, IP2I-Lyon, F-69622, Villeurbanne, France  M. Chapellier
Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France  D. Chiesa
Dipartimento di Fisica, Università di Milano-Bicocca, I-20126 Milano, Italy
INFN, Sezione di Milano-Bicocca, I-20126 Milano, Italy  I. Dafinei  INFN,
Sezione di Roma, P.le Aldo Moro 2, I-00185, Rome, Italy  F. A. Danevich
Institute for Nuclear Research of NASU, 03028 Kyiv, Ukraine  M. De Jesus  Univ
Lyon, Université Lyon 1, CNRS/IN2P3, IP2I-Lyon, F-69622, Villeurbanne, France
T. Dixon<EMAIL_ADDRESS>Université Paris-Saclay,
CNRS/IN2P3, IJCLab, 91405 Orsay, France  IRFU, CEA, Université Paris-Saclay,
F-91191 Gif-sur-Yvette, France  L. Dumoulin  Université Paris-Saclay,
CNRS/IN2P3, IJCLab, 91405 Orsay, France  K. Eitel  Karlsruhe Institute of
Technology, Institute for Astroparticle Physics, 76021 Karlsruhe, Germany  F.
Ferri  IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France  B.
K. Fujikawa  Lawrence Berkeley National Laboratory, Berkeley, California
94720, USA  J. Gascon  Univ Lyon, Université Lyon 1, CNRS/IN2P3, IP2I-Lyon,
F-69622, Villeurbanne, France  L. Gironi  Dipartimento di Fisica, Università
di Milano-Bicocca, I-20126 Milano, Italy  INFN, Sezione di Milano-Bicocca,
I-20126 Milano, Italy  A. Giuliani  Université Paris-Saclay, CNRS/IN2P3,
IJCLab, 91405 Orsay, France  V. D. Grigorieva  Nikolaev Institute of Inorganic
Chemistry, 630090 Novosibirsk, Russia  M. Gros  IRFU, CEA, Université Paris-
Saclay, F-91191 Gif-sur-Yvette, France  D. L. Helis  IRFU, CEA, Université
Paris-Saclay, F-91191 Gif-sur-Yvette, France  INFN, Laboratori Nazionali del
Gran Sasso, I-67100 Assergi (AQ), Italy  H. Z. Huang  Key Laboratory of
Nuclear Physics and Ion-beam Application (MOE), Fudan University, Shanghai
200433, PR China  R. Huang  University of California, Berkeley, California
94720, USA  L. Imbert  Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405
Orsay, France  J. Johnston  Massachusetts Institute of Technology, Cambridge,
MA 02139, USA  A. Juillard  Univ Lyon, Université Lyon 1, CNRS/IN2P3,
IP2I-Lyon, F-69622, Villeurbanne, France  H. Khalife  Université Paris-Saclay,
CNRS/IN2P3, IJCLab, 91405 Orsay, France  M. Kleifges  Karlsruhe Institute of
Technology, Institute for Data Processing and Electronics, 76021 Karlsruhe,
Germany  V. V. Kobychev  Institute for Nuclear Research of NASU, 03028 Kyiv,
Ukraine  Yu. G. Kolomensky  University of California, Berkeley, California
94720, USA  Lawrence Berkeley National Laboratory, Berkeley, California 94720,
USA  S.I. Konovalov  National Research Centre Kurchatov Institute, Institute
of Theoretical and Experimental Physics, 117218 Moscow, Russia  J. Kotilla
Department of Physics, University of Jyväskylä, PO Box 35, FI-40014,
Jyväskylä, Finland  Finnish Institute for Educational Research, University of
Jyväskylä, P.O. Box 35, FI-40014 Jyvaäskylä, Finland  Center for Theoretical
Physics, Sloane Physics Laboratory, Yale University, New Haven, Connecticut
06520-8120, USA  P. Loaiza  Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405
Orsay, France  L. Ma  Key Laboratory of Nuclear Physics and Ion-beam
Application (MOE), Fudan University, Shanghai 200433, PR China  E. P. Makarov
Nikolaev Institute of Inorganic Chemistry, 630090 Novosibirsk, Russia  P. de
Marcillac  Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France
R. Mariam  Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France
L. Marini  University of California, Berkeley, California 94720, USA  Lawrence
Berkeley National Laboratory, Berkeley, California 94720, USA  INFN,
Laboratori Nazionali del Gran Sasso, I-67100 Assergi (AQ), Italy  S. Marnieros
Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France  X.-F. Navick
IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France  C. Nones
IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France  E. B.
Norman  University of California, Berkeley, California 94720, USA  E. Olivieri
Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France  J. L.
Ouellet  Massachusetts Institute of Technology, Cambridge, MA 02139, USA  L.
Pagnanini  INFN, Gran Sasso Science Institute, I-67100 L’Aquila, Italy  INFN,
Laboratori Nazionali del Gran Sasso, I-67100 Assergi (AQ), Italy  L. Pattavina
INFN, Laboratori Nazionali del Gran Sasso, I-67100 Assergi (AQ), Italy  Physik
Department, Technische Universität München, Garching D-85748, Germany  B. Paul
IRFU, CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France  M. Pavan
Dipartimento di Fisica, Università di Milano-Bicocca, I-20126 Milano, Italy
INFN, Sezione di Milano-Bicocca, I-20126 Milano, Italy  H. Peng  Department of
Modern Physics, University of Science and Technology of China, Hefei 230027,
PR China  G. Pessina  INFN, Sezione di Milano-Bicocca, I-20126 Milano, Italy
S. Pirro  INFN, Laboratori Nazionali del Gran Sasso, I-67100 Assergi (AQ),
Italy  D. V. Poda  Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay,
France  O. G. Polischuk  Institute for Nuclear Research of NASU, 03028 Kyiv,
Ukraine  S. Pozzi  INFN, Sezione di Milano-Bicocca, I-20126 Milano, Italy  E.
Previtali  Dipartimento di Fisica, Università di Milano-Bicocca, I-20126
Milano, Italy  INFN, Sezione di Milano-Bicocca, I-20126 Milano, Italy  Th.
Redon  Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France  A.
Rojas  LSM, Laboratoire Souterrain de Modane, 73500 Modane, France  S. Rozov
Laboratory of Nuclear Problems, JINR, 141980 Dubna, Moscow region, Russia  V.
Sanglard  Univ Lyon, Université Lyon 1, CNRS/IN2P3, IP2I-Lyon, F-69622,
Villeurbanne, France  J. A. Scarpaci  Université Paris-Saclay, CNRS/IN2P3,
IJCLab, 91405 Orsay, France  B. Schmidt  Lawrence Berkeley National
Laboratory, Berkeley, California 94720, USA  Y. Shen  Key Laboratory of
Nuclear Physics and Ion-beam Application (MOE), Fudan University, Shanghai
200433, PR China  V. N. Shlegel  Nikolaev Institute of Inorganic Chemistry,
630090 Novosibirsk, Russia  V. Singh  University of California, Berkeley,
California 94720, USA  C. Tomei  INFN, Sezione di Roma, P.le Aldo Moro 2,
I-00185, Rome, Italy  V. I. Tretyak  Institute for Nuclear Research of NASU,
03028 Kyiv, Ukraine  V. I. Umatov  National Research Centre Kurchatov
Institute, Institute of Theoretical and Experimental Physics, 117218 Moscow,
Russia  L. Vagneron  Univ Lyon, Université Lyon 1, CNRS/IN2P3, IP2I-Lyon,
F-69622, Villeurbanne, France  M. Velázquez  Université Grenoble Alpes, CNRS,
Grenoble INP, SIMAP, 38420 Saint Martin d’HẀeieres, France  B. Welliver
University of California, Berkeley, California 94720, USA  L. Winslow
Massachusetts Institute of Technology, Cambridge, MA 02139, USA  M. Xue
Department of Modern Physics, University of Science and Technology of China,
Hefei 230027, PR China  E. Yakushev  Laboratory of Nuclear Problems, JINR,
141980 Dubna, Moscow region, Russia  M. Zarytskyy  Institute for Nuclear
Research of NASU, 03028 Kyiv, Ukraine  A. S. Zolotarova  Université Paris-
Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France

###### Abstract

The CUPID-Mo experiment, located at Laboratoire Souterrain de Modane (France),
was a demonstrator experiment for CUPID. It consisted of an array of 20
Li${}_{2}^{100}$MoO4 (LMO) calorimeters each equipped with a Ge light detector
(LD) for particle identification. In this work, we present the result of a
search for two-neutrino and neutrinoless double beta decays of 100Mo to the
first 0+ and $2^{+}$ excited states of 100Ru using the full CUPID-Mo exposure
(2.71 kg$\times$yr of LMO). We measure the half-life of $2\nu\beta\beta$ decay
to the $0^{+}_{1}$ state as $T_{1/2}^{2\nu\rightarrow 0_{1}^{+}}=(7.5\pm 0.8\
\text{(stat.)}\ ^{+0.4}_{-0.3}\ \text{(syst.)})\times 10^{20}\ \mathrm{yr}$.
The bolometric technique enables measurement of the electron energies as well
as the gamma rays from nuclear de-excitation and this allows us to set new
limits on the two-neutrino decay to the $2_{1}^{+}$ state of
$T^{2\nu\rightarrow 2_{1}^{+}}_{1/2}>4.4\times 10^{21}\ \mathrm{yr}\
\text{(90\% c.i.)}$ and on the neutrinoless modes of $T_{1/2}^{0\nu\rightarrow
2_{1}^{+}}>2.1\times 10^{23}\ \mathrm{yr}\ \text{(90\% c.i.)}$,
$T_{1/2}^{0\nu\rightarrow 0_{1}^{+}}>1.2\times 10^{23}\ \mathrm{yr}\
\text{(90\% c.i.)}$. Information on the electrons spectral shape is obtained
which allows us to make the first comparison of the single state (SSD) and
higher state (HSD) $2\nu\beta\beta$ decay models for the $0_{1}^{+}$ excited
state of 100Ru.

††preprint: APS/123-QED

##  I Introduction

Two-neutrino double beta decay ($2\nu\beta\beta)$ is an allowed Standard Model
process which occurs in some even-even nuclei for which single beta decays are
energetically forbidden or heavily disfavored due to large changes in angular
momentum [1, 2]. These decays have been observed to the $0^{+}$ ground states
($0^{+}_{\text{g.s.}}$) in eleven nuclei and to the first zero-plus excited
state ($0_{1}^{+}$) for two nuclei, 100Mo and 150Nd [1, 2, 3]. In addition, if
neutrinos are Majorana particles then an additional decay mode becomes
allowed: neutrinoless double beta decay ($0\nu\beta\beta$) [4, 5, 6]. Whereas
2$\nu\beta\beta$ decay conserves lepton number, 0$\nu\beta\beta$ decay would
violate this symmetry [7, 4, 8] and provide a possible path to explain the
predominance of matter over antimatter in the universe [8, 4, 9, 10, 11].
One of the most promising experimental techniques to search for
$0\nu\beta\beta$ decay are cryogenic calorimeters. These detectors provide an
excellent energy resolution, high detection efficiency and are scalable to
tonne scale arrays, such as the CUORE experiment (Cryogenic Underground
Observatory for Rare Events) [12]. The CUPID experiment [13] (CUORE Upgrade
with Particle IDentification) is an upgrade for CUORE which will use a
scintillating bolometer technology to discriminate $\beta/\gamma$ and $\alpha$
particles. CUPID-Mo [14] was a demonstrator experiment for CUPID aiming to
validate the technology of lithium molybdate (LMO) scintillating calorimeters.
Multiple models exist to describe 0$\nu\beta\beta$ decay [7, 8, 15, 16, 17,
18, 19, 20, 21], however the minimal extension to the Standard Model needed to
explain 0$\nu\beta\beta$ decay is the light Majorana neutrino exchange
mechanism [22], in which 0$\nu\beta\beta$ decay is mediated by a Majorana
neutrino. In this framework, the decay rate of 0$\nu\beta\beta$ decay is
related to the effective Majorana neutrino mass $\langle
m_{\beta\beta}\rangle$ (a weighted sum of the three neutrino masses) by:

$(T_{1/2}^{0\nu})^{-1}=G_{0\nu}g_{A}^{4}|M_{0\nu}|^{2}\frac{\langle
m_{\beta\beta}\rangle^{2}}{m_{e}^{2}},$ (1)

where $T_{1/2}^{0\nu}$ is the 0$\nu\beta\beta$ decay half-life, $G_{0\nu}$ is
the phase space factor, $M_{0\nu}$ is the nuclear matrix element (NME),
$m_{e}$ is the electron mass and $g_{A}$ is the axial vector coupling
constant. The phase space factor $G_{0\nu}$ can be computed accurately, but
the NME is the result of complicated many-body nuclear physics calculations
[23]. Several models with different approximations are used, but these only
agree to within a factor of a few on the value of the NME (cf. [24, 25, 26,
27]). At present the upper limit on the effective Majorana mass, $\langle
m_{\beta\beta}\rangle$, ranges from 60 to 600 meV [8, 28, 29, 12, 30]. Despite
its low exposure, CUPID-Mo has set the most stringent limit on
$0\nu\beta\beta$ decay in 100Mo, with $T_{1/2}^{0\nu}>1.8\times 10^{24}\
\mathrm{yr}$ or $\langle m_{\beta\beta}\rangle<(0.28$ – $0.49)\,\mathrm{eV}$
[31].

For $2\nu\beta\beta$ decay, the half-life is related to the NME as:

$(T_{1/2}^{2\nu})^{-1}=G_{2\nu}|M^{\text{eff}}_{2\nu}|^{2},$ (2)

where $G_{2\nu}$ is the phase space factor and and $M^{\text{eff}}_{2\nu}$ is
the dimensionless effective nuclear matrix element including the factor
$g_{A}^{2}$, $M^{\text{eff}}_{2}\nu=g_{A}^{2}\times M_{2\nu}$ [2].
Measurements of $T_{1/2}^{2\nu}$ for both ground states (g.s.) and excited
states (e.s.) of a daughter nucleus provide experimental data which can be
used to validate the methods to calculate the NME $M_{2\nu}$, and by
extension, provide better confidence in the calculation of $M_{0\nu}$.
Whilst the phase space, and therefore decay probability, is lower for decays
to e.s. compared to transitions to the g.s., the presence of monochromatic
$\gamma$’s can lead to a very clear detection signature. In the case of a
positive observation of $0\nu\beta\beta$ decay, measurements of the
$0\nu\beta\beta$ decay to e.s. will be a useful tool to understand the decay
mechanism [3]. Additionally, whilst the $0\nu\beta\beta$ decay to $0_{1}^{+}$
state is expected to be suppressed relative to the decay to
$0^{+}_{\text{g.s.}}$ state, it can still improve the overall sensitivity to
$\langle m_{\beta\beta}\rangle$. Finally, a $2\nu\beta\beta$ decay to
$2_{1}^{+}$ e.s. is expected to be strongly suppressed by angular momentum
conservation. However, this mode is expected to be much more likely in a
framework with Bosonic neutrinos [32, 33] and can be used to test this model.

Figure 1: Transition scheme for the decay of 100Mo to 100Ru showing the double
beta decays to the $0_{1}^{+}$ (1130.3 keV) and $2_{1}^{+}$ (539.5 keV)
states. The decay to the $0_{1}^{+}$ state proceeds via the $2_{1}^{+}$ state
with $\sim 100\%$ probability, so both a 539.5 keV and 590.8 keV $\gamma$ are
emitted in every decay, while only the $539.5$ keV $\gamma$ is emitted in the
double beta decay to the $2_{1}^{+}$ state.

We can exploit the specific combination of energy depositions associated with
these events to measure $\beta\beta$ decays to e.s. As shown in Figure 1, the
two $\beta$’s are accompanied by one or more $\gamma$ particles. These have a
much longer range ($\sim$ $10\ \mathrm{cm}$ [34] for 1 MeV $\gamma$) compared
to $\beta^{-}$ ($\sim$ 2 mm [35]) in LMO. A $\beta\beta$ e.s. event therefore
has a high probability to induce a multi-detector signature with the
$\gamma$’s escaping the detector where the decay took place and being measured
in neighboring detectors. The additional information provided by these multi-
detector events allows us to significantly reduce the background rate while
the $\gamma$ peaks provide a clean signature to robustly measure the decay
rate [36].
Since decays to higher e.s. are generally disfavored by the lower phase space
due to lower Q-value (and also by angular momentum conservation for $2^{+}$
decays) we focus on the lowest energy e.s.: the first zero plus ($0_{1}^{+}$)
and two plus ($2_{1}^{+}$) e.s. of 100Ru at 1130.3 keV and 539.5 keV (see
Figure 1).
The $2\nu\beta\beta$ decay of 100Mo to $0_{1}^{+}$ e.s. has been measured in
several experiments [37, 38, 39, 40, 41, 42], mostly using high purity
Germanium (HPGe) counters (cf. [3, 43, 2] for a recent review) with an
external Mo (enriched) source and also NEMO-3 using a tracking detector and
external sources. For the HPGe measurements given that the source and detector
did not coincide, only $\gamma$’s were detected. The $2\nu\beta\beta$ decay to
$2_{1}^{+}$ e.s. has not been observed and only lower limits on the half-life
are available [42, 39]. Placing constraints on this decay is complicated by
the overlap with the decays to the $0_{1}^{+}$ state since both decays include
a 539.5 keV $\gamma$ emission. By using a setup where the source is embedded
in the detector, such as cryogenic calorimeters, we are able to measure the
electron energies as well as the $\gamma$’s. This provides additional
information to separate the $2\nu\beta\beta$ decays to the $2_{1}^{+}$ and
$0_{1}^{+}$ e.s., and to validate that the decay is indeed
$2\nu\beta\beta\rightarrow 0_{1}^{+}$ and not a background process. This type
of detector also allows us to effectively distinguish the neutrinoless and
two-neutrino decay modes, which is not possible with HPGe detectors using an
external $\beta\beta$ source.

Figure 2: CUPID-Mo experimental geometry showing (Left) the CUPID-Mo towers
mounted in the EDELWEISS cryostat and single detector module (Right) with LMO
crystal, PTFE clamps, NTD thermistor, Si heater and Cu holder. The Ge LD is
seen through the LMO crystal.

In this work, we present an analysis of
$2\nu\beta\beta/0\nu\beta\beta\rightarrow 0_{1}^{+}/2_{1}^{+}$ e.s. which
exploits the measurement of both $\beta$ and $\gamma$ energies. In Section II
we describe the CUPID-Mo experiment; we introduce the search technique and
Geant4 Monte Carlo simulations in Section III. In Section IV, an overview of
the experimental data and treatment of multi-site events used for this
analysis are given. We describe our Bayesian analysis framework in Section V,
our treatment of systematic uncertainties in Section VI, and the obtained
results in Section VII. We conclude in Section VIII by comparing our
measurement of $2\nu\beta\beta\rightarrow 0_{1}^{+}$ state to theoretical
calculations and we discuss the prospects to further investigate the decay
mechanism as well as sensitivity to models beyond the light Majorana neutrino
exchange.

##  II Experiment

CUPID-Mo consisted of an array of twenty $\sim$ 200 g LMO detectors (Figure 2)
enriched to $\sim$ 97 % in 100Mo, with a total mass of 4.16 kg of LMO [14] and
2.26 kg of 100Mo. These detectors were operated as cryogenic calorimeters at
$\sim$ 20 mK, in the EDELWEISS [44] cryostat at the Laboratoire Souterrain de
Modane (LSM), France. This technique allows for a very high detection
efficiency ($\sim$ 75 % for $0\nu\beta\beta$ decay to g.s.) due to the
detector containing the source, and excellent energy resolution $(\sigma\sim
0.1\%)$ at 3 MeV [14]. In addition, CUPID-Mo employed a dual readout with
twenty cryogenic light detectors (LDs) consisting of Ge wafers also operated
as calorimeters. These allow for discrimination between $\alpha$’s and
$\beta/\gamma$’s since LMO is a scintillating crystal and the amount of light
produced for $\beta/\gamma$ is about five times higher compared to that
produced by $\alpha$ particles of the same energy [14, 45].

Each LMO crystal is assembled into an independent detector module with a Ge
LD, NTD-Ge [46] thermistors on both LMO and LD, Si heater, copper holder,
reflective foil (3M Vikuiti™, which guides the light to the LDs), and PTFE
clamps (shown in Figure 2). While this design allows for a modular pre-
assembly and is compatible with the operation of multiple payloads in the
EDELWEISS cryogenic system, it results in a relatively high copper/LMO mass
ratio of $\sim 1.5$. The individual detector modules are arranged in an array
of five towers (see Figure 2) so that each LMO detector faces two LDs (apart
from those on the top floor which only have a lower light detector).

CUPID-Mo took data between 2019 and 2020, collecting a total exposure of 2.71
kg$\times$yr of LMO. The data-taking was organized into 9 datasets, periods of
stable operations lasting around a month. As in [47], three short periods of
data-taking which could not be calibrated to the same accuracy were discarded.

##  III Search Technique and Simulations

In $\beta\beta$ decays to e.s. the electrons are accompanied by $\gamma$’s so
they often reconstruct in multi-site or coincident events where multiple
detectors are triggered simultaneously. The particular set of energies
obtained allows us to make a pure selection of e.s. decay events with very low
background. We define the multiplicity, $\mathcal{M}$, as the number of
crystals subject to simultaneous energy depositions above the analysis energy
threshold, chosen as 40 keV, well above the detector trigger thresholds and
within a time window of $\pm$ 10 ms well above the detector time resolution
[31].

For an event which is a time coincidence between several LMO detectors, our
minimal experimental observable information is the energy deposited in each
detector. We define energies $E_{1},E_{2},E_{3},..,E_{\mathcal{M}}$ as the
energy deposited in each detector, sorted so that:

$E_{1}>E_{2}>E_{3}...>E_{\mathcal{M}}.$ (3)

In principle, the most sensitive approach exploiting all the information
contained in the data, would be a fit directly to this multidimensional energy
spectrum. However, this analysis has significant complexity due to the
difficulty of quantifying the background shape in this high dimensional space.

An alternative approach that has been used by the CUORE-0, CUORE, and CUPID-0
experiments [48, 49, 50] is to select categories of signal signatures
consistent with being from an e.s. decay. For each category only the energy
variable containing the $\gamma$ (or two electrons - $\beta\beta$ for
$0\nu\beta\beta$ decay) peak is fitted (the peak energy $E_{p}$). The other
energy variables are referred to as “projected out” energies. This technique
has the advantage that the decay rate can be extracted as a set of one
dimensional fits to a peak over a background, which is a simple and robust
technique. We use this approach in our analysis, however we modify this
technique by dividing up categories based on their projected out energy
variables. Since the decays to different e.s. significantly overlap this
allows us to better discriminate between them. Furthermore, it avoids
integrating together regions with different background levels.

###  III.1 MC simulations

We use dedicated Monte Carlo (MC) simulations to identify the most promising
categories of events to include and optimize the associated energy cuts. The
objective is to identify categories which possess a high signal to background
ratio, to define the energy variable to project onto, to select the energy
boundaries for each signal category, and to assess the containment efficiency
therein. A Geant4 [51] model has been created which implements the CUPID-Mo
detector into the existing EDELWEISS MC software package [52]. Special care is
taken to implement the detector structure as accurately as possible, in
particular accounting for the individual dimensions of each LMO crystal, and
the passive holder material close to the detectors. We simulate both
$2\nu\beta\beta$ decays and $0\nu\beta\beta$ decays to the $0_{1}^{+}$ and
$2_{1}^{+}$ e.s.. This accounts for the angular correlation between 540 keV
and 591 keV $\gamma$’s from the $0_{1}^{+}$
($0^{+}_{1}-2^{+}_{1}-0^{+}_{\text{g.s}}$) cascade as [53]:

$\rho(\theta)=1-3\cos^{2}{(\theta)}+4\cos^{4}{(\theta)}.$ (4)

Here $\rho(\theta)$ is the probability for $\gamma$’s to be emitted at angle
$\theta$ from each other. Particles are then propagated through the
experimental geometry using the Livermore low energy physics list [54],
applicable down to very low energy (250 eV). For the $2\nu\beta\beta$ decay to
$0_{1}^{+}$ e.s. simulations are performed using both the single state
dominance (SSD) and higher state dominance (HSD) models [55]. Currently, no
experimental data favors one model over another for the $0_{1}^{+}$ state,
however the SSD model is strongly favored for the ground state decay [56, 57]
and so we use it as our default. The difference to the HSD model is then
treated as a systematic uncertainty. We combine the MC simulations with inputs
from experimental data to reproduce the detector response of each individual
detector, accounting for:

  * •

The energy threshold (set at 40 keV) [31];

  * •

The multiplicity of events;

  * •

Energy resolution of each detector-dataset pair (energy dependent);

  * •

Rejected periods of detector operations due to periods of cryogenic
instability (i.e. due to environmental disturbances);

  * •

Scintillation light, light detector energy resolution, and cuts based on the
LDs (see Section IV);

Figure 3: Distribution of multiplicities (MC simulation) for each e.s. decay.
Due to the small crystal array and significant amount of non-sensitive holder
material, both decay modes to the ($0_{1}^{+},2_{1}^{+}$) state are dominated
by single crystals hits, (10 – 20%) of events cause energy deposits in two
crystals, and $\mathcal{O}$(1%) in three or more.

Figure 4: Monte Carlo simulated two-dimensional $\mathcal{M}_{2}$ energy
distributions for decays to e.s.Top left: $2\nu\beta\beta$ decay to
$0_{1}^{+}$. Top Right: $2\nu\beta\beta$ decay to $2_{1}^{+}$. Bottom Left:
$0\nu\beta\beta$ decay to $0_{1}^{+}$. Bottom Right: $0\nu\beta\beta$ decay to
$2_{1}^{+}$. We show red boxes to highlight the $\mathcal{M}_{2}$ signal
categories used in the analysis (see Tables 1, 2) and we label the associated
category numbers.

In our MC simulations, we observe that most excited states events are in
$\mathcal{M}_{1-2}$ with a small contribution of $\mathcal{M}_{3}$ (Figure 3).
The large fraction of $\mathcal{M}_{1}$ events is caused by the significant
volume of passive material (Cu housing) and a relatively small modular
detector array. We exclude $\mathcal{M}_{1}$ 2$\nu\beta\beta$ decay events as
the continuous energy spectrum of the $\beta$’s does not produce a peak to
perform a fit. For $0\nu\beta\beta$ decays, there are mono-energetic peaks in
$\mathcal{M}_{1}$, albeit with a high background. We therefore focus on
$\mathcal{M}_{2}$ events which have the highest sensitivity, and we also
include some categories of $\mathcal{M}_{3}$ events.

###  III.2 $2\nu\beta\beta$ decay categories

We first describe categories for $2\nu\beta\beta$ decay analysis. For the
analysis of $\mathcal{M}_{2}$ events the energies can be visualized as two-
dimensional histograms. We plot these two-dimensional energy spectra for both
e.s. decays ($2\nu\beta\beta$ and $0\nu\beta\beta$ decay) to $0_{1}^{+}$ and
$2_{1}^{+}$ e.s. in Figure 4. When selecting categories of events to use in
our analysis, we require at least one peak to perform a fit. This selection is
equivalent to looking for lines or points on the two-dimensional plane.

For $2\nu\beta\beta$ decays in MC simulation we observe a very clear signature
of events where either $E_{1}$ or $E_{2}$ has a $\gamma$ energy (either
$539.5$ or $590.8$ keV) as shown in Figure 4 (top). These events can be
interpreted as the $\beta\beta$ being contained in one detector, giving rise
to a continuous distribution, and a $\gamma$ that is fully contained in a
second LMO detector. In the case of the decay to the $0_{1}^{+}$ state, the
second $\gamma$ either escapes the detector or is contained in the crystal
with the $\beta\beta$’s. We notice that the $2_{1}^{+}$ and $0_{1}^{+}$ decays
(partially) overlap and the vertical/horizontal lines cover a wide range of
projected out ($\beta\beta$) energies. This motivates us to divide the
observed signal distribution into five slices, two for the horizontal and
three for the vertical lines since the vertical lines cover a larger range in
energy. The boundaries of this categories are defined to maximise the
experimental sensitivity using MC simulations as explained in Appendix A. This
leads to the energy categories highlighted in Figure 4 (top) and listed in
Table 1 (index 0 to 4).
To identify $\mathcal{M}_{3}$ categories we note that in MC simulations
generally the electrons carry the largest kinetic energy and are thus likely
to be contained in the $E_{1}$ detector due to the relatively low $\gamma$
energies. We hence focus our search on events with $\gamma$’s in the $E_{2}$
and $E_{3}$ variables. In particular, this leads to two categories, one where
the $\gamma$ energy is split between $E_{2},E_{3}$ and another where one
$\gamma$ is fully contained in $E_{2}$ with a Compton scatter event in
$E_{3}$.
Combined with the $\mathcal{M}_{2}$ categories this leads to a total of 7
independent categories used for the $2\nu\beta\beta$ decay analysis.

Table 1: Categories for $2\nu\beta\beta$ decay to $0_{1}^{+}$, $2_{1}^{+}$ e.s., these are given a name describing the multiplicity, peak energy variable and then an index (if necessary to distinguish between categories). The table lists the energy cuts for the variable that is projected out, while the projected energy range for each category is 500 – 650 keV. For every category the $\gamma$’s are either 540 or 591 keV. Cat. | Multiplicity | Peak Variable | Energy cuts | Interpretation
---|---|---|---|---
0 | 2 | $E_{1}$ | $E_{2}\in[220,410]$ |  $\gamma$ in $E_{1}$, $\beta\beta$ in $E_{2}$
1 | 2 | $E_{1}$ | $E_{2}\in[410,500]$ |  $\gamma$ in $E_{1}$, $\beta\beta$ in $E_{2}$
2 | 2 | $E_{2}$ | $E_{1}\in[500,890]$ |  $\gamma$ in $E_{2}$, $\beta\beta$ in $E_{1}$
3 | 2 | $E_{2}$ | $E_{1}\in[890,1190]$ |  $\gamma$ in $E_{2}$, $\beta\beta$ in $E_{1}$
4 | 2 | $E_{2}$ | $E_{1}\in[1190,1880]$ |  $\gamma$ in $E_{2}$, $\beta\beta$ in $E_{1}$
5 | 3 | $E_{2}$ | - |  $\gamma$ in $E_{2}$, $\beta\beta$ in $E_{1}$, Compt. in $E_{3}$
6 | 3 | $E_{2}+E_{3}$ | $E_{2}<500$ |  $\gamma$ in $E_{2}+E_{3}$ (shared Compt.), $\beta\beta$ in $E_{1}$

###  III.3 $0\nu\beta\beta$ decay categories

Unlike for $2\nu\beta\beta$ decay, in $0\nu\beta\beta$ decay the electrons
have a fixed total energy so we are also able to search for a peak from the
summed electron energy. Therefore a different set of categories are necessary,
however the same strategy can be used. For $\mathcal{M}_{1}$ events we include
three categories corresponding to the $\beta\beta$ and possibly $\gamma$’s
being contained in the detector where the decay took place. For
$\mathcal{M}_{2}$ events, we show the two-dimensional distributions in Figure
4 (bottom two figures). We observe a variety of features, vertical, horizontal
diagonal lines and points. Each of these require a different set of energy
cuts to contain the signal, as shown in Table 2. These cuts are carefully
constructed to ensure that none of the categories overlap. Some categories
feature a mono-energetic peak in both energy variables, for these we still
choose one as the peak variable and use a range of $\pm 10$ keV for the other,
projected energy, much larger than the energy resolution. We employ a similar
strategy to the $2\nu\beta\beta$ decay categories, dividing the horizontal
lines, events with a $\beta\beta$ peak ($1904,2444,2494$ keV) in $E_{1}$ and a
Compton scatter in $E_{2}$, into 3 categories each to account for the large
changes in background level. This is not necessary for other categories due to
the low background, or a small range in the projected out energy. We neglect
some possible categories which have very small containment efficiency (events
with $E_{2}>600$ keV), that overlap significantly with $2\nu\beta\beta$ decay
(events with $E_{2}\sim 500$ – $650$ keV, $E_{1}<1904$ keV) or with 214Bi
background, diagonal lines with energy summing to $\sim$ 2400 keV. This leads
to 10 $\mathcal{M}_{2}$ $0\nu\beta\beta$ decay categories which are most
easily interpreted geometrically (see the boxes in Figure 4).

As in $2\nu\beta\beta$ decay analysis we also consider $\mathcal{M}_{3}$
events focusing on those where the electrons are contained in $E_{1}$. We
choose four categories: For events with $E_{1}\sim 1904$ keV, we split the
data into $E_{2}\in[500,650]$ keV, $E_{2}+E_{3}\in[500,650]$ keV or neither of
these two. For events with $E_{1}\sim 2400$ keV the background is expected to
be low so we just include one category. This procedure leads to 17
experimental signatures as shown in Table 2.

Table 2: Categories for $0\nu\beta\beta$ to excited states, for each category we give an index, a name, the projected energy variable and its energy window, the cuts used to define the category and a physical interpretation. Most of these cuts are more easily visualized in Figure 4.  Cat. | Multiplicity | Peak Variable | Range [keV] | Energy cuts [keV] | Interpretation
---|---|---|---|---|---
0 | 1 | $E_{1}$ |  $1864$–$1944$ |  |  $\beta\beta$ in $E_{1}$, $\gamma$ escape
1 | 1 | $E_{1}$ |  $2384$–$2534$ |  |  $\beta\beta$ in $E_{1}$, $\gamma$ escape
2 | 1 | $E_{1}$ |  $2984$–$3084$ |  |  $\beta\beta$ and $\gamma$ in $E_{1}$
3 | 2 | $E_{1}$ |  $1864$–$1944$ | $E_{2}<180$ |  $\beta\beta$ in $E_{1}$, Compt. $\gamma$ in $E_{2}$
4 | 2 | $E_{1}$ |  $1864$–$1944$ | $E_{2}\in[180,310]$ |  $\beta\beta$ in $E_{1}$, Compt. $\gamma$ in $E_{2}$
5 | 2 | $E_{1}$ |  $1864$–$1944$ | $E_{2}\in[310,500]$ |  $\beta\beta$ in $E_{1}$, Compt. $\gamma$ in $E_{2}$
6 | 2 | $E_{1}$ |  $2384$–$2534$ | $E_{2}<180$ |  $\beta\beta$ in $E_{1}$, Compt. $\gamma$ in $E_{2}$
7 | 2 | $E_{1}$ |  $2384$–$2534$ | $E_{2}\in[180,280]$ |  $\beta\beta$ in $E_{1}$, Compt. $\gamma$ in $E_{2}$
8 | 2 | $E_{1}$ |  $2384$–$2534$ | $E_{2}\in[280,500]$ |  $\beta\beta$ in $E_{1}$, Compt. $\gamma$ in $E_{2}$
9 | 2 | $E_{2}$ |  $500$–$650$ | $E_{1}\in[1894,2904]$ |  $\beta\beta$ in $E_{1}$, $\gamma$ in $E_{2}$
10 | 2 | $E_{2}$ |  $500$–$650$ |  $E_{1}\in[1944,2384]$, $E_{1}+E_{2}>2534$ |  $\beta\beta$ \+ Compt. $E_{1}$, $\gamma$ in $E_{2}$
11 | 2 |  $E_{1}$+$E_{2}$ |  $2984$–$3084$ |  $E_{1}>2534$, $E_{2}<500$ |  $\beta\beta$ \+ Compt. $E_{1}$, Compt. $E_{2}$
12 | 2 | $E_{1}-E_{2}$ |  $1734$–$2034$ |  $E_{2}<650$, $E_{1}>2384$, $E_{1}$+$E_{2}\in[3024,3044]$ |  $\beta\beta$ in $E_{1}$, $\gamma$ in $E_{2}$
13 | 3 | $E_{2}+E_{3}$ |  $500$–$650$ |  $E_{1}\in[1894,1914]$, $E_{2}<500$ |  $\beta\beta$ $E_{1}$, Compt. $E_{2}/E_{3}$
14 | 3 | $E_{2}$ |  $500$–$650$ | $E_{1}\in[1894,1914]$ |  $\beta\beta$ in $E_{1}$, $\gamma$ in $E_{2}$, Compt. $E_{3}$
15 | 3 | $E_{1}$ |  $1854$–$1954$ |  $E_{2}\not\in[500,650]$, $E_{2}+E_{3}\not\in[500,650]$ |  $\beta\beta$ in $E_{1}$ and Compt. $E_{2}/E_{3}$
16 | 3 | $E_{1}$ |  $2384$–$2534$ |  |  $\beta\beta$ in $E_{1}$

###  III.4 Determination of signal shape and efficiency

We use MC simulations to determine the signal shapes and efficiencies for our
Bayesian analysis. For each category we make the energy cuts as described in
Tables 1, 2. We also apply a set of selection cuts to the data and MC
simulations to remove events that likely arise from a known $\gamma$-ray
background peak. The cuts are chosen to correspond to roughly $\pm\Delta
E(\text{FWHM})$ with the values shown in Table 3.

Table 3: Table showing the $\gamma$ exclusion cuts used to remove events likely originating from a background $\gamma$ ray. Variable | Cut [keV] | Decay
---|---|---
$E_{1}$ | $1333\pm 5$ |  60Co
$E_{1}$ | $1173\pm 5$ |  60Co
$E_{1}+E_{2}$ | $1333\pm 5$ |  60Co
$E_{1}+E_{2}$ | $1173\pm 5$ |  60Co
$E_{1}+E_{2}$ | $1461\pm 5$ |  40K
$E_{1}$ | $2615\pm 7$ |  208Tl
$E_{1}+E_{2}$ | $2615\pm 7$ |  208Tl

After the projection onto the peak energy variable we extract one-dimensional
histograms for each decay mode ($0_{1}^{+}$ or $2_{1}^{+}$) in each category
from MC simulations.
From these histograms we obtain the containment efficiency
$\varepsilon_{\text{cont},r}^{0+/2+}$, or the probability for a decay to
reconstruct in category $r$. We also compute the total containment efficiency
(summed over all categories) as $4.4/2.0\ \%$ for $2\nu\beta\beta$ to
$0_{1}^{+}/2_{1}^{+}$ e.s. and $46/57\ \%$ respectively for $0\nu\beta\beta$
decay.
We fit these histograms to phenomenological functions to extract analytic
signal shapes to use in our analysis. These functions are described in more
detail in Appendix B, the photo-peak is approximated via three Gaussian’s to
account for the slight non-Gaussianity of the line-shape (due to each
detector-dataset pair having a different resolution), this is discussed more
in Section IV.1. Each function is normalized to unity and then we model the MC
simulated spectrum as:

$f(E)=N_{\text{tot}}\sum_{i=1}^{N}p_{i}\cdot f_{i}(E),$ (5)

where the sum runs over the normalized functions $f_{i}$, $p_{i}$ is the
probability for an event to be in a given spectral feature $i$ of the signal
distribution for the considered category, for example in one of the gamma
peaks or in the continuum, and $N_{\text{tot}}$ is the total number of events.
We show two example fits (both $0_{1}^{+}$ signal) in Figure 6,
$2\nu\beta\beta$ category 2 and Figure 6, $0\nu\beta\beta$ category 2. We see
that these models describe the MC simulations very well, in particular the
slightly non-Gaussian photo-peak shape.

Figure 5: Parameterization of $2\nu\beta\beta$ to $0_{1}^{+}$ category 2
signal. The model consists of two photo-peaks, with their respective 100Mo
X-ray escape peaks, two smeared step functions modeling Compton interactions
and a linear background. The residual defined as
$(\text{value}-\text{model})/\text{error}$ is shown in the upper figure.

Figure 6: Parameterization of $0\nu\beta\beta$ decay to $0_{1}^{+}$ category 2
signal. The model contains a photo-peak, X-ray escape peak, and linear smeared
step functions. The residual defined as
$(\text{value}-\text{model})/\text{error}$ is shown in the upper figure.

These functions encode most of our systematic uncertainties, in particular the
resolution and peak position of Gaussian and normalization of the peaks which
can then be treated as nuisance parameters of the final fit.

##  IV Data analysis

For this search we use the full CUPID-Mo data corresponding to 2.71
kg$\times$yr exposure of LMO. We process our data using a C++ software, Diana
[58, 59], developed by the Cuoricino, CUORE, and CUPID-0 collaborations and
further developed by CUPID-Mo. Most of the data processing steps are the same
as described in [47, 31]. An optimal trigger is used to identify physics
events and an optimal filter which maximizes signal-to-noise ratio is used to
compute amplitudes (both LD and LMO detectors) [60]. Thermal gain
stabilization and calibration are performed using data collected with a Th/U
calibration source. For each LMO detector we associate ”side-channels”, as the
LDs which face this crystal. For $\mathcal{M}_{1}$ data (3 $0\nu\beta\beta$
decay categories) we use the same selection cuts as in [31]: a principal
component analysis (PCA) based pulse shape discrimination (PSD) cut [61], a
normalized light distance cut, delayed coincidence (DC) cuts to remove 222Rn,
228Th decay chain events and, muon veto anti-coincidence.
However, for coincidence events ($\mathcal{M}>1$) it is beneficial to modify
several of these steps. Random coincidences between non-causally related
events (hereafter accidentals) could provide a large possible background.
Similar to CUORE [58] and as described in [31] we adjust the time difference
of events to account for the characteristic detector response based time
offsets between detector pairs. This allows us to reduce the time window used
to define coincidences to $10\,\mathrm{ms}$ (from 100 ms in previous analyses
[47]) and thus reduce the accidental coincidence background.
In order to select a clean sample of higher multiplicity e.s. candidate
events, we make use of the dual read-out and implement a light yield based cut
for $\mathcal{M}>1$ events. Unlike the $\mathcal{M}_{1}$ analysis where this
cut primarily rejects $\alpha$ backgrounds, it is designed here to tag and
reject events where an energetic electron escapes an LMO crystal, punches
through an LD, and is stopped in the adjacent crystal within a tower. In
addition, this cut very efficiently tags events from a 60Co contamination
identified on one of the LDs where the $\beta$ particle typically leaves
significant excess energy in the LD. Developing a light cut to remove these
events is less straightforward than for $\mathcal{M}_{1}$ data, since a
coincidence between two vertically adjacent crystals is accompanied by
scintillation light signal that can be absorbed in the common, intermediate
LD. For each event we compute the expected scintillation light deposited on
the LD based on the values observed in $\mathcal{M}_{1}$ events (mainly
$2\nu\beta\beta$ decays to g.s.). This accounts for these multiple
contributions to light yield. We then normalize the light energies using a
conservative estimate of the LD resolution ($\sigma^{\text{LD}}_{i,s}$)
measured in high energy background data ($\mathcal{M}_{2}$ total energy
$>1500$ keV) as:

$n_{i,s}=\frac{E^{\text{LD}}_{i,s}-E^{\text{LD,exp}}_{i,s}}{\sigma^{\text{LD}}_{i,s}},$
(6)

where $i$ refers to the pulse index (one for $E_{1}$ etc.), $s$ is the side
channel (either 0 or 1), $E_{i,s}^{\text{LD}}$ is the measured LD energy,
while $E_{i,s}^{\text{LD,exp}}$ is the expected LD energy for $\gamma/\beta$
like event and $\sigma_{i,s}^{\text{LD}}$ is the predicted energy resolution
of the LD. We then place a cut of $|n_{i,s}|<10$ for each LD and side channel
in an $\mathcal{M}>1$ event. We place a cut of $-10<n_{i,s}<3$ for the LD with
the 60Co contamination. This contamination can also lead to background events
in other LMO detectors; to remove this we make a global LD anti-coincidence
cut. We remove any LMO event (excluding those who directly face this LD) with
a trigger on this LD with LD energy $>5$ keV within a $\pm 5$ ms window.
We also adjust the PCA based PSD cuts [47, 31] to place a cut on the shape of
each pulse in a higher multiplicity event. We require that each pulse has a
normalized PCA reconstruction error of less than 23 median absolute
deviations. This was optimized by comparing the efficiency of
$\mathcal{M}_{2}$ events $\varepsilon_{s}$, obtained from events summing to
the 40K 1461 keV peak and estimating the background efficiency,
$\varepsilon_{b}$, in a side-band $E_{2}\in(450$ – $520)$ keV. This sideband,
which is only used in optimizing the PCA cut, approximates the background for
the dominant $2\nu\beta\beta$ decay categories (2, 3 and 4) whilst not
including any signal peaks. We then maximize a figure of merit
$\varepsilon_{s}/\sqrt{\varepsilon_{b}}$ which is proportional to the
experimental sensitivity. For $\mathcal{M}_{1}$ data a cut of $<9$ on the PCA
normalized reconstruction error is used (as in [31]).
Finally, we employ a data blinding by adding simulated MC events directly to
the data files. The rate is chosen by sampling randomly from a uniform
distribution with range, $\Gamma\in\Gamma_{\text{measured}}\cdot(2,10)$, where
$\Gamma_{\text{measured}}$ is the current leading limit or measurement [2] and
$(2,10)$, is a uniform distribution between 2 and 10 chosen to ensure the
injected signal is significantly larger than any possible signal in the data.
The injected rate is hidden during the analysis of the blinded data. The
blinded data are used to optimize and test the Bayesian fitting routine and
prevent biasing our results.

###  IV.1 Energy resolution and linearity

We determine the response of our detector to a monochromatic energy deposit
using both calibration sources and $\gamma$ peaks from natural radioactivity
using the same procedure described in detail in [31]. In particular, a
simultaneous fit, over each detector in each dataset, of the 2615 keV peak in
calibration data is used to extract the resolution of each detector-dataset
pair at 2615 keV $\sigma_{\text{chan,dataset}}^{2615}$.
We model the line-shape of our peaks in $\mathcal{M}_{1}$ physics data as an
exposure weighted sum of Gaussians for each channel-dataset pair. The
individual Gaussians have a common mean $\mu$ but a resolution which is a
product of $\sigma_{\text{chan,dataset}}^{2615}$ and a global scaling $R(E)$.
We fit this line-shape model to our peaks to extract $R(E)$ and $\mu(E)$ for
each peak in physics data and parameterize the energy dependence as
$R(E)=\sqrt{r_{0}^{2}+r_{1}E}$ and $\mu(E)=b_{0}+b_{1}E+b_{2}E^{2}$. The
resulting parametrization is used in the MC simulations to model the detector
response. In particular, we account for the uncertainty on the parameters of
these functions as systematic uncertainties (Section VI).

###  IV.2 Cut efficiencies

To measure $2\nu\beta\beta/0\nu\beta\beta$ decay rates it is necessary to
determine the analysis efficiency, or probability that a signal event will
pass all selection cuts. We employ the same strategy that was utilised in the
$\mathcal{M}_{1}$ 0$\nu\beta\beta$ analysis [31] in order to compute these. We
evaluate the pileup efficiency, or probability that a signal event will not
have another trigger in the same waveform using random triggers. Several other
cuts (multiplicity selection, muon veto cut, and LD anti-coincidence) induce
dead times in one or more detectors. For these we evaluate the efficiencies
using 210Po Q-value peak events by counting the number of events passing the
cuts with energy in $(5407\pm 50)$ keV. These are a proxy for physical
$\mathcal{M}_{1}$ events due to the high energy and modular structure meaning
$\alpha$ particles can only deposit energy in one crystal.
Finally, for our pulse shape and light yield cuts the efficiency is evaluated
using $\gamma$ peaks in $\mathcal{M}_{1}$ and summed $\mathcal{M}_{2}$ energy
which are a clean sample of signal like events. For each prominent $\gamma$
peak, we fit the energy distribution of events passing and failing each cut to
a Gaussian plus linear background. From the measured number of events
$N_{\text{pass}},N_{\text{fail}}$ (events in the Gaussian peak only for the
two fits) we compute the efficiency
$\varepsilon=N_{\text{pass}}/(N_{\text{pass}}+N_{\text{fail}})$. We estimate
numerically the uncertainty on $\varepsilon$ by sampling from the measured
uncertainties on $N_{\text{pass}},N_{\text{fail}}$. We observe an energy
independent efficiency for all of our cuts (shown in col. 2 of Table 4). The
energy range for the $\mathcal{M}_{1}$ 3034 keV category of $0\nu\beta\beta$
decay lies outside the range of our $\gamma$ peaks. Similar to [31] we use a
1st order polynomial to extrapolate the PCA and light yield cut efficiency to
3034 keV, accounting for the possibility of a slight energy dependence of the
cut (shown in col. 3 of Table 4). We summarise the cut efficiencies in Table
4, and use them as nuisance parameters with Gaussian priors in our Bayesian
analysis.

Table 4: Efficiency for analysis cuts, the energy independent efficiencies are computed assuming the cut efficiency does not depend on energy. For the extrapolated efficiency in $\mathcal{M}_{1}$ data at 3034 keV, we fit the efficiency as a function of energy to a first order polynomial and extrapolate this to 3034 keV, this is only relevant for the $\mathcal{M}_{1}$ signatures with 3034 keV energy.  Cut | Efficiency Energy Independent [%] | Efficiency Extrapolated [%] | Method of Evaluation
---|---|---|---
Pileup | $95.7\pm 1.0$ | - | Noise events
Multiplicity | $99.55\pm 0.07$ | - |  210Po
LD Anti-coincidence | $99.976\pm 0.017$ | - |  210Po
Delayed coincidence | $99.16\pm 0.01$ | - |  210Po
Muon Veto | $99.62\pm 0.07$ | - |  210Po
Light Cut ($\mathcal{M}_{1}$) | $99.4\pm 0.4$ | $99.7\pm 0.8$ | Gamma Peaks ($\mathcal{M}_{1}$)
Light Cut $\mathcal{M}>1$ | $97.7\pm 1.8$ | - | Gamma Peaks ($\mathcal{M}_{2}$)
PCA $<23$ ($\mathcal{M}_{1}$) | $99.2\pm 0.3$ | - | Gamma Peaks ($\mathcal{M}_{1}$)
PCA $<9$ $(\mathcal{M}1)$ | $95.2\pm 0.5$ | $94.3\pm 1.5$ | Gamma Peaks $(\mathcal{M}_{1})$

##  V Bayesian analysis

We perform a Bayesian analysis to extract decay rates. In particular, we use
an extended unbinned maximum likelihood fit implemented using the Markov Chain
Monte Carlo (MCMC) of the Bayesian Analysis Toolkit (BAT) [62] for all
categories except for two $0\nu\beta\beta$ decay categories, for which we use
a binned fit due to their exceptionally high statistics. We use a fit of the
summed data of all 19 detectors and 9 datasets.
Two separate fits are created, one to extract the decay rates of
$2\nu\beta\beta$ decay and another for the $0\nu\beta\beta$ decay. Both fits
are implemented in the same framework. In each category (index $r$) we model
our experimental data as:

$\displaystyle f_{r}(E)$ $\displaystyle=\varepsilon_{a,r}n_{2^{+}}\cdot
f_{2^{+},r}(E)+\varepsilon_{a,r}n_{0^{+}}\cdot f_{0^{+},r}(E)$ (7)
$\displaystyle+n_{r,\text{bkg}}\cdot
f_{r,\text{bkg}}(E)+\sum_{p=1}^{5}n_{r,p}f_{p}(E)$

where

  * •

$n_{0^{+}/2^{+}}$ are the observed number of counts to $0_{1}^{+}/2_{1}^{+}$
excited states.

  * •

$f_{0^{+}/2^{+},r}(E)$ are the phenomenological functions describing the
signal shape of category $r$, normalized to unity over the sum of all
categories.

  * •

$\varepsilon_{a,r}$ is the analysis efficiency for category $r$ (described in
Section VI).

  * •

$f_{r,\text{bkg}}(E)$ is a function describing the background in category $r$,
either exponential or flat for categories with low statistics.

  * •

$n_{r,p}$ are the number of background counts from known $\gamma$ lines (a
subset of 511, 583, 609, 2448 and 2505 keV) with index $p$ and $f_{p}(E)$ is
the model of this spectral shape, a single Gaussian.

We show in Fig. 7 the fit of the $2\nu\beta\beta$ decay category 3 and the
contribution from the exponential background, signal and background peaks.

Figure 7: Fit of the $2\nu\beta\beta$ decay category 3. The experimental data
are binned for visualization and we show the best fit reproduction of our
model and the contributions from the signal peaks, background peaks and the
exponential background.

Our likelihood function is:

$\displaystyle\text{log}(\mathcal{L})=$ $\displaystyle\sum_{n=1}^{N}\text{log
}(f_{r(n)}(E_{n}))$ (8) $\displaystyle+$
$\displaystyle\sum_{r=1}^{C}-\lambda_{r}+n_{r}\text{log}(\lambda_{r})-\text{log}(n_{r}!)$
$\displaystyle+$
$\displaystyle\sum_{b=1}^{B}\sum_{i=1}^{N_{\text{bin}}}-\lambda_{i}+n_{i,b}\cdot\text{log(}\lambda_{i,b})-\text{log}(n_{i,b}!),$

where the first sum is over the events in the experimental data (unbinned
categories), $r(n)$ is the category of event $n$, $C$ is the number of
categories, $\lambda_{r}$ is the total predicted number of events in category
$r$, while $n_{r}$ is the observed number of events. The last sums are a
binned likelihood where the sum $b$ runs over the $B$ categories which use a
binned likelihood. These are the categories 0 and 1 for the $0\nu\beta\beta$
decay where a large number of events make an unbinned fit computationally
unfeasible. $\lambda_{i}$ is the expectation value for bin $i$ and $n_{i,b}$
is the number of events in experimental data. We use 0.5 keV bins for these
fits, much smaller than the energy resolution over the full spectral range.
We use BAT to sample the full posterior distribution of the parameters of
interest $\vec{\theta}=(n_{0_{1}^{+}},n_{2_{1}^{+}})$ and nuisance parameters
$\vec{\nu}$, which include all of the background components (both the
exponential/flat background and the $\gamma$ peaks) and also parameters of
systematic uncertainties (see Section VI)

$p(\vec{\theta},\vec{\nu}|\mathcal{D})=\frac{P(\mathcal{D}|\vec{\theta},\vec{\nu})P(\vec{\theta},\vec{\nu})}{P(\mathcal{D})},$
(9)

where $\mathcal{D}$ represents the data. We use uninformative (flat) priors on
all background model components and parameters of interest $n_{2^{+}/0^{+}}$.
We define observables for the decay rates of $0_{1}^{+}/2_{1}^{+}$ decay
($\Gamma_{0^{+}/2^{+}}$)

$\displaystyle\Gamma_{0^{+}/2^{+}}$ $\displaystyle=\frac{n_{0^{+}/2^{+}}\cdot
W}{\varepsilon_{\text{cont},0^{+}/2^{+}}\cdot\eta\cdot Mt\cdot N_{a}\cdot
1000},$ (10)

where $W$ is the molecular mass for enriched Li${}_{2}\,^{100}$MoO4, $N_{a}$
is the Avogadro constant, $Mt$ is the exposure (in
$\mathrm{kg}\times\mathrm{yr}$), $\eta$ is the isotopic abundance of 100Mo
and, $\varepsilon_{\text{cont},0^{+}/2^{+}}$ is the containment efficiency for
$0_{1}^{+}/2_{1}^{+}$ signal, i.e. the total probability a simulated event is
in one of the categories. The half-life is then given by $\ln{(2)}/\Gamma$.
For each step of the Markov chains the values of $\Gamma_{0^{+}/2^{+}}$ are
stored and used to compute the marginalized posterior distributions. We
include all of our systematic uncertainties as nuisance parameters in our
analysis as described in Section VI.

##  VI Systematic uncertainties

###  VI.1 Energy resolution and bias

We account for the uncertainty in the energy resolution of our peaks. Both the
terms $r_{1}$ and $r_{0}$ in the resolution scaling function ($R(E)$, see
Section IV.1) are given Gaussian priors based on the measured uncertainty from
the fit to $\gamma$ peak resolution. A multivariate prior is not necessary
since the correlation is very low ($-0.2$%) because the $r_{0}$ term is fixed
very well by noise events (random triggers used to estimate energy resolution
at 0 keV). At each stage in the Markov chain this resolution is applied to all
the signal components and background peaks, by adjusting the resolution of the
functions $f_{2^{+}/0^{+},r}(E),f_{\text{r,bkg}}(E)$.

We also account for the uncertainty in the energy scale. Our energy bias is
parameterized as a second order polynomial as described earlier in Section
IV.1. We adjust the mean position of all peaks in our signal and background
template functions ($f_{0^{+}/2^{+}},f_{\text{bkg}}$) by this bias. We add a
nuisance parameter to the model accounting for the uncertainty on the bias,
which is given a Gaussian prior with $0.07$ keV $\sigma$, and vary the
position of each peak (in a correlated way) by this amount. For the
$0\nu\beta\beta$ decay analysis the peaks cover a wide range in energy and we
take a conservative approach varying the position of all the peaks by the
largest uncertainty of $0.3\ \mathrm{keV}$ (from 3034 keV).

###  VI.2 Analysis efficiency

We account for the uncertainty in the analysis efficiency, or the probability
that a signal like event would pass all cuts. These cut efficiencies are
included as nuisance parameters in our model which are given Gaussian priors
based on the estimated uncertainty (see Section IV.2). A nuisance parameter is
included for each cut, including two separate parameters for the constant and
extrapolated light cut, and PCA$<9$ cut. We then predict the efficiency in
category $r$ as:

$\varepsilon_{r}=\prod_{c=1}^{11}\varepsilon_{c}^{p(\mathcal{M}(r),c)},$ (11)

where the product $c$ runs over all cuts (both energy independent and
extrapolated, see Table 4) and $p(\mathcal{M}(r),c)$ is a power which
represents how many times a cut was applied to an event. For example, for PSD
the cut is applied for both pulses so the efficiency is raised to the power of
the multiplicity.

###  VI.3 Containment Efficiency

The final systematic uncertainty we account for is the containment efficiency
from MC simulations. In particular, this is related to the accuracy of our
Geant4 MC simulations which can broadly be divided into two parts: the
accuracy of the simulated geometry and the implemented physics process. For
the simulated geometry we vary the amounts of the three main materials in the
experimental setup (see Figure 2):

  * •

The dimensions of the LMO crystals (and therefore density since the mass is
very well known);

  * •

The thickness of the copper holders;

  * •

The thickness of the Ge LDs.

We run a set of simulations varying the LMO dimensions (by $\pm 100/200\
\mathrm{\mu m}$ or around $\pm$0.2/0.4%), the Cu holder thickness (by $\pm
100/200\ \mathrm{\mu m}$ or $\pm$2.5/5%) and the Ge LD thickness (by $\pm
10/20\ \mathrm{\mu m}$ or around $\pm$6/12%). We also account for the possible
inaccuracy of Geant4 physics cross sections by running MC simulations varying
both the Compton effect and photo-effect cross sections by $\pm 5$% which is a
conservative estimate of Geant4 cross section accuracy [63, 64]. For each
systematic effect (Cu thickness, LMO dimension, Ge thickness, Compton or
photo-effect cross section scale) we introduce a parameter $\nu_{i}$. For
every effect, decay, category (index $r$) and peak (index $p$) we parameterize
the variation in peak containment efficiency. This is the product of the
overall containment efficiency and the Gaussian peak probability (see Section
III.4). We use first order polynomial fits to obtain
$d\varepsilon_{\text{cont},r,p}/d\nu_{i}$. We show the results of all these
tests in Appendix C. We observe that the most significant effects are the
photo-effect cross section (higher scale tends to increase containment) and
the Compton cross section (higher tends to decrease containment). The
geometrical uncertainties provide a much smaller effect. Effects where the
variation in containment efficiency is consistent with 0 (within $2\sigma$)
are not included. Based on the parameters $\vec{\nu}$ we predict the
efficiency for category $r$ and peak $p$ as:

$\varepsilon_{\text{cont},r,p}=\varepsilon^{0}_{\text{cont},r,p}+\sum_{c=1}^{5}\nu_{c}\cdot
d\varepsilon_{\text{cont},r,p}/d(\nu_{c}),$ (12)

where the sum $c$ runs over systematic effects, and
$\varepsilon^{0}_{\text{cont},r,p}$ is the default efficiency (for the MC
simulation without varying any parameters). In our fit, the nuisance
parameters accounting for uncertainties of copper thickness and LMO dimension
are given Gaussian priors with a standard deviation of 100$\ \mathrm{\mu m}$;
for the Ge thickness we use 20$\,\mathrm{\mu m}$, and we use $\ 5\%$ for the
Compton and photo-effect physics cross sections. At each step of the Markov
chain the amplitude of each peak is adjusted according to these parameters.
The normalization of the fit functions are adjusted accordingly and we
recompute the containment efficiency for both $0_{1}^{+}/2_{1}^{+}$ signals.

###  VI.4 Choice of $2\nu\beta\beta$ decay to $0_{1}^{+}$ state model

The final systematic uncertainty we consider is the $2\nu\beta\beta$ to
$0_{1}^{+}$ state decay model. By default we use the SSD model, however we
also run MC simulations and a fit using the HSD model. We then consider the
difference between these two results as a systematic on the $2\nu\beta\beta$
decay to $0_{1}^{+}$ state decay rate.

##  VII Results

###  VII.1 Toy Monte-Carlo sensitivity and bias

To validate our fitting routine and predict the median exclusion sensitivity
we use an ensemble of pseudo-experiments (toy MC). We generate 2000 toy
datasets for both $0\nu\beta\beta$ and $2\nu\beta\beta$ decay analysis based
on the background parameters from the fit to data (see Sections
VII.2.1,VII.2.2). We assume the $2\nu\beta\beta\rightarrow 0_{1}^{+}$ signal
from the best fit to data and zero signal of $2\nu\beta\beta\rightarrow
2_{1}^{+}$ or $0\nu\beta\beta$ decay modes. We sample from a Poisson
distribution for the number of counts of each component in each toy
experiment. We fit each dataset to determine the distribution of possible
limits and therefore the median sensitivity for each decay. This is shown in
Figures 9 and 9. We extract median exclusion sensitivities, $\hat{T}$, at 90%
credibility interval (c.i.) for the three decays which have not yet been
measured of:

$\displaystyle\hat{T}_{1/2}^{2\nu\rightarrow 2_{1}^{+}}$
$\displaystyle=3.8\times 10^{21}\ \mathrm{yr},$ (13)
$\displaystyle\hat{T}_{1/2}^{0\nu\rightarrow 2_{1}^{+}}$
$\displaystyle=3.1\times 10^{23}\ \mathrm{yr},$ (14)
$\displaystyle\hat{T}_{1/2}^{0\nu\rightarrow 0_{1}^{+}}$
$\displaystyle=9.8\times 10^{22}\ \mathrm{yr}.$ (15)

We also use these toys to verify the fitting procedure is not biased. We fit
our toy datasets for $2\nu\beta\beta$ decay analysis with both a model
including the $2_{1}^{+}$ signal called $H_{0}=H(0_{1}^{+}+2_{1}^{+}+B)$ and
one including only $0_{1}^{+}$ signal, $H_{1}=H(0_{1}^{+}+B)$. We compute

$r_{i}=\frac{n_{0^{+}}^{i}-\bar{n_{0^{+}}}}{\sigma(n_{0^{+}})},$ (16)

where $\bar{n_{0^{+}}}$ is the marginalized mode from the fit to real data
used to generate the toys, $n_{0^{+}}^{i}$ is the marginalized mode in toy
$i$, and $\sigma(n_{0^{+}})$ is the estimated error (central 68% c.i.). The
distribution of $r_{i}$ for both models is shown in Fig. 10. We observe a
clear bias in $r_{i}$ for the $H_{0}$ model with a mean of $-0.58\sigma$,
while $r_{i}$ is unbiased for the $H_{1}$ model (mean of $-0.005\sigma$).
Hence, we decide to use the $H_{1}$ model for the case that no evidence of a
decay to $2_{1}^{+}$ state is found. This bias is due to non-negative Bayesian
priors on the rate of $2_{1}^{+}$ decay which allow the rate of $2_{1}^{+}$ to
fluctuate only up. Since the $0_{1}^{+}$ and $2_{1}^{+}$ rates are anti-
correlated this causes the $0_{1}^{+}$ rate to be biased.

Figure 8: Sensitivity for $0\nu\beta\beta$ decay to both $0_{1}^{+}$ and
$2_{1}^{+}$ excited states obtained from an ensemble of 2000 toy experiments.
We show the median sensitivity (dashed line) and the limit obtained on real
data (solid lines). In both cases the limits on real data are consistent with
those expected from toy experiments.

Figure 9: Sensitivity for $2\nu\beta\beta$ decay to $2_{1}^{+}$ excited state
obtained from an ensemble of 2000 toy experiments. As with $0\nu\beta\beta$
decay we show the median sensitivity (dashed line) and the limit obtained on
real data (solid line). Again the limit on real data is consistent with the
toy experiments. Figure 10: The distribution of $r_{i}$ (see Eq. 16) for both
$H_{0}$ and $H_{1}$ models for our fits to toy data. We observe a clear bias
for the $H_{0}$ model while the $H_{1}$ model is unbiased.

Figure 11: Posterior probability distributions for the decay rates,
$\Gamma^{2\nu}$ to $2_{1}^{+}$ state (top left) and $0_{1}^{+}$ state (top
right). The mode for $2_{1}^{+}$ decay is at zero rate and so we find no
evidence for this decay and we extract the measurement of $0_{1}^{+}$ mode
using the $H(0_{1}^{+}+B)$ model. Bottom left/right are the posteriors of
$\Gamma^{0\nu}$ to $2_{1}^{+}$ and $0_{1}^{+}$ state respectively, again
showing no evidence of the decay.

###  VII.2 Fit of $2\nu\beta\beta$ with $H_{0}=H(0_{1}^{+}2_{1}^{+}+B)$ model

We run the $2\nu\beta\beta$ decay fit with both $2\nu\beta\beta$ decay to
$2_{1}^{+}$ and $0_{1}^{+}$ e.s. signals and background components, i.e. the
$H_{0}$ model. We can use this fit to determine if we observe evidence for the
decay to $2_{1}^{+}$ state and to set a limit if no evidence is found. We
consider first a fit including all possible background peaks, the maximal
model. We then repeat the fit removing any background peak where the central
68% confidence interval contains 0 counts and replacing the exponential
background with a constant if the slope is compatible with 0 (within 1
$\sigma$). We call this the minimal model and it is this model we use for
statistical inference. We show the marginalized posterior distribution of the
decay rate $\Gamma$ for $2\nu\beta\beta$ decay to $2_{1}^{+}$ state (including
all systematics) in Fig. 11 (left). We see that the mode is at 0 rate and
therefore find no evidence for $2\nu\beta\beta$ to $2_{1}^{+}$ state. We
correspondingly set a limit (including all systematics) of:

$T_{1/2}^{2\nu\rightarrow 2_{1}^{+}}>4.4\times 10^{21}\ \mathrm{yr}\
\text{(90\% c.i.)}$ (17)

This is the most stringent constraint on this process in 100Mo, an improvement
of $\sim 80\%$ over the previous constraint, $2.5\times 10^{21}\,\mathrm{yr}$
[42].

####  VII.2.1 Fit of $2\nu\beta\beta$ with $H_{1}=H(0_{1}^{+}+B)$ model

Since we find no evidence of the $2\nu\beta\beta$ decay to $2_{1}^{+}$ state
and toy experiments show that including this parameter in our fit would bias
our measurement, we run the $2\nu\beta\beta$ decay $H_{1}=H(0_{1}^{+}+B)$ fit
without the $2_{1}^{+}$ decay contribution. We find that this model is able to
describe our data in all seven categories very well (see Fig. 7 for an example
and all fits in Appendix D). We verify the goodness of fit using our ensemble
of toy experiments. For each fit we extract the global mode or the set of
parameters $\vec{\theta}$ with the maximum probability. We use the posterior
probability of these parameters:

$k=p(\vec{\theta}|\mathcal{D}),$ (18)

as a test statistic to quantify the quality of the fit by comparing the
probability distribution obtained in toy experiments $p(k^{{}^{\prime}})$ to
the fit on real data. We extract the $p$-value:

$p=\int_{0}^{k}p(k^{{}^{\prime}})dk^{{}^{\prime}}=0.49$ (19)

indicating that the model describes the data well.
We extract the posterior distribution on $\Gamma^{2\nu\rightarrow 0_{1}^{+}}$
from the fit on the data and we extract the central 68.3% c.i. on the decay
rate as:

$\Gamma_{2\nu\rightarrow 0_{1}^{+}}=(9.2\pm 1.0\ \text{(stat.+syst.)})\times
10^{-22}\,\mathrm{yr^{-1}}.$ (20)

As the systematics are allowed to float freely in our fit this error is the
sum of the statistical and systematic components. To evaluate the statistical
error, we repeat the fit fixing all nuisance parameters connected to
systematics, obtaining:

$\Gamma_{2\nu\rightarrow 0_{1}^{+}}=(9.1\pm 1.0\ \text{(stat.)})\times
10^{-22}\,\mathrm{yr^{-1}}.$ (21)

A difference in the error compared to the statistical plus systematic fit is
only observed after the first digit. We also run a fit using the HSD model,
which leads to:

$\Gamma_{2\nu\rightarrow 0_{1}^{+}}^{\text{HSD}}=(8.9\pm 1.0\
\text{(stat.+syst.)})\times 10^{-22}\,\mathrm{yr^{-1}},$ (22)

a $4\%$ decrease in the decay rate. Under the assumption that statistical and
systematic errors add in quadrature we obtain:

$\Gamma_{2\nu\rightarrow 0_{1}^{+}}=(9.2\pm 1.0\ \text{(stat.)}\
^{+0.3}_{-0.5}\ \text{(syst.)})\times 10^{-22}\,\mathrm{yr^{-1}},$ (23)

or converting to half life:

$T_{1/2}^{2\nu\rightarrow 0_{1}^{+}}=(7.5\pm 0.8\ \text{(stat.)}\
^{+0.4}_{-0.3}\ \text{(syst.)})\times 10^{20}\ \mathrm{yr}.$ (24)

This is a new independent measurement of this decay, with total uncertainty
consistent with the previous leading measurement [42].

####  VII.2.2 $0\nu\beta\beta$ decay fit

From the fit for $0\nu\beta\beta$ decay we find no evidence of either the
decay to $0_{1}^{+}/2_{1}^{+}$ excited states. The best fit reproductions of
the experimental data are shown in Appendix D. The posterior distributions of
$\Gamma^{0\nu\rightarrow 2_{1}^{+}/0_{1}^{+}}$ are shown in Fig. 11 (bottom
left and right), this leads to limits (both at 90% c.i. including all
systematic uncertainties) of:

$\displaystyle T_{1/2}^{0\nu\rightarrow 2_{1}^{+}}>2.1\times 10^{23}\
\mathrm{yr}\ \text{(90\% c.i.)},$ (25) $\displaystyle T_{1/2}^{0\nu\rightarrow
0_{1}^{+}}>1.2\times 10^{23}\ \mathrm{yr}\ \text{(90\% c.i.)}.$ (26)

These are new leading limits on these processes, a factor of 1.3 stronger than
previous limits from [39] in both cases, despite a factor of $\sim 6$ lower
exposure of 100Mo.

####  VII.2.3 $2\nu\beta\beta$ decay spectral shape

The CUPID-Mo source equals detector geometry also allows us to investigate the
$2\nu\beta\beta\rightarrow 0_{1}^{+}$ spectral shape. This provides a concrete
demonstration of a method for this analysis which can be applied to future
experiments with large statistics. Our analysis described in Section III does
not reconstruct the $\beta\beta$ spectral shape directly, so we develop a
procedure to extract this shape. For simplicity this analysis is only applied
to $\mathcal{M}_{2}$ data which is the main contribution to the sensitivity.

We select all events that have either $E_{1}$ or $E_{2}$ in the range
$[539.5-1.75,539.5+1.75]$ keV and $[590.8-1.75,590.8+1.75]$ keV. This choice
of the interval is the FWHM energy resolution, roughly optimal for maximising
the sensitivity $S/\sqrt{B}$. We note there is some ambiguity in this
selection when $E_{1}\sim E_{2}$ (when the $\beta\beta$ has the same energy as
the $\gamma$) however the excellent CUPID-Mo energy resolution makes this
negligible. We then define $E_{\beta\beta,540},E_{\beta\beta,591}$ as the
energy which does not contain the peak for both 540 and 591 keV peaks
respectively and we construct histograms of these energies. This procedure is
applied to both the data, signal MC simulations, and to the preliminary
background model reconstruction of the data.

To compare the SSD and HSD models we use a simultaneous binned Bayesian
likelihood fit of the two $\beta\beta$ spectra. We use three components;
background, $0_{1}^{+}$ signal and $2_{1}^{+}$ signal and 50 keV bins. Using
the SSD model, this fit reconstructs the half-life as:

$\displaystyle T_{1/2}^{2\nu\rightarrow 0_{1}^{+}}$ $\displaystyle=(7.3\pm
1.0\ \text{(stat.)})\times 10^{20}\,\mathrm{yr},$ (27) $\displaystyle
T_{1/2}^{2\nu\rightarrow 2_{1}^{+}}$ $\displaystyle>3.6\times
10^{21}\,\mathrm{yr}.$ (28)

These are consistent with the values from Sections VII.2,VII.2.1 and should be
considered a cross check of this more robust analysis which does not depend on
the quality of the background model. We show this fit reproduction in Fig. 12.
Repeating the fit using MC simulations obtained with the HSD model, we extract
the evidence for both models $p(\mathcal{D}|\text{SSD/HSD})$ and therefore
probabilities for the HSD or SSD models of:

$p(\text{SSD}|\mathcal{D})=0.51,\ p(\text{HSD}|\mathcal{D})=0.49.$ (29)

This indicates that the CUPID-Mo data is not able to differentiate between the
two models. However, this analysis provides a method which could be used to
differentiate these two models in a future experiment with larger statistics.

Figure 12: Plots showing the $2\nu\beta\beta$ decay to e.s. $\beta\beta$
spectral shape for 540 keV peak (left) and 591 keV peak (right). We overlay
the best fit prediction. We find that in the region $>600$ keV, the spectrum
is dominated by $2\nu\beta\beta\rightarrow 0_{1}^{+}$.

##  VIII Discussion

####  VIII.0.1 Matrix element for $2\nu\beta\beta$ decay to $0_{1}^{+}$

From our measurement of $2\nu\beta\beta\rightarrow 0_{1}^{+}$ excited state we
can measure experimentally the nuclear matrix element for this process based
on:

$\Gamma^{2\nu}/\ln{(2)}=G_{2\nu}\cdot|M^{\text{eff}}_{2\nu}|^{2}.$ (30)

The phase space factor is given by $(60.55$ – $65.18)\times 10^{-21}\
\mathrm{1/yr}$ depending on whether the HSD or SSD model is used [65].
Therefore the matrix element (assuming the SSD model) is given by:

$\displaystyle M^{\text{eff}}_{2\nu\rightarrow 0_{1}^{+}}$
$\displaystyle=\sqrt{\frac{\Gamma}{(\ln{(2)}\cdot
G)}}\pm\frac{1}{2\sqrt{\ln{(2)}G}}\frac{\Delta\Gamma}{\sqrt{\Gamma}}$ (31)
$\displaystyle=0.143\pm 0.008\ \text{(stat.+syst.)}.$

We compare the theoretical predictions of dimensionless
$M^{\text{eff}}_{2\nu}=g_{A}^{2}M_{2\nu}$ calculated assuming and unquenched
value of $g_{A}=1.27$. These are 0.395, 0.595, 0.185 for the shell model [66],
microscopic interacting Boson model [67], and quasi-particle random phase
approximation [68], respectively. This shows the decay rate is quenched
relative to theoretical predictions.

####  VIII.0.2 Effective Majorana neutrino mass, $m_{\beta\beta}$

The limits on $0\nu\beta\beta\rightarrow 0_{1}^{+}$, despite the lower phase
space, can be used to set a limit on $m_{\beta\beta}$. We use the phase space
factor ($3.162\times 10^{-15}\,\mathrm{yr^{-1}}$ from [65]) and the NMEs from
[66, 67] to obtain:

$m_{\beta\beta}<(11-15)\,\mathrm{eV},$ (32)

depending on the NMEs used. Whilst this is still several orders of magnitude
above the constrains from the g.s. decay this can be improved significantly in
future for a large experiment with minimized dead material.

####  VIII.0.3 Bosonic neutrinos

Double beta decays to $J\neq 0$ excited states can be sensitive to a Bosonic
contribution to the neutrino wave-function [33]. In particular, under the SSD
hypothesis, the predicted half-lives for $2\nu\beta\beta$ to the $2_{1}^{+}$
e.s. are $2.4\times 10^{22},\ 1.7\times 10^{23}\,\mathrm{yr}$ for Bosonic and
Fermionic neutrinos respectively. The current limit of $4.4\times
10^{21}\,\mathrm{yr}$ is still below these predictions.

####  VIII.0.4 SSD vs HSD models

Both NEMO-3 and CUPID-Mo have demonstrated that SSD hypothesis describes the
experimental data very well for the $2\nu\beta\beta$ to $0^{+}_{\text{g.s.}}$
decay in 100Mo [56, 57]. However, for the decay to $0_{1}^{+}$ state the
mechanism (SSD or HSD/closure approximation) is not known. CUPID-Mo is only
the second experiment to reconstruct the $2\nu\beta\beta$ to $0_{1}^{+}$
$\beta\beta$ spectral shape. However, the present level of statistics is
insufficient to distinguish the two modes. The excellent energy resolution and
low background rates mean this will be possible in a future experiment. In
this paper, we demonstrate a method to quantify numerically this and obtained
the first result on the compatibility of these two models with data.

##  IX Conclusion

In this paper we have presented a new analysis of $0\nu\beta\beta$ and
$2\nu\beta\beta$ transitions of 100Mo to the first two ($2_{1}^{+}/0_{1}^{+}$)
excited states of 100Ru using the full exposure of CUPID-Mo. This analysis
exploits the information available for a source equals detector experiment,
where both $\beta$ and $\gamma$’s can be measured. A measurement of 100Mo
$T_{1/2}^{2\nu\rightarrow 0_{1}^{+}}$ was obtained:

$T_{1/2}^{2\nu\rightarrow 0_{1}^{+}}=(7.5\pm 0.8\ \text{(stat.)}\
^{+0.4}_{-0.3}\ \text{(syst.)})\times 10^{20}\ \mathrm{yr}.$ (33)

For the other three decay modes, no evidence was found and we extract the
limits:

$\displaystyle T_{1/2}^{2\nu\rightarrow 2_{1}^{+}}$ $\displaystyle>4.4\times
10^{21}\ \mathrm{yr}\ \text{(90\% c.i.)},$ (34) $\displaystyle
T_{1/2}^{0\nu\rightarrow 0_{1}^{+}}$ $\displaystyle>1.2\times 10^{23}\
\mathrm{yr}\ \text{(90\% c.i.)},$ (35) $\displaystyle T_{1/2}^{0\nu\rightarrow
2_{1}^{+}}$ $\displaystyle>2.1\times 10^{23}\ \mathrm{yr}\ \text{(90\%
c.i.)}.$ (36)

These are the leading limits on these processes. The sensitivity was limited
by the small size of the array and large amount of dead material which results
in a low containment efficiency. Future experiments such as CUPID [13] or
CROSS [69] will feature a more tightly packed array, much lower amounts of
dead material, and a much larger exposure. This will lead to a significantly
improved sensitivity with a precision measurement of $2\nu\beta\beta$ decay to
$0_{1}^{+}$ state, distinction between the SSD and HSD mechanisms of the decay
and the possibility to measure the $2\nu\beta\beta$ decay to $2_{1}^{+}$
state.

##  X Acknowledgements

This work has been partially performed in the framework of the LUMINEU
program, a project funded by the Agence Nationale de la Recherche (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 CEA/SPEC for their valuable contribution
in the detector conception and of LAL (now IJCLab) for the detector holders
fabrication. 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”. A.S.
Barabash, S.I. Konovalov, I.M. Makarov, V.N. Shlegel and V.I. Umatov were
supported by the Russian Science Foundation under grant No. 18-12-00003. We
acknowledge 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
the Agence Nationale de la Recherche (ANR), France. 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. J.
Kotila is supported by Academy of Finland (Grant Nos. 3314733, 320062, 345869)
This research used resources of the National Energy Research Scientific
Computing Center (NERSC). This work makes use of the Diana data analysis
software which has been 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 . Majority of the work
described here was completed before February 24, 2022.

## References

  * Saakyan [2013] R. Saakyan, Two-Neutrino Double-Beta Decay, Annu. Rev. Nucl. Part. Sci. 63, 503 (2013).
  * Barabash [2020] A. S. Barabash, Precise Half-Life Values for Two-Neutrino Double-$\beta$ Decay: 2020 Review, Universe 6, 159 (2020).
  * Barabash [2017] A. S. Barabash, Double beta decay to the excited states: Review, AIP Conf. Proc. 1894, 020002 (2017), arXiv:1709.06890 [nucl-ex] .
  * Dell’Oro  _et al._ [2016] S. Dell’Oro, S. Marcocci, M. Viel, and F. Vissani, Neutrinoless Double Beta Decay: 2015 Review, Advances in High Energy Physics 2016, 2162659 (2016).
  * Schechter and Valle [1982] J. Schechter and J. W. F. Valle, Neutrinoless double-$\beta$ decay in SU(2)×U(1) theories, Phys. Rev. D 25, 2951 (1982).
  * Furry [1939] W. H. Furry, On Transition Probabilities in Double Beta-Disintegration, Phys. Rev. 56, 1184 (1939).
  * Bilenky and Giunti [2015] S. M. Bilenky and C. Giunti, Neutrinoless double-beta decay: A probe of physics beyond the Standard Model, International Journal of Modern Physics A 30, 1530001 (2015).
  * Dolinski  _et al._ [2019] M. J. Dolinski, A. W. Poon, and W. Rodejohann, Neutrinoless Double-Beta Decay: Status and Prospects, Annu. Rev. Nucl. Part. Sci. 69, 219 (2019).
  * Fukugita and Yanagida [1986] M. Fukugita and T. Yanagida, Barygenesis without grand unification, Phys. Lett. B 174, 45 (1986).
  * Davidson  _et al._ [2008] S. Davidson, E. Nardi, and Y. Nir, Leptogenesis, Physics Reports 466, 105 (2008).
  * Deppisch  _et al._ [2018] F. F. Deppisch _et al._ , Neutrinoless Double Beta Decay and the Baryon Asymmetry of the Universe, Phys. Rev. D 98, 055029 (2018), arXiv:1711.10432 [hep-ph] .
  * Adams  _et al._ [2022] D. Q. Adams _et al._ (CUORE), Search for Majorana neutrinos exploiting millikelvin cryogenics with CUORE, Nature 604, 53 (2022), arXiv:2104.06906 [nucl-ex] .
  * Armstrong  _et al._ [2019] W. R. Armstrong _et al._ (CUPID Interest Group), CUPID pre-CDR, arXiv:1907.09376 (2019).
  * Armengaud  _et al._ [2020a] E. Armengaud _et al._ (CUPID-Mo), The CUPID-Mo experiment for neutrinoless double-beta decay: performance and prospects, Eur. Phys. J. C 80, 44 (2020a).
  * Deppisch  _et al._ [2012] F. F. Deppisch, M. Hirsch, and H. Päs, Neutrinoless double-beta decay and physics beyond the standard model, J. Phys G: Nucl. Part. Phys. 39, 124007 (2012).
  * Rodejohann [2012] W. Rodejohann, Neutrinoless double-beta decay and neutrino physics, J. Phys G: Nucl. Part. Phys. 39, 124008 (2012).
  * Prézeau  _et al._ [2003] G. Prézeau, M. Ramsey-Musolf, and P. Vogel, Neutrinoless double $\beta$ decay and effective field theory, Phys. Rev. D 68, 034016 (2003).
  * Atre  _et al._ [2009] A. Atre, T. Han, S. Pascoli, and B. Zhang, The search for heavy Majorana neutrinos, Journal of High Energy Physics 2009, 030 (2009).
  * Blennow  _et al._ [2010] M. Blennow, E. Fernandez-Martinez, J. Lopez-Pavon, and J. Menéndez, Neutrinoless double beta decay in seesaw models, Journal of High Energy Physics 2010, 96 (2010).
  * Mitra  _et al._ [2012] M. Mitra, G. Senjanović, and F. Vissani, Neutrinoless double beta decay and heavy sterile neutrinos, Nuclear Physics B 856, 26 (2012).
  * Cirigliano  _et al._ [2018] V. Cirigliano, W. Dekens, J. de Vries, M. L. Graesser, and E. Mereghetti, A neutrinoless double beta decay master formula from effective field theory, Journal of High Energy Physics 2018, 97 (2018).
  * Benato [2015] G. Benato, Effective Majorana Mass and Neutrinoless Double Beta Decay, Eur. Phys. J. C 75, 563 (2015), arXiv:1510.01089 [hep-ph] .
  * Engel and Menéndez [2017] J. Engel and J. Menéndez, Status and Future of Nuclear Matrix Elements for Neutrinoless Double-Beta Decay: A Review, Rept. Prog. Phys. 80, 046301 (2017), arXiv:1610.06548 [nucl-th] .
  * Rath  _et al._ [2013] P. K. Rath, R. Chandra, K. Chaturvedi, P. Lohani, P. K. Raina, and J. G. Hirsch, Neutrinoless $\beta$$\beta$ decay transition matrix elements within mechanisms involving light Majorana neutrinos, classical Majorons, and sterile neutrinos, Phys. Rev. C 88, 064322 (2013), arXiv:1308.0460 [nucl-th] .
  * Šimkovic  _et al._ [2013] F. Šimkovic, V. Rodin, A. Faessler, and P. Vogel, 0$\nu$$\beta$$\beta$ and 2$\nu$$\beta$$\beta$ nuclear matrix elements, quasiparticle random-phase approximation, and isospin symmetry restoration, Phys. Rev. C 87, 045501 (2013), arXiv:1302.1509 [nucl-th] .
  * Vaquero  _et al._ [2013] N. L. Vaquero, T. R. Rodríguez, and J. L. Egido, Shape and pairing fluctuation effects on neutrinoless double beta decay nuclear matrix elements, Phys. Rev. Lett. 111, 142501 (2013).
  * Barea  _et al._ [2015a] J. Barea, J. Kotila, and F. Iachello, 0$\nu\beta\beta$ and 2$\nu\beta\beta$ Nuclear Matrix Elements in the Interacting Boson Model With Isospin Restoration, Phys. Rev. C 91, 034304 (2015a).
  * Agostini  _et al._ [2020] M. Agostini _et al._ (GERDA), Final Results of GERDA on the Search for Neutrinoless Double-$\beta$ Decay, Phys. Rev. Lett. 125, 252502 (2020).
  * Gando  _et al._ [2016] A. Gando _et al._ (KamLAND-Zen), Search for Majorana Neutrinos Near the Inverted Mass Hierarchy Region with KamLAND-Zen, Phys. Rev. Lett. 117, 082503 (2016).
  * Azzolini  _et al._ [2018a] O. Azzolini _et al._ (CUPID-0), First Result on the Neutrinoless Double-$\beta$ Decay of 82Se with CUPID-0, Phys. Rev. Lett. 120, 232502 (2018a).
  * Augier  _et al._ [2022] C. Augier _et al._ , Final results on the $0\nu\beta\beta$ decay half-life limit of 100Mo from the CUPID-Mo experiment, (2022), arXiv:2202.08716 [nucl-ex] .
  * Dolgov and Smirnov [2005] A. Dolgov and A. Smirnov, Possible violation of the spin-statistics relation for neutrinos: Cosmological and astrophysical consequences, Physics Letters B 621, 1 (2005).
  * Barabash  _et al._ [2007] A. S. Barabash, A. D. Dolgov, R. Dvornicky, F. Simkovic, and A. Y. Smirnov, Statistics of neutrinos and the double beta decay, Nucl. Phys. B 783, 90 (2007), arXiv:0704.2944 [hep-ph] .
  * Hubbell and Seltzer [2004] J. H. Hubbell and S. M. Seltzer, X-Ray mass attenuation coefficients, NIST Standard Reference Database 126 (2004).
  * Berger  _et al._ [2017] M. Berger _et al._ , Stopping-Power & Range Tables for Electrons, Protons, and Helium Ions, NIST Standard Reference Database 124 https://dx.doi.org/10.18434/T4NC7P (2017).
  * Barabash [1990] A. S. Barabash, A Possibility for experimentally observing two neutrino double beta decay, JETP Lett. 51, 207 (1990).
  * Barabash  _et al._ [1995] A. S. Barabash _et al._ , Two neutrino double beta decay of 100Mo to the first excited $0^{+}$ state in 100Ru, Phys. Lett. B 345, 408 (1995).
  * Barabash  _et al._ [1999] A. S. Barabash, V. I. Umatov, R. Gurriaran, F. Hubert, and P. Hubert, 2$\nu\beta\beta$ decay of 100Mo to the first 0+ excited state in 100Ru, Phys. Atom. Nucl. 62, 2039 (1999).
  * Arnold  _et al._ [2007] R. Arnold _et al._ (NEMO), Measurement of double beta decay of 100Mo to excited states in the NEMO 3 experiment, Nucl. Phys. A 781, 209 (2007), arXiv:hep-ex/0609058 .
  * Kidd  _et al._ [2009] M. F. Kidd, J. H. Esterline, W. Tornow, A. S. Barabash, and V. I. Umatov, New Results for Double-Beta Decay of 100Mo to Excited Final States of 100Ru Using the TUNL-ITEP Apparatus, Nucl. Phys. A 821, 251 (2009), arXiv:0902.4418 [nucl-ex] .
  * Belli  _et al._ [2010] P. Belli _et al._ , New observation of $2\beta 2\nu$ decay of 100Mo to the $0_{1}^{+}$ level of 100Ru in the ARMONIA experiment, Nuclear Physics A 846, 143 (2010).
  * Arnold  _et al._ [2014] R. Arnold _et al._ , Investigation of double beta decay of 100Mo to excited states of 100Ru, Nuclear Physics A 925, 25 (2014).
  * Belli  _et al._ [2020] P. Belli _et al._ , Double Beta Decay to Excited States of Daughter Nuclei, Universe 6, 239 (2020).
  * Armengaud  _et al._ [2017a] E. Armengaud _et al._ (EDELWEISS), Performance of the EDELWEISS-III experiment for direct dark matter searches, JINST 12 (08), P08010 (2017), arXiv:1706.01070 [physics.ins-det] .
  * Armengaud  _et al._ [2017b] E. Armengaud _et al._ (LUMINEU), Development of 100Mo-containing scintillating bolometers for a high-sensitivity neutrinoless double-beta decay search, Eur. Phys. J. C 77, 785 (2017b).
  * Haller  _et al._ [1984] E. E. Haller, N. P. Palaio, M. Rodder, W. L. Hansen, and E. Kreysa, NTD Germanium: A Novel Material for Low Temperature Bolometers, in  _Neutron Transmutation Doping of Semiconductor Materials_ , edited by R. D. Larrabee (Springer US, Boston, MA, 1984) pp. 21–36.
  * Armengaud  _et al._ [2021] E. Armengaud _et al._ (CUPID), New Limit for Neutrinoless Double-Beta Decay of 100Mo from the CUPID-Mo Experiment, Phys. Rev. Lett. 126, 181802 (2021), arXiv:2011.13243 [nucl-ex] .
  * Adams  _et al._ [2021] D. Q. Adams _et al._ (CUORE), Search for double-beta decay of $\mathrm{{}^{130}Te}$ to the $0^{+}$ states of $\mathrm{{}^{130}Xe}$ with CUORE, Eur. Phys. J. C 81, 567 (2021), arXiv:2101.10702 [nucl-ex] .
  * Alduino  _et al._ [2019] C. Alduino _et al._ (CUORE), Double-beta decay of ${}^{130}\hbox{Te}$ to the first $0^{+}$ excited state of ${}^{130}\hbox{Xe}$ with CUORE-0, Eur. Phys. J. C 79, 795 (2019), arXiv:1811.10363 [nucl-ex] .
  * Azzolini  _et al._ [2018b] O. Azzolini _et al._ (CUPID-0), Search of the neutrino-less double beta decay of 82Se into the excited states of 82Kr with CUPID-0, Eur. Phys. J. C 78, 888 (2018b).
  * Allison  _et al._ [2016] J. Allison _et al._ , Recent developments in Geant4, Nucl. Instrum. Meth. A 835, 186 (2016).
  * Armengaud  _et al._ [2013] E. Armengaud _et al._ (EDELWEISS), Background studies for the EDELWEISS dark matter experiment, Astropart. Phys. 47, 1 (2013), arXiv:1305.3628 [physics.ins-det] .
  * Hamilton [1940] D. R. Hamilton, On Directional Correlation of Successive Quanta, Phys. Rev. 58, 122 (1940).
  * Chauvie  _et al._ [2004] S. Chauvie _et al._ , Geant4 low energy electromagnetic physics, in  _IEEE Symposium Conference Record Nuclear Science 2004._, Vol. 3 (2004) pp. 1881–1885.
  * Simkovic  _et al._ [2001] F. Simkovic, P. Domin, and S. V. Semenov, The Single state dominance hypothesis and the two neutrino double beta decay of 100Mo, J. Phys. G 27, 2233 (2001), arXiv:nucl-th/0006084 .
  * Arnold  _et al._ [2019] R. Arnold _et al._ (NEMO-3), Detailed studies of 100Mo two-neutrino double beta decay in NEMO-3, Eur. Phys. J. C 79, 440 (2019).
  * Armengaud  _et al._ [2020b] E. Armengaud _et al._ (CUPID-Mo), Precise measurement of $2\nu\beta\beta$ decay of 100Mo with the CUPID-Mo detection technology, Eur. Phys. J. C 80, 674 (2020b).
  * Alduino  _et al._ [2016] C. Alduino _et al._ (CUORE), Analysis Techniques for the Evaluation of the Neutrinoless Double-Beta Decay Lifetime in 130Te with CUORE-0, Phys. Rev. C 93, 045503 (2016).
  * Azzolini  _et al._ [2018c] O. Azzolini _et al._ (CUPID-0), Analysis of cryogenic calorimeters with light and heat read-out for double beta decay searches, Eur. Phys. J. C 78, 734 (2018c).
  * Gatti and Manfredi [1986] E. Gatti and P. Manfredi, Processing the signals from solid-state detectors in elementary-particle physics, Riv. Nuovo Cim. 9, 1 (1986).
  * Huang  _et al._ [2021] R. Huang _et al._ (CUPID), Pulse Shape Discrimination in CUPID-Mo using Principal Component Analysis, JINST 16 (03), P03032 (2021), arXiv:2010.04033 [physics.data-an] .
  * Caldwell  _et al._ [2009] A. Caldwell, D. Kollár, and K. Kröninger, BAT – The Bayesian analysis toolkit, Comput. Phys. Commun. 180, 2197 (2009).
  * Weidenspointner  _et al._ [2013] G. Weidenspointner _et al._ , Validation of Compton scattering Monte Carlo simulation models, in  _2013 IEEE Nuclear Science Symposium and Medical Imaging Conference and Workshop on Room-Temperature Semiconductor Detectors_ (2013).
  * Cirrone  _et al._ [2010] G. A. P. Cirrone, G. Cuttone, F. Di Rosa, L. Pandola, F. Romano, and Q. Zhang, Validation of the Geant4 electromagnetic photon cross-sections for elements and compounds, Nucl. Instrum. Meth. A 618, 315 (2010).
  * Kotila and Iachello [2012] J. Kotila and F. Iachello, Phase-space factors for double-$\beta$ decay, Phys. Rev. C 85, 034316 (2012).
  * Coraggio  _et al._ [2022] L. Coraggio _et al._ , Shell-model calculation of 100Mo double-$\beta$ decay, Phys. Rev. C 105, 034312 (2022), arXiv:2203.01013 [nucl-th] .
  * Barea  _et al._ [2015b] J. Barea, J. Kotila, and F. Iachello, $0\nu\beta\beta$ and $2\nu\beta\beta$ nuclear matrix elements in the interacting boson model with isospin restoration, Phys. Rev. C 91, 034304 (2015b).
  * Pirinen and Suhonen [2015] P. Pirinen and J. Suhonen, Systematic approach to $\beta$ and $2\nu\beta\beta$ decays of mass $A=100$–$136$ nuclei, Phys. Rev. C 91, 054309 (2015).
  * Bandac  _et al._ [2020] I. C. Bandac _et al._ (CROSS), The $0\nu 2\beta$-decay CROSS experiment: preliminary results and prospects, JHEP 01, 018 (2020), arXiv:1906.10233 [nucl-ex] .

##  Appendix A Optimization of categories

As explained in Section III, some categories are divided up by their projected
out energies. In particular, for the $2\nu\beta\beta$ decay analysis the
difference in change of spin and Q-values, $Q_{2\beta,0_{1}^{+}}=1904.1$ keV
and $Q_{2\beta,2_{1}^{+}}=2494.9$ keV, leads to an expected difference in the
$E_{\beta\beta}$ spectral shape. This can be used to help reduce correlation
between the decays to these two states.

We use a preliminary background model fit to optimize these choices. We
emphasize that this fit is only used for the optimization of the choice of
categories and not in the Bayesian analysis. We first optimize a choice of
$\beta\beta$ energies to minimize the expected measurement error
($\sqrt{S+B}/S$) for $2\nu\beta\beta$ to $0_{1}^{+}$ e.s. decay. This
optimization is performed on the 540 keV peak, however a consistent result is
also obtained with 591 keV. A separate optimization is performed for the
vertical line where the peak is in $E_{2}$ and a horizontal where the peak is
in $E_{1}$ (see Fig. 4). This leads to accepting only events with $E_{2}>220$
keV in the first case, and $E_{1}<1900$ keV in the second.

Next we optimize the choice of energies, $E_{a}$, such that the beta energy
with a $\gamma$ in $E_{1}$ (horizontal line) is divided into two slices
separated by $E_{a}$:

$\displaystyle E_{2}\in$ $\displaystyle[220\ \mathrm{keV}-E_{a}],\ [E_{a}-500\
\mathrm{keV}].$ (37)

Similarly for the events with the $\gamma$ in $E_{2}$ (vertical line) we
divide into slices separated by the energies $E_{b},E_{c}$:

$\displaystyle E_{1}\in$ $\displaystyle[500\ \mathrm{keV}-E_{b}],\
[E_{b}-E_{\text{c}}],\ [E_{c}-1900\ \mathrm{keV}].$ (38)

This optimization maximizes the limit setting sensitivity for $2\nu\beta\beta$
to $2_{1}^{+}$ e.s. We estimate this as $S/\sqrt{B}$ using the 540 keV peak,
this is a proxy for the correlation between the two decays. These
optimizations lead to $E_{a},E_{b},E_{c}=410,890,1190$ keV respectively.

We use a similar optimization scheme for the 0$\nu\beta\beta$ decay analysis.
In this case the relevant categories are horizontal lines with peaks in
$E_{1}$ (see Fig. 4). We therefore make cuts only on the $E_{2}$ variable.
This is performed separately for the peaks at 1900 keV and $\sim$ $2400$ keV,
and results in the categories shown in the Table 2. Here we optimize by
maximizing the approximate sensitivity to $0_{1}^{+}$ decay for the 1904 keV
peak, and the $2_{1}^{+}$ decay for the 2400 keV peaks.

##  Appendix B Signal model functions

To model the signal shape as described in Section III we use phenomenological
functions which are often used for modeling the shape of signals in cryogenic
calorimeters (for example see [58, 57]). In particular we use a linear
combination of:

$\displaystyle
f_{\text{peak}}(E;\vec{p},\mu,\vec{R},\sigma)=\sum_{i=1}^{3}p_{i}\cdot\mathcal{N}(E,\mu,R_{i}\cdot\sigma(E)),$
(39) $\displaystyle f_{\text{xray}}(E;\mu,E_{X},\sigma)=\mathcal{N}(E,\mu-
E_{X},\sigma(E)),$ (40) $\displaystyle
f_{\text{step}}(E;E_{min},\mu,s,\sigma)=N\cdot\text{Erfc}\Big{(}\frac{\mu-E}{\sqrt{2}\sigma}\Big{)},$
(41) $\displaystyle f_{\text{step-
neg}}=N\cdot\text{Erfc}\Big{(}\frac{E-\mu}{\sqrt{2}\sigma}\Big{)},$ (42)
$\displaystyle f_{\text{bkg}}(E;s;E_{min},E_{max})=N\cdot(1+s\cdot E),$ (43)
$\displaystyle f_{\text{bump}}=\mathcal{N}(E,\mu,\sigma).$ (44)

Here $\mathcal{N}(E,\mu,\sigma)$ is a Gaussian with mean $\mu$ and standard
deviation $\sigma$, $N$ is a normalization coefficient and Erfc is the
complementary error function. The first function models a peak, either
$\gamma$ or $\beta\beta$, the second models events where a Mo X-ray escapes
the crystal. The third and fourth functions (step/step-neg.) are step
functions to account for Compton scattering events. The step function accounts
for the Compton scatter of a single $\gamma$ leading to a partial energy
deposition and $E<E_{\text{photopeak}}$ and the step negative accounts for a
combination of photo-absorption of one $\gamma$ and Compton scatter of the
other and $E>E_{\text{photopeak}}$. We include a linear background and a
single Gaussian to model some features in the data where a diagonal line
crosses the projected out energy leading to a very small but broader peak
(bump). For example see the category of $0\nu\beta\beta$ decay to $2_{1}^{+}$
e.s. with $E_{1}=1904$ keV and $E_{2}\in[500,650]$ keV, here a diagonal line
crosses the box at around 591 keV. This line is caused by a 2494 keV
$\beta\beta$ in one crystal and a 540 keV $\gamma$ shared between this crystal
and another.

##  Appendix C Containment efficiency uncertainty

As mentioned in Section VI.3, several sets of MC simulations were used to
estimate the uncertainty in the containment efficiency. For each peak,
category and systematic test the percentage change in the containment
efficiency was extracted. These values are shown in Fig. 13.

Figure 13: Plots showing the systematic uncertainty on containment efficiency
for each peak, category, decay mode, and systematic effect. The z-axis scale
is the percentage effect. Top left this is $2\nu\beta\beta$ decay to
$0_{1}^{+}$ state and on the right for decay to $2_{1}^{+}$ state. Bottom left
is $0\nu\beta\beta$ decay to $0_{1}^{+}$ state, while bottom right is
$0\nu\beta\beta$ decay to $2_{1}^{+}$ state.

##  Appendix D Fits to each category

We show the best fit reproduction of each category for both $2\nu\beta\beta$
and $0\nu\beta\beta$ decay analysis in Figs. 14 and 15. Note that the spectra
are binned for visualization purposes and we exclude plots with 0 counts in
the experimental data. For the $2\nu\beta\beta$ decay fits we use the
$H_{1}=(0_{1}^{+}+B)$ model.

Figure 14: Fits of the seven $2\nu\beta\beta$ decay regions of interest, the
experimental data is shown in grey while the best fit (global mode) is shown
in blue. We see clearly a signal in each of the regions of interest.

Figure 15: Fits of the $0\nu\beta\beta$ decay regions of interest, the
experimental data is shown in black while the best fit (global mode) is shown
in blue. We only show categories with non-zero counts in the fit region.