$\nu$DoBe -- A Python Tool for Neutrinoless Double Beta Decay

We present $\nu$DoBe, a Python tool for the computation of neutrinoless double beta decay ($0\nu\beta\beta$) rates in terms of lepton-number-violating operators in the Standard Model Effective Field Theory (SMEFT). The tool can be used for automated calculations of $0\nu\beta\beta$ rates, electron spectra and angular correlations for all isotopes of experimental interest, for lepton-number-violating operators up to and including dimension 9. The tool takes care of renormalization-group running to lower energies and provides the matching to the low-energy effective field theory and, at lower scales, to a chiral effective field theory description of $0\nu\beta\beta$ rates. The user can specify different sets of nuclear matrix elements from various many-body methods and hadronic low-energy constants. The tool can be used to quickly generate analytical and numerical expressions for $0\nu\beta\beta$ rates and to generate a large variety of plots. In this work, we provide examples of possible use along with a detailed code documentation. The code can be accessed through: GitHub: https://github.com/OScholer/nudobe Online User-Interface: https://nudobe.streamlit.app


Introduction
Neutrinoless double beta decay (0νββ) searches are the most sensitive probe of the violation of lepton number (LNV). The process takes place in a nucleus and converts two neutrons to two protons and two electrons, but no outgoing anti-neutrinos. So far only upper limits have been set on 0νββ rates [1][2][3][4][5][6][7][8][9][10][11], but next-generation experiments hope to make a first detection by probing 0νββ lifetimes reaching 10 27 -10 28 years [12][13][14][15][16]. A nonzero signal would have profound implications. It would imply lepton number is violated by two units. Following the famous black-box theorem [17] a positive observation of 0νββ would indicate that neutrinos are indeed Majorana mass eigenstates. Additionally, it would give important clues towards the mechanism of the neutrino mass and also be a strong hint for leptogenesis scenarios to explain the absence of anti-matter in our universe. Detailed reviews about 0νββ can be found in Refs. [18][19][20]. While tremendous experimental effort is going towards the first 0νββ detection, we must keep in mind that it is a complicated process involving particle, hadronic, nuclear, and atomic physics. The interpretation of a signal (or lack thereof) requires care. First of all, even if a nonzero decay rate is measured this does not immediately point towards the underlying source of LNV. Typically, 0νββ experiments are interpreted in terms of the effective neutrinos mass, m ββ , that enters through the exchange of light Majorana neutrinos (either active or potential sterile neutrinos). In a broad class of BSM models there can be competing mechanisms from the exchange of heavy particles such as heavy right-handed neutrinos or double charged scalars [21,22]. Studies of 0νββ-decay in multiple isotopes can be utilized to identify or reject certain types of BSM models [23][24][25][26].
As 0νββ decay is a low-energy process, the typical Q value of the reaction is a few MeV, it can be efficiently described through the use of effective field theory (EFT) techniques. In particular, EFTs provide a useful method to systematically categorize and order the possible LNV mechanisms [27][28][29][30][31][32]. If the origin of LNV happens at scales well above the scale of 0νββ decay, this allows for a model-independent description. EFT techniques are Figure 1: The EFT approach applied in νDoBe . Given a certain lepton number violating BSM-model, all 0νββ observables of interest can be calculated by following a subsequent chain of EFTs. In νDoBe the user is required to provide the matching conditions at either SMEFT or LEFT level. The calculations of 0νββ observables, plots, etc. are then automated by νDoBe . Figure modified from Ref. [32]. also very powerful in the description of the associated hadronic and nuclear physics that is required to connect data to the underlying LNV mechanism that is typically written at the level of elementary particles.
νDoBe is based on an end-to-end EFT framework, from particle to nuclear physics, that has been developed in Refs. [31,32] and is illustrated in Fig. 1. Assuming that LNV originates at scales well above the electroweak scale 1 , v ≃ 246 GeV, possible LNV mechanisms can be described by local effective LNV operators among Standard Model fields. The resulting EFT, often called the Standard Model Effective Field Theory (SMEFT), contains an infinite number of operators that can be ordered by their dimension. The higher the dimension of the operator, the more suppressed its low-energy effects are by additional powers of Q/Λ where Q is a low-energy scale and Λ ≫ v the scale where the SMEFT operators are generated. LNV operators appear with odd dimension [34] and νDoBe includes operators with dimension up-to-and-including dimension nine [35][36][37].
The SMEFT operators are then evolved from Λ to the electroweak scale, where heavy Standard Model fields (W, Z, Higgs, top) are integrated out and the theory is matched to an EFT that is often-called the low-energy EFT (LEFT). The LEFT operators are then evolved to the GeV scale and matched to hadronic LNV operators by using chiral EFT which provides an expansion in m π /Λ χ where m π denotes the pion mass, around 100 MeV, and Λ χ ∼ 2 GeV, the chiral-symmetry-breaking scale. The chiral perturbation theory Lagrangian is used to calculate the so-called neutrino transition operator that mitigates nn → pp + ee transitions. This step is non-trivial as strong interactions among nucleons lead to deviations from naive dimensional analysis expectations typically used to power count different contributions [38,39]. The neutrino transition operator is then inserted into nuclear wave functions to obtain the relevant nuclear matrix elements (NMEs) that can be combined with atomic calculations to obtain 0νββ decay rates.
The EFT approach calculates 0νββ decay rates in a systematic expansion (v/Λ) α (Λ χ /v) β (m π /Λ χ ) γ , where α, β, and γ are exponents that depend on the LNV mechanism under consideration. The decay rate is expressed in terms of phase space factors (PSFs), NMEs, low-energy constants (LECs), evolution factors, and the Wilson coefficient of the SMEFT operators, as outlined in detail in Refs. [31,32]. The actual procedure is, however, somewhat challenging and, frankly, tedious for the expert and non-expert alike because of the large number of steps and the sizeable number of QCD, nuclear, and atomic matrix elements. νDoBE is constructed to automate this procedure: the user chooses which LNV SMEFT operators appear at the high-energy scale Λ (or alternatively which LEFT operators at the electroweak scale), and the tool presents the resulting 0νββ decay observables. In this way, any high-scale model of LNV can be quickly analyzed.
In this work we describe νDoBe (Neutrinoless DOuble-BEta calculator) in detail and give some explicit examples of its use. Before doing so let us list the most important applications of νDoBe: • An automatic calculation of 0νββ decay rates and electron kinematics of all isotopes of experimental interest in terms of LNV SMEFT operators up to and including dimension nine. Alternatively, the decay rate can be expressed in terms of LNV LEFT operators. • The tool computes differential decays rates based on state-of-the-art information regarding renormalization-group factors, hadronic and nuclear matrix elements, and atomic phase space factors. • It is possible to choose various sets of NMEs (based on different many-body nuclearstructure methods) to assess the theoretical uncertainty of the predictions. Users can define their own set of NMEs if necessary. • The expression for 0νββ-decay rates includes recently identified short-range contributions associated to hard-neutrino-exchange processes [38,39]. The code uses current best estimates for the associated LECs but different values can be specified. • νDoBe can be used to generate interesting plots of half-lives or various different distributions as function of neutrino masses or LNV Wilson coefficients. This can be used to quickly analyse specific BSM scenarios. • For a given BSM model or a set of Wilson coefficients νDoBe can generate analytical expressions of the decay half-life and output them in latex (or html) form.
• Given an experimental bound on the half-life of certain isotope, νDoBe can be used to extract limits on the effective Majorana neutrino mass or higher-dimensional LNV LEFT or SMEFT operators.  mass m ββ using νDoBe. In particular, we will apply different sets of NMEs and study the potential impact of the recently introduced short-range contributions originating from hard neutrino-exchange [38,39] which were not included in Ref. [10]. We briefly describe the necessary inputs to the code here while a Jupyter notebook with all details is accessible on GitHub.
If only the standard mass mechanism is considered, the half-life is expressed as where the various M i denote NMEs, while g N N ν is an LEC associated with hard-neutrino exchange. To replicate the results from [10] we apply the same value for the phase space factor V 4 ud G 01 ( 136 Xe) = 1.458 × 10 −14 y −1 (the factor V 4 ud is included because of a slightly different definition of G 0k ). νDoBe has the NMEs split into various components, but most literature is based on the effective combination while setting g N N ν = 0. In Table 1 we show the different sets of NMEs studied here. 2 2 For some NME calculation the magnetic factor has been calculated incorrectly through gM = κ1gV , -5 - Figure 2: Limits on the effective Majorana mass m ββ obtained from the recent half-life limit by the KamLAND-Zen experiment [10]. On the left, the limits are shown without including the short-range contribution proportional to g N N ν as in [10], while on the right we do include it (see text for details). In the lower plots we show the current limit on the sum of the neutrino masses i m i ≲ 0.12 eV [65] as a gray band. The allowed parameter regions for normal (NO) and inverse (IO) mass ordering are displayed as red and green bands, respectively. They are obtained by variation of the unknown Majorana phases utilizing the plot_m_eff() function of νDoBe .
In general, to apply νDoBe we need the individual NMEs M F,GT,T,F sd . However, often only the combined NME M ν = − 1 g A M F + M GT + M T is given in the literature. Within νDoBe we can work around this by simply defining M GT = M ν while putting all other NMEs where κ1 is the isovector anomalous magnetic moment of the nucleon, instead of the correct gM = (1+κ1)gV . For more details see the appendices or Refs. [31,32]. If the individual NMEs M M M GT,T are given explicitly in the original publication, we correct Mν appropriately. In case that both CD Bonn and Argonne potentials are given we take the values corresponding to the CD Bonn potential. Table 2: Groups of LEFT operators that result in the same half-lives.
(including M F,sd ) to zero. By utilizing the internal get_limits() function we then arrive at the limits on m ββ in the range [36,146] meV. The individual limits on m ββ derived from each NME method are listed in Table 1 and the left column of Figure 2. Whenever NMEs calculated using both CD Bonn and Argonne nuclear potentials are available we choose to show the CD-Bonn ones. If, instead, the Argonne potential NMEs are considered the upper limit moves to 156 meV as in [10]. Now we want to study the effect of the additional short-range contribution proportional to g N N ν which has not been considered in [10]. For those NME calculations where the corresponding short-range NME M F,sd is available we simply add this contribution. However, for many NMEs no short-range contributions have been calculated. Here, we take a very simplistic approach by approximating the short-range NME M F,sd from the ratios M F,sd /M ν found in those NME calculations with M F,sd available. In 136 Xe we find that |M F,sd /M ν | > 0.2, hence, we take M F,sd = − 1 5 M ν as a lower limit approximation of the resulting half-life (Note, again, that g N N ν is negative and, hence, the contribution of a negative M F,sd is constructive). In this approximation we arrive at limits on m ββ in the range [25,68] meV which are more stringent, illustrating the importance of fixing the shortdistance contributions. The limits from all different NMEs with short-range contributions are shown in the right column of Figure 2 and Table 1.
We stress that these limits are only indicative as the short-distance contributions must be calculated for each nuclear many-body method consistently with the long-distance contributions, see Sect. (5.3.5) for a more detailed discussion. This has recently been done in Refs. [66,67].

