1 Introduction

Spontaneous nuclear double beta decay is a second order weak interaction process that was theoretically considered for the first time by Goeppert-Mayer [1]. It can occur in some even-even nuclei when two bound neutrons simultaneously undergo beta decay and are transformed into two bound protons emitting two electrons and two (anti)neutrinos. Two-neutrino double beta decay, \(2\nu \beta \beta \), is one of the rarest directly observed radioactive processes with half-lives ranging from \(7\times 10^{18}\) to \(2\times 10^{21}\) years [2, 3].

The decay rate of \(2\nu \beta \beta \) decay can be expressed as

$$\begin{aligned} 1/T^{2\nu }_{1/2} =g_{A}^4 G^{2\nu } |M^{2\nu }|^2, \end{aligned}$$
(1)

where \(g_{A}\) is the axial-vector coupling constant, \(G^{2\nu }\) is a phase space factor, and \(M^{2\nu }\) is a nuclear matrix element (NME). Measurement of the \(2\nu \beta \beta \) half-life gives direct access to the value of the NME for this process and therefore provides experimental input into nuclear models that are used to evaluate NMEs. Moreover, \(2\nu \beta \beta \) may provide answers to the question of \(g_{A}\) quenching in nuclear matter that is currently being actively discussed [4,5,6,7,8]. Detailed studies of \(2\nu \beta \beta \) may therefore be useful to improve NME calculations for the neutrinoless mode of double beta decay, \(0\nu \beta \beta \), the process which violates total lepton number and is one of the most sensitive probes of physics beyond the Standard Model. A recent review of the \(0\nu \beta \beta \) NME calculation methods, challenges and prospects can be found in [9].

Previous measurements have shown that the \(^{100}\)Mo \(2\nu \beta \beta \) half-life is shorter compared to other \(\beta \beta \) isotopes [10,11,12,13,14,15,16,17], and it is therefore a promising nucleus for precise studies of the process. Here we present the most accurate to date study of \(^{100}\)Mo \(2\nu \beta \beta \) decay including single electron energy and angular distributions of the electrons emitted in the decay with an unprecedented statistics of \(5\times 10^5\) events. The impact of the single electron energy spectra on nuclear models that are used to calculate the NME is also presented.

Searches for most commonly discussed \(0\nu \beta \beta \) mechanisms (exchange of a light Majorana neutrino, right-handed currents, super-symmetry) with NEMO-3 have been reported earlier in [18, 19]. In this paper we present results obtained for \(^{100}\)Mo \(0\nu \beta \beta \) decay accompanied by the emission of Majoron bosons with spectral indices \(n\ge 2\), as well as constraints on contributions from bosonic neutrinos and from Lorentz invariance violation to \(2\nu \beta \beta \) spectra of \(^{100}\)Mo.

2 The NEMO-3 detector

The NEMO-3 detector, its calibration and performance are described in detail in [20] and more recently in [19]. A combination of tracking and calorimetric approaches allows for a full reconstruction of \(\beta \beta \) event topology. A tracking chamber is used to reconstruct electron tracks, their origin and end points. The electron energies and arrival times are measured with a plastic scintillator calorimeter. The cylindrical detector measuring 3 m in height and 5 m in diameter is made up of 20 wedge-shaped sectors of identical size. Each sector hosts 7 thin foil strips containing a \(\beta \beta \) isotope. The source foils are positioned in the middle of the tracking detector at a radius of 1 m and have a height of 2.48 m.

The tracking detector is based on a wire chamber made of 6180 open drift cells operating in Geiger mode with helium as the main working gas with the addition of ethanol (4%), argon (1%) and water vapour (0.15%). The wire cells are strung vertically parallel to the source foils and have average transverse and longitudinal resolutions of 0.5 mm and 0.8 cm (\(\sigma \)) respectively. The tracking volume is surrounded by a segmented calorimeter composed of 1940 optical modules made of 10 cm thick polystyrene scintillator blocks coupled to low radioactivity photomultiplier tubes (PMT). The energy resolution of optical modules for 1 MeV electrons ranges from 5.8 to 7.2% and the time resolution is 250 ps (\(\sigma \)). The detector was calibrated by deploying \(^{207}\)Bi, \(^{90}\)Sr and \(^{232}\)U sources during the course of data collection. The stability of the PMT gains was monitored by a dedicated light injection system that was run every 12 hours.

The NEMO-3 detector is supplied with a solenoid which generates a 25 G magnetic field parallel to the tracking detector wires and provides charge identification by track curvature. The detector is surrounded by passive shielding consisting of a 19 cm thick iron plates to suppress the external gamma ray flux, and of borated water, paraffin and wood to moderate and absorb environmental neutrons.

One of the unique advantages of the NEMO-3 technology is the ability to unambiguously identify electrons, positrons, gamma- and delayed alpha-particles. This approach leads to a strong suppression of backgrounds by eliminating events that do not exhibit a \(\beta \beta \) topology. In addition, it allows for an efficient background evaluation by selecting event topologies corresponding to specific background channels. An electron is identified by a reconstructed prompt track in the drift chamber matching to a calorimeter deposit. Extrapolating the track to the foil plane defines the event vertex in the source. The track extrapolation to the calorimeter identifies the impact point of the electron track with the corresponding optical module and is used to correct the reconstructed energy of the electron deposited in the scintillator. The track curvature in the magnetic field is used to distinguish electrons from positrons. A \(\gamma \)-ray is identified as an energy deposit in the calorimeter without an associated track in the drift chamber. An \(\alpha \)-particle is identified by a short straight track delayed with respect to the prompt electron in order to tag \(^{214}\)Bi \(\rightarrow \) \(^{214}\)Po delayed coincidences.

The NEMO-3 detector took data at the Modane Underground Laboratory (LSM) in the Frejus tunnel at a depth of 4800 m w.e. enabling the cosmic muon flux suppression by a factor of > 10\(^{6}\). The detector hosted source foils of 7 different \(\beta \beta \) isotopes. The two isotopes with the largest mass were \(^{100}\)Mo (6.914 kg) [19] and \(^{82}\)Se (0.932 kg) [21] with smaller amounts of \(^{48}\)Ca, \(^{96}\)Zr, \(^{116}\)Cd, \(^{130}\)Te and \(^{150}\)Nd [22,23,24,25,26].

Two types of purified molybdenum foils were installed in NEMO-3, metallic and composite. Both foil types were enriched in \(^{100}\)Mo with the isotopic enrichment factor ranging from \(95.14 \pm 0.05\) to \(98.95 \pm 0.05\%\). The average enrichment factor was 97.7% for metallic foils and 96.5% for composite foils. The metallic foils contained \(2479\pm 5\) g of \(^{100}\)Mo. The mean metallic foil density is 58 mg/\(\hbox {cm}^2\) with a total foil surface of 43,924 \(\hbox {cm}^2\). The composite foils contained \(4435\pm 13\) g of \(^{100}\)Mo. They were produced by mixing a fine molybdenum powder with polyvinyl alcohol (PVA) glue and deposited between Mylar foils of 19 \(\upmu \)m thickness. The average surface density of the composite foils is 66 mg/\(\hbox {cm}^2\) and the total foil surface area is 84,410 \(\hbox {cm}^2\).

Monte Carlo (MC) simulations are performed with a GEANT-3 based [27] program using the DECAY0 [28] event generator. The time-dependent status and performance of the detector are taken into account in modelling the detector response.

The data presented here were collected between February 2003 and October 2010 with a live time of 4.96 years and a total exposure of 34.3 kg year of \(^{100}\)Mo. This is the same exposure as that used for \(0\nu \beta \beta \) results published earlier [19].

3 Background model

Trace quantities of naturally-occurring radioactive isotopes can occasionally produce two-electron events and thus can mimic \(\beta \beta \)-decay events. The largest contributions come from isotopes that are progenies of \(^{238}\)U (\(^{234m}\)Pa, \(^{214}\)Pb, \(^{214}\)Bi, \(^{210}\)Bi) and of \(^{232}\)Th (\(^{228}\)Ac, \(^{212}\)Bi, \(^{208}\)Tl), as well as \(^{40}\)K.

The background is categorised as internal if it originates from radioactive decays inside the \(\beta \beta \) source foils, see Fig. 1a. Two electrons can be produced via \(\beta \)-decay followed by a Møller scattering, \(\beta \)-decay to an excited state with the subsequent internal conversion or due to Compton scattering of the de-excitation photon.

Fig. 1
figure 1

Mechanisms of internal (a) and external (b) background production in the source foil

Fig. 2
figure 2

Single electron events energy spectra for metallic and composite molybdenum. The error bars correspond to statistical uncertainty only

Decays inside the tracking detector volume form a separate background category. The main source of this background is radon, \(^{222}\)Rn. The decay of radon progenies near the source foil can produce signal-like events in an analogous manner to internal background decays.

The last background category is due to the external \(\gamma \)-ray flux produced by decay of radioactive isotopes in detector components, the surrounding area and due to neutron interactions in the shield and material of the detector. The PMT glass is the main source of these \(\gamma \)-rays. They can produce two-electron events due to \(e^+ e^-\) pair creation in the source foil and subsequent charge misidentification, double Compton scattering or Compton scattering followed by Møller scattering, see Fig. 1b.