Higher Dimensional Mechanisms
LEFT: We can use νDoBe to put limits on the higher-dimensional LNV operators. Considering operators from the low-energy EFT we can put these into 13 different groups of operators that result in the same half-lives (see Table 5 in Appendix A) because of parity conservation in QCD. For convenience we summarize these operator groups here again in Table 2. We require the corresponding NMEs to calculate the limits on the higher dimensional (d ≥ 6) operators of interest, which means we can only use a few NME calculation methods namely those from the shell model [64], QRPA [51] and IBM2 [59]. The resulting limits on the higher dimensional operators are given in Figure 3. The limits are extracted at the SMEFT → LEFT matching scale m W and acquire a broad range of values from V , and C (9) V ), thus illustrating the fact that different operators can have very different 0νββ efficiencies.
We can estimate the scale of new physics Λ through with d ∈ [6,7,9] being the operator's LEFT dimension and the Higgs vacuum expectation value (vev) v ≃ 246 GeV. This translates to limits on the scale of new physics in the range of ∼ 2 TeV -20 PeV. SMEFT: Alternatively, we can use the new limit on the half-life in 136 Xe from KamLAND-Zen to provide limits on the dimensionfull Wilson coefficients of the SMEFT operators. These limits are given in Fig. 4 and correspond to scales in the range of ∼ 1 TeV -400 TeV For more details, we refer to the provided Jupyter notebook on GitHub as well as Section 5.10.1.

The Minimal Left-Right Symmetric Model
The minimal left-right symmetric model (mLRSM) [68][69][70][71] is a well-studied extension to the Standard Model. It enlarges the Standard Model gauge group by adding the right-handed SU (2); hence, the mLRSM gauge group reads SU At the same time, this additional symmetry requires the existence of new fermionic and bosonic degrees of freedom. Typically, the mLRSM incorporates two scalar triplets ∆ L ∈ (1, 3, 1, 2) and ∆ R ∈ (1, 1, 3, 2) and a scalar bidoublet Φ ∈ (1, 2, 2 * , 0). Fermions are grouped into leftand right-handed SU (2) L,R doublets requiring the introduction of right-handed neutrinos In this section, we revisit the mLRSM as studied in [25,32] (conventions are equal to [25]) utilizing νDoBe and the shell model NMEs [64] provided therein. Lepton number violation in the mLRSM is introduced via Yukawa interactions, When the neutral components of the scalar multiplets acquire non-zero vacuum expectation values, they give rise to neutrino mass matrices (3.8) The ratio of the bidoublet's vevs ξ = κ ′ κ describes the mixing between the left-handed and right-handed W -bosons. When the right-handed triplet gains a non-zero vev, the mLRSM -9 -gauge group is broken down to the Standard Model gauge structure. Matching the mLRSM onto the SMEFT results in the following LNV Lagrangian, where H is the Standard Model Higgs with the vev v = κ 2 + κ ′ 2 and the SMEFT Wilson coefficients are given by (3.10) The Wilson coefficients are completely described by fixing the vevs of the triplet scalars v L,R , the heavy right-handed triplet mass m ∆ R , the heavy neutrino masses m ν R,i , i = 1...3, the lightest active neutrino mass m ν min , the neutrino mixing matrices U L,R , the complex vev phases θ L and α as well as the left-right mixing parameter ξ.
To apply νDoBe , we need to translate these Wilson coefficients to the correct operator basis (see Appendix A) via LeudΦ , C With the help of νDoBe we can use these Wilson coefficients to generate a SMEFT model class to study different parametric scenarios of mLRSM realizations in Python. As an example,  (3.13). The shaded areas show the Majorana mass m ββ arising from the light-neutrino-exchange mechanism alone, while the scattered points show the resulting effective mass for the different mLRSM realizations when varying the unknown Majorana phases, θ L and α. Both the normal (blue) and the inverted (red) mass orderings are displayed. The gray area shows the current limit given by KamLAND-Zen.
we revisit four different realizations of the mLRSM studied in Ref. [32], given by and by setting the left-handed triplets vev v L and the ratio of the bidoublet vevs ξ to one  In Figure 5 we show the evolution of the effective neutrino mass parameter, where M (ν) 3 is the NME for the light-neutrino-exchange mechanism (LνEM)(see [32]). Figure 6 shows similar information but now we normalized to m ββ , the effective mass in the -12 -LνEM. The scattered points show variations of the unknown Majorana phases in the neutrino mixing matrix as well as θ L and α. We observe that in all four parameter settings the case of inverted mass ordering is hardly influenced by the additional higher dimensional interactions introduced in the mLRSM. This would change for smaller values of the left-right symmetry breaking scale. On the contrary, for normal mass ordering the higher dimensional d ≥ 7 contributions do increase the expected decay rate when compared to the LνEM. In addition, the case of normal mass ordering does no longer show the typical funnel.
We refrain from a more detailed analysis of the parameter space as our main goal is to illustrate how νDoBE can be readily applied for 0νββ analyses of complicated models. Again, a detailed Jupyter notebook is accessible via GitHub. In this section we revisit the 0νββ decay induced by new leptoquark (LQ) interactions first studied in Ref. [72]. Assuming the Standard Model without right-handed neutrinos, one can add up to 10 scalar or vector LQs with non-trivial couplings to the Standard Model at tree level. These are summarized in Table 3. LQs with Q = ±1/3 and Q = ±2/3 then generate new decay channels for 0νββ decay different from the usually considered LνEM with the corresponding Feynman diagrams displayed in Figure 7. The matching procedure involves a minor subtlety when compared to the previous case study of the mLRSM as the LNV interactions relevant for 0νββ decay are generated when electroweak symmetry breaking introduces non-diagonal mixing between different LQs. The analysis is more straightforward if one matches the BSM model directly onto LEFT instead of SMEFT.