A detailed discussion of the NEMO-3 background model is presented in [29] and results of screening measurements can be found in [19, 20, 29]. Here we follow the same background model as that presented for the \(^{100}\)Mo \(0\nu \beta \beta \) analysis [19]. However, radioactive isotopes contributing to the low energy region of the \(^{100}\)Mo \(2\nu \beta \beta \) spectrum were not relevant for the \(0\nu \beta \beta \) analysis in [19] and are therefore discussed in more detail below. The background in question comes from traces of \(\beta \)-decaying isotopes \(^{210}\)Bi, \(^{40}\)K and \(^{234m}\)Pa in \(^{100}\)Mo foils. In addition, \(^{100}\)Mo \(2\nu \beta \beta \) decay to the \(0^+_1\) excited state of \(^{100}\)Ru is also taken into account as a source of internal background. The experimental half-life value of \(T_{1/2} = 6.7^{+0.5}_{-0.4} \times 10^{20}\) years [3] is used to evaluate this contribution.

The activities of \(\beta \)-emitters in \(^{100}\)Mo foils are determined from the fit to the electron energy distribution for a single electron event sample, which is shown in Fig. 2 separately for metallic and composite foils. To disentangle the \(^{210}\)Bi contribution from the source foils and the surface of the tracker wires the activity measured in [29] is used for the latter. Figure 2 shows the sum of both contributions. Secular equilibrium is assumed between \(^{214}\)Pb and \(^{214}\)Bi. The same is done between \(^{228}\)Ac, \(^{212}\)Bi and \(^{208}\)Tl, where the branching ratio of 35.94% is taken into account. There is sufficiently good agreement between data and MC for the single electron energy spectrum. The observed deviations of MC from data are within 6% and are not significant when the systematic uncertainty on the external background is taken into account.

The results of the internal \(^{100}\)Mo foil contamination measurements carried out with the NEMO-3 detector are shown in Table 1.

Table 1 \(^{100}\)Mo source foil contamination activities measured with the NEMO-3 detector. Activities of \(^{214}\)Bi and \(^{208}\)Tl are from [19]

4 Two-neutrino double beta decay of \(^{100}\)Mo

Candidate \(\beta \beta \) events are selected by requiring two reconstructed electron tracks, each associated with an energy deposited in an individual optical module. The energy deposited by the electron in a single optical module should be greater than 300 keV. Each PMT must be flagged as stable according to the light injection survey [19]. The tracks must both originate from the \(^{100}\)Mo source foil, and their points of intersection with the plane of the source foil must be within 4 cm transverse to and 8 cm along the direction of the tracker wires, in order to ensure that the two tracks are associated to a common event vertex. The track curvatures must be consistent with electrons moving outwards from the source foil. The timing and the path length of the electrons must be consistent with the hypothesis of simultaneous emission of two electrons from a common vertex in the \(^{100}\)Mo source foil [19]. There should be no \(\gamma \)-ray hits and \(\alpha \)-particle tracks in the event.

After the above event selection there are 501,534 \(^{100}\)Mo two-electron candidate events, with 193,699 coming from the metallic foils and 307,835 from the composite foils. Table 2 shows the number of expected background and candidate signal events in \(^{100}\)Mo foils. The number of \(2\nu \beta \beta \) events is obtained from a binned log-likelihood fit to the two-electron energy sum distribution under the single state dominance (SSD) nuclear model, as detailed below. The average signal-to-background ratio is \(\mathrm{S/B}=79\), with \(\mathrm{S/B}=63\) for the metallic foils and \(\mathrm{S/B}=94\) for the composite foils. The detector acceptance and selection efficiency for \(2\nu \beta \beta \) \(^{100}\)Mo events calculated using MC simulations is \(\epsilon = (2.356\pm 0.002)\%\), with \(\epsilon _{met} = (2.472\pm 0.003)\%\) and \(\epsilon _{com} = (2.292\pm 0.002)\%\) for the metallic and composite molybdenum foils respectively. Using the above values gives the \(^{100}\)Mo \(2\nu \beta \beta \)-decay half-life of \(T_{1/2} = (6.65 \pm 0.02 ) \times 10^{18}\) year for the metallic foils and \(T_{1/2} = (6.91 \pm 0.01 ) \times 10^{18}\) year for the composite foils. The difference between the two sample measurements may be explained by inaccuracy of the thin foil modelling and is taken into account in estimation of the systematic uncertainty in Sect. 4.2. We consider the mean value over the two data samples as the more reliable half-life estimation

$$\begin{aligned} T_{1/2} = (6.81 \pm 0.01 ) \times 10^{18}~\text{ year }. \end{aligned}$$
(2)
Table 2 Expected number of background events in the two-electron channel and the number of \(^{100}\)Mo \(2\nu \beta \beta \) candidate events in molybdenum foils

The two-electron energy sum spectra and the distributions of cosine of the angle between two electrons emitted from \(^{100}\)Mo foil are shown in Fig. 3, separately for the metallic and composite foils as well as for the total \(^{100}\)Mo sample.

Fig. 3
figure 3

Distributions of two-electron summed kinetic energy and the opening angle between two electron tracks in \(^{100}\)Mo foils after an exposure of 34.3 kg year. Data are compared to the MC prediction of the SSD model (see text), where the resulting event numbers are taken from a binned log-likelihood fit

The electron energy measured in the calorimeter is smaller than the energy at the point of origin due to energy losses in the foil and in the drift chamber. For instance in the case of \(^{100}\)Mo \(2\nu \beta \beta \) decay the mean electron track length from the source foil to the calorimeter is 75 cm and the mean energy loss of electrons in the drift chamber is 43 keV. The single and summed electron energy distributions are presented for the measured values of the electron kinetic energy \(E_e\) and sum of the measured electron kinetic energies \(E_{SUM}\), respectively, i.e., without correction for the energy loss.

The angular distribution is corrected with the well-measured distribution of the opening angle between two electrons emitted in \(^{207}\)Bi decay. The MC distribution of the cosine of the angle between two electron tracks has been reweighted based on data collected in the regular energy calibration runs performed with \(^{207}\)Bi sources. The correction is biggest for small opening angles, and is at the level of 4% on average.

4.1 Role of intermediate nuclear states in \(^{100}\)Mo \(2\nu \beta \beta \) transition

The nuclear \(\beta \beta \) decay (A,Z) \(\rightarrow \) (A,Z+2) is realized via two subsequent virtual \(\beta \) transitions through the complete set of states of intermediate nucleus (A,Z+1). In the case of \(^{100}\)Mo \(2\nu \beta \beta \) transition between the ground states of the parent (\(^{100}\)Mo) and daughter (\(^{100}\)Ru) nuclei with spin-parity \(0^+\) the process is governed by two Gamow-Teller transitions through \(1^+\) states of \(^{100}\)Tc. Nuclear theory does not predict a priori whether there is a dominance of transition through the \(1^+\) ground state (SSD hypothesis [30,31,32]) or through higher lying excited states, namely from the region of the Gamow-Teller resonance (HSD hypothesis). The SSD versus HSD analysis is feasible as the ground state of \(^{100}\)Tc has spin-parity \(J^P=1^+\) and is lying close to the ground state of \(^{100}\)Mo.

The evidence in favour of SSD in \(^{100}\)Mo \(2\nu \beta \beta \) decay was already observed at the beginning of NEMO-3 data analysis [33]. Further hints for the SSD model in the \(^{100}\)Mo \(2\nu \beta \beta \) decay were obtained in charge-exchange experiments by observing a strong Gamow-Teller transition to the \(1^+\) ground state of \(^{100}\)Tc in the \(^{100}\)Mo(\(^{3}\)He,t)\(^{100}\)Tc reaction [34]. It was estimated that this transition could contribute as much as 80% to the total value of the \(^{100}\)Mo \(2\nu \beta \beta \) matrix element.

It was shown in [31, 32] that SSD and HSD models can be directly distinguished by making high precision kinematics measurements of \(2\nu \beta \beta \) decay products. The distribution of the individual electron energies was shown to have the most discriminating power, especially in the low energy part of the spectrum. Figure 4 shows the individual electron energy spectra for three nuclear models, with SSD-3 being a modification of the SSD model where a finer structure of intermediate states is accounted for [35, 36].

Figure 5 shows the energy sum and angular distribution of the final state electrons where the data are fitted with the HSD model. The tension between the data and the model is evident already from these distributions with \(\chi ^2/\mathrm{ndf}=4.57\) (\(p \,\mathrm{value}=5.3\times 10^{-12}\)) and \(\chi ^2/\mathrm{ndf}=1.98\) (\(p \,\mathrm{value}=0.007\)) for the energy sum and angular distributions respectively. However, the strongest evidence comes from the single electron energy distributions shown in Fig. 6 for the three models, HSD, SSD and SSD-3, fitted to the data. It is clear from the distributions and \(\chi ^2\) values that the HSD model can be ruled out with high confidence while SSD and SSD-3 provide a fairly good description of the data.