A Leptoquark Mechanism
After doing so, the parts of the low-energy Lagrangian relevant for 0νββ decay are given by [72] [72].
Hence, we obtain the following LEFT Wilson coefficients Now we can use the get_limits_LEFT() function of νDoBe to evaluate the limits on the parameters ϵ i , α S,V i , i ∈ [L, R] imposed by the recent KamLAND-Zen results when assuming only one parameter at a time to be non-vanishing. In doing so, using the IBM2 NMEs [59] provided with νDoBe we find that In comparison with the limits obtained in the original work [72] the above bounds are somewhat more stringent, which is mainly due the improved KamLAND-Zen results [10]. This aspect is partially compensated by using a different set of NMES and PSFs reflecting the improved status of theoretical calculations. A more detailed comparison with Ref. [72] can be found in the provided Jupyter notebook.
-14 -Neutrinoless double beta decay experiments provide the most stringent tests of lepton number violation. The experimental prospects are excellent and it will play an important role in the resolution of the neutrino mass puzzle and the search for beyond-the-Standard-Model physics in general.
That being said, 0νββ is a complicated process involving particle, nuclear, and atomic physics and its interpretation requires a great care. The main goal of νDoBe is to help the community with these difficulties. Being a low-energy process, effective field theory methods can be used efficiently to describe 0νββ for a large class of LNV scenarios stemming from possibly rich physics living at UV scales. νDoBe computes 0νββ (differential) decay rates in terms of Wilson coefficients of LNV operators constructed within the Standard Model Effective Field Theory (SMEFT). It uses the state-of-the-art renormalization-group evolution factors and QCD, nuclear, and atomic matrix elements and comes with a large number of built-in plotting and analysis tools. Naturally, νDoBe can help to easily extract limits on the m ββ parameter associated with the LνEM, but also on all the other beyond-the-Standard-Model Wilson coefficients of interest, imposed by current and upcoming experimental results. We hope it will be used by the community to quickly and accurately analyze 0νββ predictions and constraints related to beyond-the-Standard-Model models involving lepton number violation. We encourage users of νDoBe to share their BSM analyses with the community by adding a Jupyter notebook to νDoBe 's GitHub.
In the upcoming updates we aim to extend νDoBe in several directions, some of which are work in progress. Most importantly, we will include the possibility of relatively light sterile neutrinos, which means addition of effective operators arising in the SMEFT framework extended by three right-handed singlet neutrinos, so called νSMEFT. Several works have computed the impact of light sterile neutrinos on 0νββ observables, but the analysis is complicated and essentially done on a case-by-case basis [22,33,[73][74][75][76][77]. An automated tool would be very helpful in this regard. In addition to the currently available approximate calculations of the phase-space factors described in Appendix C we plan to add the option of employing the exact numerical solutions for the radial electron wave functions using a broader variety of nuclear potentials and including other subtle effects such as the electron screening along with and beyond the treatments adopted in Refs. [78,79]. Last but not least, we envision to incorporate predictions for and constraints from other lepton-numberviolating processes, such as kaon decays.

Online Tool
Finally, for quick analyses we created an online user-interface using Streamlit. The online user-interface is aimed towards delivering an easy and fast to use possibility for the most relevant use-cases like, e.g., calculating the decay observables (half-life, spectra, angular correlation) from a given model or studying limits on the different Wilson coefficients given half-life limits from experiments.

Documentation: νDoBe -A Python Tool for Neutrinoless Double Beta Decay
In the following we provide a detailed description of all the different functionalities of νDoBe including various explicit code examples. required) Additionally, Python 3 is required while 3.6 or higher is recommended. If you use Python 2 the code will give wrong results! Besides internal usage of pandas, νDoBe generally uses pandas' DataFrame class to output results in a table format. Pandas DataFrames provide a convenient framework for further analyses similar to numpy with additional features like, e.g., latex exports or plotting functions.

List of Functions and Classes
To ensure that there are no conflicts with other third-party python modules we recommend setting up a dedicated virtual environment when using νDoBe .

Setting up νDoBe
To use νDoBe simply download the code from GitHub and copy it into your project's main directory. Your projects directory structure should look something like this:

Physical Constants
Physical constants such as particles masses and mixing angles are defined in the constants.py file and can be changed if that is required.

Units
Generally, all dimensionful parameters in νDoBe are defined in GeV, i.e., we set GeV = 1, numerically. The only exception to this is the minimal neutrino mass in plotting functions which is taken in eV. To make it easier to not get confused with units, the numerical values of different energy units can be accessed in the nudobe.constants module.

Nuclear Matrix Elements
0νββ decay rates depend on a set of nuclear matrix elements that involve complicated nuclear structure calculations. As of today, the NMEs are computed within different nuclear many-body methods and the results, unfortunately, tend to differ. νDoBe comes with three sets of NMEs that one can choose from; the interacting boson model 2 ("IBM2") [59], the quasi random-phase approximation ("QRPA") [51] and the shell model ("SM") [64]. These particular sets of NMEs are chosen because within the corresponding method all NMEs required to analyze 0νββ from dimension-9 SMEFT operators have been computed. The corresponding files are stored in the NMEs/ folder as .csv files. Other NME approximation methods can be studied by, simply, adding a new .csv file to the NMEs/ folder. If you do so, please make sure to follow the definitions of the NMEs described in appendix D.
Within νDoBe there are two possibilities to set NME methods: 1. You can set the method right at the start when initiating a model as will be shown in the following section. Every calculation you will do within this model will then use the corresponding NMEs.
-19 -2. Alternatively, most functions have a parameter called method that will reset the NMEs temporarily when calling the function (see the functions definitions).
Additionally, within the model classes the NMEs can be accessed and changed in the model.NMEs dictionary.

Phase-Space Factors
νDoBe includes two approximation schemes for the calculation of PSFs and electron wave functions [86]. The PSF-scheme is defined when initiating a model class via the PSF_scheme parameter. Scheme "A" includes PSFs calculated from approximate wave functions in a uniform charge distribution while scheme "B" gives PSFs calculated exactly assuming a point-like nucleus. The Phase-Space factors for the approximation schemes A and B are stored in the PSFs/ folder as PSFs_A.csv and PSFs_B.csv. Again, similar to the NMEs you can use custom PSFs by replacing entries in the .csv files accordingly. Note, however, that this will only change the overall magnitude of the PSFs while the spectra and angular correlations are still calculated from the electron wave functions defined in the PSFs.py file and by the choice of the PSF-scheme. Additional methods for calculating the electron wave functions will be added in a future update. Note, however, that while different approximation schemes do have a noticeable effect on the magnitude of the PSFs, the general shape of the electron spectra and angular correlation is not expected to be influenced largely by a different choice of wave functions [86]. Within a model class, you can also access and change the PSFs by replacing the corresponding entries in the model.PSFpanda DataFrame.

Low Energy Constants
The choice of low energy constants (LECs) strongly influences the resulting half-lives. The values of the known LECs as well as order of magnitude estimates have been summarized in [32] and are shown in Table 4 for convenience.
The preset values of the LECs are defined in the constants.py file. Global changes of the LECs are best taken care of in this file. Locally, within the model classes we have implemented two generic settings that can be applied when doing calculations.

1.
In the first setting we use the known LECs and set the unknown LECs to null except keeping g N N 6,7 = g πN V =g πN V = 1 such that we still allow for short-range vector operators in the LEFT framework to contribute to the overall amplitude. This setting can be applied by calling unknown_LECs = False when initiating a model.

2.
In the second setting we take the unknown LECs to be equal to their positive order of magnitude estimate. This option can be chosen by applying unknown_LECs = True.
You can also manually adjust each LEC. For each initiated model class, the LECs are only meaningful within a given renormalization scheme that also affects the corresponding nuclear forces and thus the long-and short-distance NMEs. Refs. [66,90] outlined a strategy how this can be achieved. The idea is that the amplitude of the process nn → pp + ee, A ν , is observable and should therefore be regulator independent. While the process itself cannot be measured in any practical way, the amplitude has been computed in [66,90] by relating it to momentum integral of a known kernel (proportional to the neutrino propagator) multiplied by the generalized forward Compton scattering amplitude nnW + → ppW − . The lowand high-energy regime of this integral can be described model independently through, respectively, chiral EFT and the operator product expansion for perturbative QCD. The full amplitude is then obtained by interpolating between the two regimes using appropriate form factors for single-nucleon and nucleon-nucleon interactions. Once A ν is obtained at some kinematic point, it becomes possible to fix g N N ν for an EFT of nucleon-nucleon interactions using any regularization scheme. Crucially, the value of g N N ν is not fixed but depends on the applied EFT and the associated regularization scheme, however, the observable amplitude A ν is always correctly described. This procedure was followed in Ref. [67] which applied chiral EFT nucleon-nucleon interactions to extract g N N ν and then performed ab initio computations of the NMEs related to 0νββ of 48 Ca. The value of g N N ν in the table is the value taken from Ref. [67]. However, strictly speaking this value is not compatible with -21 -the NME sets of νDoBe for heavier isotopes which were obtained with different many-body methods. For that reason we have assigned a 50% uncertainty on the value of g N N ν .

Setting Up a Model
The goal of the tool is to directly connect BSM models that contain additional LNV sources to 0νββ decay rates and electron-kinematics measured in experiments. The main task of the user ("you", from now on) is to provide the matching relations between the specific BSM model to the LNV SMEFT or LEFT operators. Currently, the tool contains all ∆L = 2 operators relevant for 0νββ involving first-generations quarks and leptons up-toand-including dim-9 operators (all ∆L = 2 SMEFT operators have odd dimension [34]). The basis of dim-5, dim-7, and dim-9 SMEFT operators has been derived in Refs. [35][36][37]45] and we follow the notation of these references. A full list of the relevant operators is given in Appendix A.

EFT Model Classes
The code allows you to set up an EFT model consisting of LNV operators that trigger scale float Optional -Sets the scale of new physics Λ for the WCs. Needs to be set larger than or equal to 80 i.e. higher than the mass of the W boson. If scale>m_W νDoBe will automatically run the provided WCs down to m_W when initiating a model. name string Optional -Defines a name for the model. This will show up in plots.
unknown_LECs bool Optional -If set to True the unknown LECs will be set to their NDA estimates (see Table 4). If set to False the unknown LECs will be turned off i.e. set to 0.

Importing the Model Classes
To get started you first need to import the EFT classes which allow you to set up a model in either SMEFT or LEFT: [ ]: from nudobe import EFT

Initiating a Model
You can set up a model by defining the corresponding non-vanishing Wilson coefficients in a dictionary style.

List of Wilson coefficients
The relevant Wilson coefficients under consideration are given in Tables 5 (LEFT) and 6 (SMEFT). Within nudobe they are defined as dictionaries. The dictionary key for each WC is given in the column "Code Label" in each Alternatively, the list of WCs are also stored in the nudobe.constants module and can be accessed there.