The difference between SSD and SSD-3 in describing the data is maximised with a cut on the electron energy sum of \(E_{SUM} > 1.4\) MeV as shown in Fig. 7, which also increases the signal-to-background ratio. There is a slight preference of the SSD-3 model over SSD in this case, contrary to the results obtained without this cut demonstrated at Fig. 6. Due to systematic effects connected to the energy reconstruction and electron energy loss simulations discussed below these two models cannot be discriminated against each other. The SSD is chosen as the baseline model and is used to estimate the \(^{100}\)Mo \(2\nu \beta \beta \) half-life (see Sect. 4 and Fig. 3). We note that differences in the low energy part of the single electron spectra (Fig. 4) affect the selection efficiency of \(^{100}\)Mo \(2\nu \beta \beta \) events. Consequently, the measured half-life for the SSD model is 14% shorter than the analogous result for the HSD model. The SSD-3 model would give a 1.8% shorter half-life than that of the SSD model.

Fig. 4
figure 4

Theoretical distributions of the individual electron kinetic energy for three models of \(^{100}\)Mo \(2\nu \beta \beta \) decay: HSD, SSD and SSD-3

Fig. 5
figure 5

Two-electron events. Energy sum and cosine of the angle between the two electrons for HSD model

Fig. 6
figure 6

Distribution of individual electron kinetic energy in the \(\beta \beta \) channel from \(^{100}\)Mo foils compared with MC spectra under the HSD, SSD and SSD-3 nuclear models. The HSD hypothesis is excluded (\(\chi ^{2}/\text {ndf} = 1159/27\)) while the data are consistent with the SSD and SSD-3 models (\(\chi ^{2}/\text {ndf} = 41.5/27\) and \(\chi ^{2}/\text {ndf} = 49.7/27\) respectively)

Fig. 7
figure 7

Distribution of individual electron kinetic energy in the \(\beta \beta \) channel from \(^{100}\)Mo foils with the cut on the summed electron energy \(E_{SUM} > 1.4\) MeV to maximise the signal-to-background ratio. The data are compared with MC spectra under the HSD, SSD and SSD-3 nuclear models. The HSD hypothesis is excluded (\(\chi ^{2}/\text {ndf} = 1508/27\)) while the data are consistent with the SSD and SSD-3 models (\(\chi ^{2}/\text {ndf} = 39/27\) and \(\chi ^{2}/\text {ndf} = 30.6/27\) respectively)

4.2 Systematic uncertainties on \(^{100}\)Mo \(2\nu \beta \beta \) half-life

Apart from the statistical uncertainties on the fitted number of signal events, the measurement of the \(2\nu \beta \beta \) decay half-life is subject to a number of systematic uncertainties.

The uncertainty on the reconstruction and selection efficiency including the detector acceptance effects is evaluated by carrying out dedicated calibrations with \(^{207}\)Bi sources whose activities were known with a 5% uncertainty. Consequently, the systematic error on the signal efficiency is taken to be 5%.

Limited precision of MC simulation program in modelling of multiple scattering processes and electron energy losses in molybdenum \(\beta \beta \) source foils also contribute to the total systematic error. Corresponding uncertainty is evaluated as the difference between the mean half-life value and the values obtained with metallic (\(-2.3\)%) and composite (\(+1.5\)%) foils.

The 1.8% half-life value difference between the SSD and SSD-3 nuclear models is taken as a systematic error due to the \(^{100}\)Mo \(2\nu \beta \beta \) decay model.

The uncertainty on the energy scale translates into an error on the half-life measurement of 0.6%.

The \(^{100}\)Mo mass uncertainty gives directly the corresponding uncertainty of the half-life value and is estimated to be 0.2%.

The error on the activities of external backgrounds, radon and the foil contamination with \(^{214}\)Bi and \(^{208}\)Tl is 10% as shown in [19]. The uncertainty on the backgrounds from \(^{40}\)K in the source foils as well as from \(^{210}\)Bi is estimated to be 4%. The observed discrepancy in the \(^{234m}\)Pa decay scheme reported in [37] and [38] lead to a 30% normalisation uncertainty on the activity from this isotope. The 7.5% error on the rate of the \(^{100}\)Mo \(2\nu \beta \beta \) decay to the excited states [3] is also taken into account. Overall, due to a high signal-to-background ratio the uncertainty on all background contributions produces only a 0.2% systematic uncertainty on the \(^{100}\)Mo \(2\nu \beta \beta \) half-life determination.

The systematic uncertainties on the measured \(2\nu \beta \beta \) \(^{100}\)Mo half-life are summarised in Table 3. The individual sources of the systematic error are assumed to be uncorrelated and the total uncertainty is obtained to be [\(+5.6,-5.8\)]%.