Example: Light neutrino exchange mechanism in LEFT
You can generate a model class that represents the 0νββ-decay induced by the exchange of light Majorana neutrinos with m ββ = 1 eV using NMEs calculated in the IBM2 framework via inplace bool Optional -If True the resulting WCs after running will replace the models WCs.
The RGE running of the LEFT operators is given in [32] (Eqs. 14-16) while the RGEs for the dimension 7 SMEFT operators can be found in [42] (Eqs. [21][22][23][24]. The RGEs of the SMEFT dimension 9 operators will be included in a future update. Currently, νDoBe does not evolve the SMEFT dimension 9 operators when using the run() function. Therefore, for reasons of consistency, we recommend defining SMEFT models directly at the matching scale m W if both dimension 7 and dimension 9 operators are to be studied. Alternatively, you can use the SMEFT models .run() function to run the dimension 7 SMEFT operators from any arbitrary scale to the LEFT matching scale.
Example: RGE evolution of a specific SMEFT model.
[1]: from nudobe import EFT, constants from EFT import SMEFT from constants import * After importing the SMEFT class we define a model with two non-vanishing dimension 7 operators. To show the impact of the RGE running we first define the model class without setting a scale explicitly such that the SMEFT class assumes the scale is m W . This time we use the unit parameters defined in the constants module to define scales and masses.  method string Optional -Sets the NME calculation method. You can choose from "IBM2", "SM" and "QRPA". If None the models method will be used.
vary_LECs bool Optional -If set to True the unknown LECs will be varied within their NDA estimates (see Table 4). If set to False the unknown LECs will stay fixed. method string Optional -Sets the NME calculation method. You can choose from "IBM2", "SM" and "QRPA". If None the models method will be used.
Example: Light neutrino exchange mechanism in LEFT.
As an example we study once more the standard mass mechanism induced by the exchange of light Majorana neutrinos. We want to find the expected half-lives in the case of an effective Majorana mass m ββ = 100 meV: [1]: from nudobe import EFT, constants from EFT import LEFT from constants import * The resulting range of values can be used as an estimate for the theoretical uncertainty of the predictions arising from QCD matrix elements (that is, the LECs).
method string Optional -Sets the NME calculation method. You can choose from "IBM2", "SM" and "QRPA". If called via the functions module the preset value is "IBM2". If called via a model class he models method will be used if None.
decimal integer Optional -Sets the numbers of decimals used in rounding.
output string Optional -Get output in "latex" or "html" format.
unknown_LECs bool Optional -If set to True the unknown LECs will be set to their NDA estimates (see Table 4). If set to False the unknown LECs will be turned off i.e. set to 0.
PSF_scheme string Optional -Choose PSFs and electron wave functions -"A": approximate solution to a uniform charge distribution. "B": exact solution to a point-like charge to generate an expression for the decay rate either in latex or html format. This equation represents a matrix equation of the form and you can obtain the matrix coefficients from nudobe.functions.generate_matrix(WC, isotope = "136Xe", method = "IBM2", unknown_LECs = False, PSF_scheme ="A") model.generate_matrix(isotope, WC = None, method = None, unknown_LECs = False, PSF_scheme = None) Parameter Type Description isotope string Defines the isotope to be studied.

WC list
Optional -List of non-zero Wilson coefficients that should contribute to the half-life ["WCname1" ..., "WCnameN"]. If None the models non-zero WCs will be used.
method string Optional -Sets the NME calculation method. You can choose from "IBM2", "SM" and "QRPA". If called via the functions module the preset value is "IBM2". If called via a model class he models method will be used if None.
unknown_LECs bool Optional -If set to True the unknown LECs will be set to their NDA estimates (see Table 4). If set to False the unknown LECs will be turned off i.e. set to 0.
PSF_scheme string Optional -Choose PSFs and electron wave functions -"A": approximate solution to a uniform charge distribution. "B": exact solution to a point-like charge Both functions will only consider those Wilson coefficients which are non-zero in your model.
Example: The light-neutrino-exchange mechanism with an additional LNV right-handed current Let's assume a LEFT model with two lepton number violating terms in the Lagrangian [1]: import numpy as np from nudobe import EFT, constants from EFT import LEFT from constants import * Comparing the results on the decay rate T −1 1/2 R and R2 from the functions model.half_lives() and model.generate_matrix(), respectively, we see that they are equal up to some small accumulated rounding errors. yaxis string Optional -Choose from "t", "1/t", "m_eff" to get the half-life, inverse half-life or effective neutrino mass.
isotope string Optional -Defines the isotope to be studied. xscale string Optional -Sets scaling of x-axis e.g. "log" or "lin".
yscale string Optional -Sets scaling of y-axis e.g. "log" or "lin". experiments dictionary Optional -The same as limits but instead of giving the y-axis limit you set the half-life limit.
ordering string Optional -"both", "NO" or "IO" sets neutrino mass ordering.   In our half-life figure we want to show the current experimental limit from GERDA [2] in yellow. Hence, we choose 76 Ge as the isotope and define GERDA in the parameter [4]: isotope = "76Ge" experiments = {"GERDA" : {"limit" : 1.8e+26, "color" : "y", "linestyle" : "-", "linewidth" : 1, "alpha" : 0.5, "fill" : True, "label" : "GERDA"}} The gray area now corresponds to the plot above with just m ββ , i.e., it shows the expected half-life from the LνEM. We see that due to the new LNV operator, the 0νββ rates are significantly enhanced for the normal ordering for small m min but this effect is less pronounced in the inverted ordering. The yellow area shows the regions of the parameter space excluded by the GERDA experiment [2].
Additionally, instead of plotting the mass mechanism alongside our BSM model we can normalize the half-life generated by our model to the LνEM and plot it via -39 -From these plots we can easily see that while the case of inverse mass hierarchy is almost unaffected by the additional higher dimensional contribution, the normal mass hierarchy case results in significantly higher decay rates for minimal neutrino masses below 0.1 eV.

Scatter Plots
The above plots show all allowed values when varying the unknown phase(s) of the Wilson coefficient on the x-axis. If, instead, you want to additionally vary the unknown LECs, get a graphical representation of the probability distribution when doing a variation, or simply want a different design you can generate scatter plots through -40 - yaxis string Optional -Choose from "t", "1/t", "m_eff" to get the half-life, inverse half-life or effective neutrino mass.
isotope string Optional -Defines the isotope to be studied.
vary_phases string Optional -If True, vary the unknown phases of the x-axis WC.
vary_LECs string Optional -If True, vary the unknown LECs within their NDA estimate (see Table 4). xscale string Optional -Sets scaling of x-axis e.g. "log" or "lin".
limits dictionary Optional -Plots limits from experiments {Name : {limit, color, linestyle, linewidth, alpha, fill, la-bel}}. The Isotope is assumed to be the same as the isotope parameter for the plot.
experiments dictionary Optional -The same as limits but instead of giving the y-axis limit you set the half-life limit.
ordering string Optional -"both", "NO" or "IO" sets neutrino mass ordering. alpha_mass string Optional -Set alpha for mass mechanism if show_mbb=True.
alpha_cosmo string Optional -Set alpha for mass mechanism if cosmo=True. Example: light-neutrino-exchange mechanism with variation of g N N ν Again, we import the module and define the model by setting the NME method and WCs

Half-life Ratios
To distinguish different models among each other it can be interesting to look at the ratios of half-lives compared to a reference isotope [25], e.g. 76 Ge, You can return the ratios in all isotopes of interest in a tabular form by using -44 -model.ratios(reference_isotope = "76Ge", normalized = True, WC = None, method=None, vary_LECs = False, n_points = 100)

Parameter
Type Description reference_isotope string Optional -Sets the reference isotope for the half-life ratios.
normalized bool Optional -If True the ratios will be normalized to the standard mass mechanism R/R m ββ .
If None the models WCs will be used.
method string Optional -Sets the NME calculation method. You can choose from "IBM2", "SM" and "QRPA". If None the models method will be used.
vary_LECs bool Optional -If set to True the unknown LECs will be varied within their NDA estimates (see Table 4). If set to False the unknown LECs will stay fixed Here, the x-axis shows the range of ratio values while the y-axis gives the isotope mass number A.

Phase-Space Observables
The code can generate both the single electron spectra dΓ/dϵ defined in (C.3) as well as the energy dependant angular correlation coefficients a 1 /a 0 defined in (C.5) for any given model. To calculate the numerical values for a given normalized electron energy isotope string Optional -Defines the isotope to be studied.
method string Optional -Sets the NME calculation method. You can choose from "IBM2", "SM" and "QRPA". If None the models method will be used.

Plots are then generated through
Parameter Type Description isotope string Optional -Defines the isotope to be studied.
method string Optional -Sets the NME calculation method. You can choose from "IBM2", "SM" and "QRPA". If None the models method will be used. Example: The light-neutrino-exchange mechanism with an additional LNV right-handed current Let's look at our previous example again which combines a lepton number violating right-handed vector current with the standard mass mechanism. This time we choose m ββ = 0.1 eV and C (6) VR = 3 × 10 −7 . We initiate the model through [1]: import numpy as np from nudobe import EFT from EFT import LEFT [2]: WC = {"m_bb" : 0.1e-9, "VR(6)" : 3e-7 } model = LEFT(WC, method = "SM", name = r"$m_{\beta\beta} + C_\mathrm{VR}^{(6)}$") Then we define an array Ebar that represents the normalized kinetic energy of the -50 -outgoing electrons via E = (E e − m e )/Q. Additionally, we need to define a parameter e which helps us to avoid poles in the calculation of phase-space quantities g 0k and h 0k (see Appendix C) [4]: e = 1e-5 #avoid poles Ebar = np.linspace(0+e, 1-e, 10) We can then calculate the single electron spectrum for the angular correlation. In the latter case we refrained from normalizing the x-axis to present both options. While in case of the LνEM, the energy distribution for a single electron peaks, as expected, at Q/2, the exotic contribution proportional to C (6) VR features a dip around that value (see e.g. Ref. [25]), and therefore modifies the combined spectrum substantially, resulting in two equally high peaks on the sides of the distribution, i.e. for small and large energies of the electron, while there is a local minimum at Q/2. The negative angular correlation of the electrons obtained for the mass mechanism corresponds to the fact that they are preferably emitted back-to-back. The additional contribution from C (6) VR distorts the shape, such that there are small ranges of energies, for which the correlation becomes positive and the dip around the middle of the distribution is more pronounced. This is caused, again, by a fairly distinctive shape of the correlation for this higher dimensional 0νββ mechanism, see e.g. Ref. [25].
Measuring the spectra of the individual electrons as well as the angular correlation between the emitted electrons can help to identify BSM scenarios with non-standard vector interactions [25]. Experiments like NEMO-3 [4,11] and its next-generation upgrade Su-perNEMO [91] are equipped to resolve the individual electron kinematics such that both the spectra and the angular correlation can be studied. In fact, these observables are available for the ordinary double beta decay observed by NEMO-3 and also in this case they can be used to probe new physics [92].

Single operator dominance
Given experimental limits on the half-life of the 0νββ-decay one can translate these into limits on the different Wilson coefficients assuming the contributions from a single operator dominates the 0νββ amplitude. Under this assumption and using the LEFT operators, the decay rates scale as We can calculate limits on the different Wilson coefficients from a given half-life limit, t lim 1/2 ( A Z), through .
Additionally, one can translate limits on the different Wilson coefficients into limits on the expected new-physics scale Λ. Assuming that LEFT operators of dimension d scale as we can estimate the scale of new physics as Because of operator mixing induced by the running of the Wilson coefficients from m W → Λ χ it is important to clearly define the scale at which one wants to apply the assumption of single operator dominance. Similarly, when considering SMEFT the limits on the WCs can be obtained directly at the matching scale m W . Here, limits are taken with respect to the dimensionful WC and correspondingly the new physics scale can be estimated assumingC i = 1 as method string Optional -Sets the NME calculation method. You can choose from "IBM2", "SM" and "QRPA". If called via the functions module the preset value is "IBM2". If called via a model class he models method will be used if None.
groups bool Optional -If set to True operators with equal limits will be put into groups. Limits are then only shown for these operator groups. If False all operator limits will be shown.
scale string Optional -Only for LEFT models. Use "up" to get limits at m W or "down" to get limits at Λ χ .
unknown_LECs bool Optional -If set to True the unknown LECs will be set to their NDA estimates (see Table 4). If set to False the unknown LECs will be turned off i.e. set to 0. string Optional -Sets the NME calculation method. You can choose from "IBM2", "SM" and "QRPA". The preset value is "IBM2".
unknown_LECs bool Optional -If set to True the unknown LECs will be set to their NDA estimates (see Table 4). If set to False the unknown LECs will be turned off i.e. set to 0.
PSF_scheme string Optional -Choose PSFs and electron wave functions -"A": approximate solution to a uniform charge distribution. "B": exact solution to a point-like charge groups bool Optional -If set to True operators with equal limits will be put into groups. Limits are then only shown for these operator groups. If False all operator limits will be shown.

Two-operator scenario
If we drop the assumption of a single LNV operator the analysis becomes a bit more complicated. Assuming two non-vanishing real WCs C x,y with a relative phase ϕ the halflife is given as Now, instead of calculating the half life in dependency on C x,y , M xx,xy,yy and ϕ we can also solve for C y Given an experimental limit on the half-life, we can use this relation to find the allowed parameter space. Within nudobe this is done with the function functions.get_contours(WCx, WCy, half_life, isotope, method = "IBM2", phase = 3/4*np.pi, n_points = 5000, x_min = None, x_max = None, unknown_LECs = False, PSF_scheme ="A")

Parameter Type Description
WCx string Independent Wilson coefficient. Can be from SMEFT or LEFT.
WCy string Dependent Wilson coefficient to be calculated from WCx and the half-life limit. Should be from the same EFT as WCx.
half_life float Experimental half-life limits isotope string Defines the corresponding isotope the half-life limit is obtained from.
n_points float Optional -Sets the number of points to be calculated.
phase float Optional -Sets the relative phase ϕ.
x_min float Optional -Sets the minimum value of WCx to be studied.
x_max float Optional -Sets the maximum value of WCx to be studied.
unknown_LECs bool Optional -If set to True the unknown LECs will be set to their NDA estimates (see Table 4). If set to False the unknown LECs will be turned off i.e. set to 0.
PSF_scheme string Optional -Choose PSFs and electron wave functions -"A": approximate solution to a uniform charge distribution. "B": exact solution to a point-like charge The resulting contours can be plotted via

Parameter Type Description
WCx string Independent Wilson coefficient. Can be from SMEFT or LEFT.
WCy string Dependent Wilson coefficient to be calculated from WCx and the half-life limit. Should be from the same EFT as WCx.
limits dictionary Experimental half-life limits. Given as a dictionary of the type {Name : [half-life limit, isotope, color, label, linewidth, linealpha, linestyle, alpha]}. The color will decide the color of the plot while the label sets the label for the legend. Linealpha sets the alpha channel for the outer line of the contour. If it is not defined it will be set equal to alpha which defines the alpha channel of the filled area. If alpha is not set or set to None it is set automatically.

A List of Operators
In total, νDoBe contains 32 different LEFT operators and 36 different SMEFT operators that contribute to 0νββ. The relevant LEFT operators up to dimension 9 have been obtained in [31,32] and are summarized in Table 5. The general LEFT Lagrangian can be written as In the SMEFT we consider all LNV operators up to dimension 7 and all dimension 9 operators that match onto LEFT dimension 9 [37]. The relevant SMEFT operators are summarized in Table 6 and the SMEFT Lagrangian can be written as Therefore, we define the Wilson coefficients of LEFT as dimensionless numbers while the Wilson coefficients of the SMEFT are dimensionful, i.e., C d i,SMEFT ∝ Λ 4−d . We choose this somewhat awkward convention because it is done in most of the 0νββ literature.
In defining the operator basis, we use the following conventions: The partial derivative acting to the left and right is defined as The covariant derivative used to define the SMEFT operators (ignoring SU (3)) is given by with the Pauli matrices labeled via τ I , I = 1, 2, 3. The generators of SU (3), that are present in the definition of several LEFT operators, are denoted by where λ a , a = 1, ..., 8 represent the Gell-Mann matrices.

Class Op. Name Code Label
Op. Structure Half-Life Group SMEFT Dim.
V L , O V L , O SR , O V L , O SR , O V L , O SR , O T , O V L , O (9) 1L Table 6: Here we show the different SMEFT dimension 9 operators that we considered. Again, the "Code Label" column shows the name of the operators as used in νDoBe . The "LEFT Matching" column provides a list of LEFT operators each SMEFT operator matches onto with the precise matching given in Appendix B.

B Operator Matching
The transition amplitudes and correspondingly all 0νββ observables calculated by νDoBe are evaluated at LEFT level. Hence, models which are generated at SMEFT level first have to be matched down onto LEFT. At SMEFT level lepton number violation by two units without baryon number violation occurs only at odd dimensions [34]. Here we provide the explicit matching conditions at the matching scale E = m W for the mass-mechanism and the long-range operators and for the short-range LEFT operators. The dim-9 SMEFT operator O deQLH 2 D , additionally, matches onto the low-energy ∆ L = 2 operator O To study the relevance of O V R and an additional tensor contribution With the chiral power counting ϵ χ = m π /Λ χ and Λ χ ∼ m N ∼ 1 GeV we find that the tensor operator's ∂ ν e R γ µ ν C L + e R γ µ ∂ ν ν C L u L σ µν d R leading contribution to the decay -70 -amplitude is at O(ϵ 5 χ ) as it is proportional to the lepton momenta k ∼ ϵ 3 χ Λ χ , the neutrino momentum q ∼ ϵ χ Λ χ and requires the inclusion of NLO nuclear currents ∼ ϵ χ . We can therefore savely ignore it. With m d ∼ m u ∼ ϵ 2 χ Λ χ and the leading order contributions from C (6) V L ∼ Λ χ ϵ 2 χ and C (6) V R ∼ Λ χ ϵ 3 χ , the leading order contribution from the remaining parts of O Therefore, compared to the other dim-7 and dim-9 operators [32] O (7) V R ′ is suppressed by an additional factor of ϵ 2 χ and we ignore it.

C.1 Calculation of PSFs
The different phase-space factors G 0k that enter the half-life formula are defined as × h 0k (ϵ 1 , ϵ 2 , R) cos θ + g 0k (ϵ 1 , ϵ 2 , R) × p 1 p 2 ϵ 1 ϵ 2 dϵ 1 dϵ 2 d(cos θ), where ϵ i , p i denote the energy and momentum of the i−th. single electron, θ is the angle between the emitted electrons, we defined and h 0k , g 0k were introduced in Ref. [86]. The constants C k take care of the different conventions used in Refs. [86] and [32]. νDoBe takes care of this re-scaling internally, such that no action is required from the user and C k = 1 can be chosen when calculating PSFs to use in νDoBe . The quark-mixing parameter V ud is not part of G 0k in contrast to the usual literature convention. This is to account for the fact that higher dimensional typically do not exchange two W bosons. Instead, factors of V ud enter the sub-amplitudes A k (see Ref. [32]). νDoBe includes two approximation schemes for the calculation of PSFs and electron wave functions, i.e., an approximate solution to a uniform charge distribution (Scheme A) as well as an exact solution to a point-like nuclear potential (Scheme B).

C.2 Related Observables
Apart from the total decay rate, we are also interested in the differential rates, i.e., the energy and angular electron spectra [86]. The normalized single electron spectrum can be determined through

D Nuclear Matrix Elements
The different NMEs that enter the calculation of the sub-amplitudes A k are defined via g M (q 2 ) = (1 + κ 1 )g V (q 2 ), g P (q 2 ) = − 2m N g A (q 2 ) q 2 + m takes care of different conventions used in the common literature and [32]. νDoBe takes care of this re-scaling internally when importing NMEs from the corresponding .csv files. That is, the NMEs defined in the .csv files are multiplied by C sd when initiating a model class such that, similar to the case of C k in the PSFs, users can choose C sd = 1 when calculating NMEs to use in νDoBe .