The final value of the half-life for the \(2\nu \beta \beta \) decay of \(^{100}\)Mo under the SSD model is:

$$\begin{aligned} T_{1/2} = \left[ 6.81 \pm 0.01\,\left( \text{ stat }\right) ^{+0.38}_{-0.40}\,\left( \text{ syst }\right) \right] \times 10^{18}~\text{ year }. \end{aligned}$$
(3)

This value is in good agreement with the world average value of \((7.1 \pm 0.4) \times 10^{18}\) year [3] and with a recent result obtained using low-temperature scintillating bolometers (\(\hbox {Li}_{{2}}^{100}\hbox {MoO}_4\)), \([6.90 \pm 0.15(\text{ stat }) \pm 0.37(\text{ syst })]\times 10^{18}\) year [17].

5 Search for new physics with continuous \(^{100}\)Mo \(\beta \beta \) energy spectra

Deviations in the shape of the \(2\nu \beta \beta \) energy spectra can provide hints of new physics. Below we report on results of searches for physics beyond the Standard Model that can modify the two-electron energy sum distribution of the \(^{100}\)Mo \(2\nu \beta \beta \) decay due to emission of Majoron bosons, the existence of a bosonic component in the neutrino states and possible Lorentz invariance violation.

The shape of the two-electron energy sum distribution in various types of decays is characterized by the spectral index n [39], being determined by the phase space \(G \sim (Q_{\beta \beta }-T)^n\), where \(Q_{\beta \beta }\) is the full energy released in the decay minus two electron masses and T is the sum of kinetic energies of two emitted electrons. The ordinary \(2\nu \beta \beta \) decay has a spectral index of \(n=5\). Any modification from this functional form can be an indication of new physics.

A number of grand unification theories predict the existence of a massless or light boson which couples to the neutrino. Neutrinoless \(\beta \beta \) decay can proceed with the emission of one or two Majoron bosons resulting in a continuous energy sum spectrum with spectral index \(n \ne 5\). The decay accompanied by a single Majoron emission has \(n=1,2\) and 3, while models with two Majoron emissions predict \(n=3\) and 7  (see [40] and references therein). The results for the neutrinoless \(\beta \beta \) decay with the emission of a Majoron corresponding to the spectral index \(n=1\) have already been published in [18, 19]. The Majoron-accompanied \(0\nu \beta \beta \) decay modes with spectral indices \(n=2,3\) and 7 are considered here.

It was noted in [41] that violation of the Pauli exclusion principle resulting in a bosonic component in the neutrino states can be tested by looking at the shape of the energy and angular distributions of the electrons emitted in \(\beta \beta \) decay. For the two-electron energy sum distribution the corresponding index would be \(n=6\).

Lorentz invariance is a fundamental symmetry. However, new physics at very high energies close to the Planck scale can manifest itself in small effects at low energies, including Lorentz invariance violation. Consequently, searches for non-Lorentz invariant effects have attracted active theoretical and experimental effort [42,43,44,45]. The possibility to test Lorentz invariance with \(\beta \beta \) decay was discussed in [46, 47]. In case of \(2\nu \beta \beta \) decay the Lorentz invariance violation may be manifested as a modification of the conventional electron sum spectrum due to an additional contribution of the Lorentz-violating perturbation with a spectral shape of \(n=4\).

Table 3 Summary of systematic uncertainties on the measured \(2\nu \beta \beta \) \(^{100}\)Mo half-life
Fig. 8
figure 8

Spectrum of the of kinetic energy sum of two electrons for the standard \(^{100}\)Mo \(2\nu \beta \beta \) decay (spectral index \(n=5\)) compared to the spectra for neutrinoless \(\beta \beta \) decay with the emission of one or two Majorons \(0\nu Mn\) (\(n = 2, 3, 7\)); shape of the perturbation to the standard \(2\nu \beta \beta \) decay due to Lorentz invariance violation \(2\nu \)-LIV (\(n=4\)) and spectrum for \(2\nu \beta \beta \) decay with bosonic neutrino \(2\nu \)-Boson (\(n=6\))

The theoretical distributions of the two-electron energy sum for different modes of \(^{100}\)Mo \(\beta \beta \) decay discussed above are shown in Fig. 8. The difference in the shape of the distributions due to different spectral indices n is used to evaluate possible contributions from physics beyond the Standard Model. No significant deviations from the expected \(^{100}\)Mo \(2\nu \beta \beta \) spectral shape (\(n=5\)) have been observed and therefore limits on new physics parameters have been set using the full energy sum spectrum of the full \(^{100}\)Mo data set. The contributions of the \(\beta \beta \) decay modes with spectral indices \(n=2,3,6,7\) are constrained with a modified frequentist \(CL_s\) method [48, 49] using a profile likelihood fitting technique (COLLIE software package [50]). A profile likelihood scan is used for the distribution with the spectral index \(n=4\) in order to explore possibility of negative as well as positive Lorentz-violating perturbation.

The systematic uncertainties on background contributions discussed in Sect. 4.2, the 5% uncertainty on the detector acceptance and selection efficiency for signal, a possible distortion in the shape of the two-electron energy sum spectrum due to the energy calibration accuracy, as well as a 5% error on the modelling of the energy loss of electrons are taken into account in limit setting without imposing a constraint on the normalization of standard \(2\nu \beta \beta \) contribution.

The limits on the half-lives for different \(0\nu \beta \beta \) modes with Majoron(s) emission, and for the bosonic neutrino admixture obtained with the \(CL_s\) method are given in Table 4.

Table 4 Lower bounds on half-lives (\(\times 10^{21}\) year) at \(90\%\) C.L. from \(0\nu \beta \beta \) searches with Majoron emission (spectral indices \(n=2,3,7\)), and searches for the bosonic neutrino admixture. The ranges in the expected half-life limits are from the \(\pm 1\sigma \) range of the systematic uncertainties on the background model, signal efficiency and distortions in the shape of the energy spectrum

The half-life limits on the Majoron \(0\nu \beta \beta \) modes are translated into the upper limits on the lepton number violating parameter \(\langle g_{ee}\rangle \), which is proportional to the coupling between the neutrino and the Majoron boson, using the relation,

$$\begin{aligned} 1/T_{1/2} = |\langle g_{ee}\rangle |^m G |M|^2, \end{aligned}$$
(4)

where G is the phase space (which includes the axial-vector coupling constant \(g_A\)), M is the nuclear matrix element, and \(m=2(4)\) is the mode with the emission of one (two) Majoron particle(s). The M and G values are taken from [51]. For the single Majoron emission and \(n=3\), M and G are taken from [52]. There are no NME and phase space calculations available for \(n=2\).

The upper limits on the Majoron-neutrino coupling constant \(\langle g_{ee}\rangle \) are shown in Table 5. One can see that the NEMO-3 results presented here are the current best limits for \(n=3\) and the single Majoron emission mode and are comparable with the world’s best results from the EXO-200 [53] and GERDA [54] experiments for the other two modes.

Table 5 Upper limits on the Majoron-neutrino coupling constant \(\langle g_{ee}\rangle \) from NEMO-3 (\(^{100}\)Mo, this work) and EXO-200 (\(^{136}\)Xe) [53] and GERDA (\(^{76}\)Ge) [54] experiments. All limits are at 90% C.L. The ranges are due to uncertainties in NME calculations

The contribution of bosonic neutrinos to the \(2\nu \beta \beta \)-decay rate can be parametrised as [41]:

$$\begin{aligned} W_{tot} = \cos ^4 \chi W_f + \sin ^4 \chi W_b, \end{aligned}$$
(5)

where \(W_f\) and \(W_b\) are the weights in the neutrino wave-function expression corresponding to the two fermionic and two bosonic antineutrino emission respectively. The purely fermionic, \(T_{1/2}^{f}\), and purely bosonic, \(T_{1/2}^{b}\), half-lives are calculated under the SSD model to be [41] :

$$\begin{aligned} T_{1/2}^f(0^+g.s. )= & {} 6.8\times 10^{18}~\text{ year },\nonumber \\ T_{1/2}^b(0^+g.s.)= & {} 8.9\times 10^{19}~\text{ year }. \end{aligned}$$
(6)

Using the NEMO-3 half-life limit of \(T_{1/2}^b(0^+g.s.) > 1.2\times 10^{21}\) year (Table 4) an upper limit on the bosonic neutrino contribution to the \(^{100}\)Mo \(2\nu \beta \beta \) decay to the ground state can be evaluated as:

$$\begin{aligned} \sin ^2 \chi < 0.27 \,(90\%\,\text{ C.L. }). \end{aligned}$$
(7)

Although this limit is stronger than the bound obtained earlier in [41], the \(2\nu \beta \beta \) transition of \(^{100}\)Mo to the ground state is not very sensitive to bosonic neutrino searches due to a small value of the expected bosonic-to-fermionic decay branching ratio \(r_0 (0^+g.s. ) = 0.076\). The \(^{100}\)Mo \(2\nu \beta \beta \) decay to the first excited \(2^+_1\) state has a branching ratio of \(r_0 (2^+_1 ) = 7.1\) [41] and is therefore potentially more promising despite a lower overall decay rate. The current best experimental limit for this process is \(T_{1/2}(2^+_1) > 2.5\times 10^{21}\) year [55]. This bound is still an order of magnitude lower than the theoretically expected half-life value of \(T_{1/2}^b(2^+_1) = 2.4\times 10^{22}\) year for purely bosonic neutrino, and two orders of magnitude lower than the corresponding expected value for purely fermionic neutrino, \(T_{1/2}^f(2^+_1) = 1.7\times 10^{23}\) year [41].

The Standard Model Extension (SME) provides a general framework for Lorentz invariance violation (LIV) [42]. In this model, the size of the Lorentz symmetry breakdown is controlled by SME coefficients that describe the coupling between standard model particles and background fields. Experimental limits have been set on hundreds of these SME coefficients from constraints in the matter, photon, neutrino and gravity sectors [42]. The first search for LIV in \(2\nu \beta \beta \) decay was carried out in [56]. The two-electron energy sum spectrum of \(^{136}\)Xe was used to set a limit on the parameter \(\mathring{a}^{(3)}_{of}\), which is related to a time-like component of this LIV operator. The value of this parameter was constrained to be \(-2.65\times 10^{-5}\) GeV \(< \mathring{a}^{(3)}_{of} < 7.6\times 10^{-6}\) GeV by looking at deviations from the predicted energy spectrum of \(^{136}\)Xe \(2\nu \beta \beta \) decay [56].

In this work we adopt the same method, using the phase space calculations from [57], and perform a profile likelihood scan over positive and negative contributions of LIV to two-electron events by altering the \(^{100}\)Mo \(2\nu \beta \beta \) energy sum spectrum with positive and negative values of \(\mathring{a}^{(3)}_{of}\). The result of this scan is shown in Fig. 9.

Fig. 9
figure 9

Profile likelihood scan over observed two-electron LIV counts in \(^{100}\)Mo \(2\nu \beta \beta \) energy sum spectrum. The 90% CL exclusion limit is shown with the dashed line

The minimum of the profile log-likelihood function corresponds to \(-135\) counts and is not statistically significant even at 1\(\sigma \) level. The 90% CL exclusion limit is shown in Fig. 9 with the dashed line and gives \(-1798\) and 1527 events for negative and positive contributions to the deviation from the \(^{100}\)Mo \(2\nu \beta \beta \) energy sum spectrum respectively. The corresponding constraint on \(\mathring{a}^{(3)}_{of}\) is calculated using equations (2)–(6) in [56]. The result for \(^{100}\)Mo obtained with a full set of NEMO-3 data is

$$\begin{aligned} -4.2\times 10^{-7}~\text{ GeV }< \mathring{a}^{(3)}_{of} < 3.5\times 10^{-7}~\text{ GeV }\,(90\%\,\text{ C.L. }).\nonumber \\ \end{aligned}$$
(8)

A summary of the best available constraints on LIV and CPT violation parameters can be found in compilation [42].

6 Summary

The results of the \(2\nu \beta \beta \) decay of \(^{100}\)Mo with the full data set of the NEMO-3 experiment corresponding to a 34.3 kg\(\times \)year exposure are presented. The summed energy of two electrons, the single electron energy and the angular distributions between the two electrons have been studied with an unprecedented statistical precision (\(5\times 10^5\) events). The single electron energy distribution has been used to discriminate between different nuclear models providing direct experimental input into NME calculations. The HSD model is excluded with high confidence, while the SSD model is consistent with the NEMO-3 data. The corresponding half-life for the \(2\nu \beta \beta \) decay of \(^{100}\)Mo is found to be

$$\begin{aligned} T_{1/2} = \left[ 6.81 \pm 0.01\,\left( \text{ stat }\right) ^{+0.38}_{-0.40}\,\left( \text{ syst }\right) \right] \times 10^{18}~\text{ year }. \end{aligned}$$
(9)

Deviations from the expected shape of the \(^{100}\)Mo \(2\nu \beta \beta \) energy sum spectrum have been studied to obtain constraints on parameters for physics beyond the Standard Model. The most stringent upper limit to date has been obtained for the Majoron-neutrino coupling parameter \(\langle g_{ee}\rangle \) for the decay mode with a single Majoron particle emission and the spectral index \(n=3\). For other \(0\nu \beta \beta \) modes with two Majoron bosons emission a comparable sensitivity with the world’s best limits has been achieved. The most stringent constraints on the bosonic neutrino admixture and Lorentz invariance violation in \(2\nu \beta \beta \) decay have been set.