This white paper is an outgrowth of a National Science Foundation (NSF) Project Scoping Workshop, the purpose of which was to assess the current status of calculations for the nuclear matrix elements governing neutrinoless doublebeta decay and determine if more work on them is required. The recent e↵ort to define the United States' role in ton-scale experiments, together with the conclusion in 2021 of the Department of Energy (DOE) Topical Collaboration on neutrinoless double-beta decay, made such an exercise extremely timely.
The main conclusions of the workshop can be summarized as follows:
• Neutrinoless double-beta decay of nuclei is an important window into the physics of neutrinos. It could be the first lepton-number violating process ever observed. As such, it would provide key insights into physics beyond the Standard Model-in particular how the matter-anti-matter asymmetry of the universe arose.
• Much progress on the theory of neutrinoless double-beta decay of nuclei has been made over the last five years.
An end-to-end set of e↵ective field theories (EFTs) that shows how to evolve the physics of the decay from the scale at which lepton number is violated (possibly much larger than TeV) down to the scales relevant for nuclei has been developed. Lattice quantum chromodynamics (LQCD) calculations of the process in the two-nucleon (NN) system are being set up, and the first ab initio many-body calculations of neutrinoless double-beta decay matrix elements M 0⌫ have been carried out.
• Much remains to be done before theory can successfully complement the large experimental e↵ort to observe neutrinoless double-beta decay and measure its rate. Both the accuracy and precision of LQCD and ab initio nuclear many-body calculations need to be improved if crisp statements about experimental observation or non-observation are to be made. The uncertainty in nuclear many-body calculations remains largely unquantified, making it di cult to interpret the significant di↵erences predicted by di↵erent approaches for the rate of neutrinoless double-beta decay in candidate nuclei.
• Uncertainty quantification is thus crucial to future progress. Better assessment of both the parametric uncertainty and the model uncertainty in predictions of neutrinoless double-beta decay matrix elements is needed.
The tools for the former exist, and methodology for the latter is under development. The ab initio calculation of a variety of nuclear observables related to neutrinoless double-beta decay can help establish and reduce the uncertainty in M 0⌫ that arises from the complexity of the nuclear many-body problem.
In recent years the search for new fundamental physics, for the forces and particles that underlie the Standard Model, for the explanation of the excess of matter over antimatter and similar mysteries, and for the sources both of symmetries and their violation, has moved increasingly to low-energy experiments. Among the most visible and promising are those that seek to observe neutrinoless double-beta (0⌫ ) decay, a process in which two neutrons inside an atomic nucleus turn into protons, emitting two electrons and no neutrinos. An observation of this process would show that lepton number (L) is not conserved and that the neutrino mass has a Majorana component, implying that the mass eigenstates are self-conjugate [1]. Observation of 0⌫ decay would thus provide crucial information about neutrino mass generation [2][3][4], and suggest that the matter-antimatter asymmetry in the universe originated in leptogenesis [5]. The major implications of an observation made the construction of a ton-scale 0⌫ -decay experiment the top priority for new projects in the 2015 NSAC Long Range Plan [6], which set the decadal priorities for nuclear physics. The anticipated investment is in the range of 250-400 million dollars.
Smaller experiments already put stringent limits on the decay rate [7][8][9][10][11][12][13][14][15][16][17][18], e.g. T 0⌫ 1/2 > 2.3 ⇥ 10 26 yr for the decay of 136 Xe [19]. The next very few years will see stricter limits from experiments-such as LEGEND-200, CUORE, KamLAND-Zen 800, and SNO+-that are currently operating or under construction. On a slightly longer time scale, ton-scale experiments [19][20][21][22][23][24] based on 76 Ge, 100 Mo, 136 Xe, and perhaps other isotopes will come on line. The goal of these large experiments is the ability to detect any decay caused by the exchange of light Majorana neutrinos if the neutrino mass hierarchy is inverted (i.e., if the neutrino with the largest electron-flavor component is the heaviest), as well as increased sensitivity to decay caused primarily by the exchange of other still-hypothetical particles.
In order to extract the e↵ective light-neutrino Majorana mass m ⌘ | P i U 2 ei m i | (with m i the mass of the neutrino mass eigenstate i and U ei the elements of the Pontecorvo-Maki-Nakagawa-Sato matrix) from any of these impressive experiments, one needs nuclear matrix elements (denoted by M 0⌫ ) of the decay operators. The degree to which ton-scale experiments will be sensitive to decay caused by the exchange of inverted-hierarchy light neutrinos depends on these nuclear matrix elements, as does the the extent to which experiments in more than one isotope will prove useful. The nuclear matrix elements su↵er at present from sizable uncertainties [25]. Their accurate computation, with a quantified uncertainty, is therefore a task of the greatest importance.
The need for precise nuclear matrix elements is in fact more general than the notion that the exchange of light Majorana neutrinos causes 0⌫ decay. That idea is based on the assumption that lepton-number violation (LNV) originates at very high energies and manifests itself in the decay through the "high-scale seesaw," which leaves Majorana neutrino masses as its only remnant at low energies. If that is indeed the case, 0⌫ decay and neutrinooscillation experiments will together tell us most of what we can learn. High-scale LNV is only one scenario, however, and even in the restricted class of seesaw models, it applies only if the new particles are very heavy right-handed neutrinos. In many Beyond-the-Standard-Model (BSM) scenarios, other lower-scale sources of LNV can also induce 0⌫ decay. In left-right symmetric models, for example, heavy neutrinos and charged scalars with TeV-scale masses can be exchanged. In other scenarios there may be light right-handed (sterile) neutrinos with masses much lower than the electroweak scale. The large number of ways in which lepton number could be violated (see, e.g., Ref. [26] for a review) means that ton-scale searches for 0⌫
decay have a significant discovery potential beyond the invertedhierarchy high-scale seesaw. Each kind of LNV leads to its own set of transition operators, the nuclear matrix elements of which must be calculated. If the calculations are su ciently accurate, we can assess the sensitivity of the generation of experiments now coming online to various kinds of LNV. We can also provide a subsequent generation of experiments with information on how best to narrow the range of possibilities for LNV and neutrino-mass generation through measurements of single-electron spectra, electron angular distributions, and the isotopic dependence of the decay rate.
The ability to compute all the relevant nuclear matrix elements requires work at widely separated energy scales, from the high energies at which LNV originates all the way down to nuclear energies, and the ability to bridge the scales. EFT provides the bridge by expanding observables and Lagrangians in the ratios of the important energy scales. In reality the calculation is done via a series of EFTs-a connected set of bridges rather than a single one; see Fig. 1 for an illustration. The SM EFT allows us to encode the e↵ects of di↵erent LNV mechanisms in operators involving neutrinos, electrons, and d and u quarks, thereby taking us from the TeV scale to the scale of quark confinement at around 1 GeV. Converting these operators into hadronic operators that are organized through chiral perturbation theory requires non-perturbative input from LQCD. Chiral-perturbation-theory operators are then used to derive operators in a nucleons-only Hilbert space; following that step, the operators can be used in many-body calculations of nuclei. In combination, the bridges deliver us a separate set of chiral EFT nn ! pp transition operators for each LNV source that are to be used in nuclear many-body calculations. The combination of SM EFT, LQCD, chiral EFT, and ab initio (first-principles) nuclear many-body methods, each of which has the ability to control uncertainty, therefore provides a path-the only path, in fact-toward the reliable estimation of uncertainties in M 0⌫ .
Chiral EFT is key to the progress made to this point, and to future e↵orts to quantify uncertainties. Chiral EFT [27][28][29][30] is the extension of chiral perturbation theory to the few-nucleon problem. Just as with chiral perturbation theory, chiral EFT is organized as an expansion in powers of p/⇤ or m ⇡ /⇤, where p is a typical nucleon momentum, m ⇡ is the pion mass and ⇤ is the theory's "breakdown scale" of about 500 MeV. But chiral EFT is not a perturbative theory, because it has to account for nuclear binding. Although discussions of exactly how to do that continue (see, e.g., [31]) chiral EFT has the virtue of delivering consistent nuclear forces and 0⌫ operators up to a given order in the chiral EFT expansion. Even better, these operators include the consequences of QCD's chiral symmetry, e.g., connections between pionic operators and the axial current that governs beta decay. Perhaps most significantly for the purposes of this document, chiral EFT permits estimation of the uncertainty associated with the model of the nuclear force and the interactions that govern 0⌫ decay. A kth order chiral EFT calculation should have a fractional error of O({p, m ⇡ } k+1 /⇤ k+1 ). It follows that di↵erent implementations of chiral EFT-di↵erent orders of the calculation, di↵erent regulator choices-should give answers that are consistent with one another once this error estimate is taken into account. Bayesian techniques have recently been employed to quantify this error [32], and show that-provided the chiral EFT calculation is implemented carefully-the error estimate provides a good account of the predictive accuracy of chiral EFT in light nuclei [33,34]. Chiral-EFT forces and operators therefore provide the starting point for ab initio calculations that use the nuclear many-body methods described below.
The nuclear-theory community has made significant progress, at all the levels in this tower of EFTs, toward more accurate calculations of M 0⌫ . But this progress has in part served to confirm that there are O(1) uncertainties in M 0⌫ . These uncertainties (unless reduced) will prevent us from learning about the sources of LNV, even if several experiments detect the process.
There is therefore still much to do. In particular:
• The 0⌫ transition operators used in nuclear-structure physics are now written in terms of "low-energy constants" (LECs) that multiply terms in the chiral-EFT Lagrangian that is used at the hadronic scale. In chiral EFT, the LECs multiplying the terms at the lowest orders are thus the most important. Previously unrecognized LECs associated with zero-range nn ! pp transition operators appear even at leading order in the 0⌫ piece of the chiral Lagrangian, for both light-Majorana neutrino exchange [36] and TeV scale LNV [35]. We must improve our knowledge of these LECs, both by relating 0⌫ decay to other I = 2 processes and by direct calculation within LQCD.
• To use the results of EFT and LQCD in the computation of nuclear matrix elements-that is, to use the chiral-EFT Hamiltonians and transition operators that these methods supply in many-body calculations-we need to improve ab initio methods. The improvement will involve an increase in accuracy, the use of a wide range of chiral-EFT Hamiltonians (to allow uncertainty quantification), and a careful analysis of the way such methods employ the EFT operators. The first two of these will require, in addition to analytic work, more e cient use of our best supercomputing resources. Existing codes and their extensions will need to be reworked to leverage accelerators such as GPUs. Benchmarking with methods that are known to give very accurate results (so-called "quasi-exact" methods that have thus far been restricted by complexity to light nuclei) is also important.
• At both the hadronic and nuclear scales, we need a consistent and unified quantification of uncertainties. We must be able to both propagate parameter uncertainties to observables and to account for and disentangle deficiencies in our calculations. The innovative use of Bayesian methods will be essential.
In short, the framework developed in the last few years to combine LQCD, EFT, and ab initio nuclear structure is not yet e cient enough to allow a genuine assessment of uncertainty. To be of real use in the search for new physics, all three ingredients must be improved in the coming decade and made more computationally e cient; their uncertainty also needs to be reliably addressed. But these kinds of intelligent improvements will not, on their own, be enough: increased access to computing resources and dedicated exascale allocations will also be important.
The NSF Project Scoping Workshop that led to this white paper was organized by Jon Engel, Witek Nazarewicz, and Daniel Phillips, and held virtually on January 31 and February 1, 2022. After reviewing the experimental and theoretical status of the field and discussing recent developments, the attendees set forth the challenges to interpreting experimental results, discussed ways to address those challenges, mapped out a path forward, and planned this report. More details about the workshop (list of participants, program, presentations) can be found on its website.
A recent Snowmass white paper [37] provides a particle-physics perspective on these issues. The Project Scoping Workshop, and this report, are focused more on the nuclear-theory aspects of 0⌫ calculations. The nuclear-and particle-theory for 0⌫ decay is also reviewed-together with the experimental situation-in Ref. [38]. Our workshop and this report di↵er from these papers in the detail in which they discuss modern methods and the accompanying uncertainty quantification that will allow a meaningful error bar in the prediction of M 0⌫ for any particular BSM mechanism.
In the next section we provide a summary of the current state-of-the-art in both the nuclear-physics aspects of M 0⌫ (Sec. II A) and uncertainty quantifcation (UQ) in nuclear theory (Sec. II B). Section III then discusses the innovations and calculations that are needed to advance the nuclear theory of M 0⌫ , while Sec. IV describes a plan to quantify uncertainty in those calculations. Because of the significant emphasis on UQ for 0⌫ matrix elements in this report Secs. IV and II B are quite detailed and explicit about how we think that UQ can be carried out. We close in Sec. V with a summary of the theory advances and collaborative structures that are needed in order to establish precise and accurate calculations of neutrinoless double-beta decay.
Much of the current state of the art in the computation of M 0⌫ arose from work in the DBD Topical Theory Collaboration. LQCD, EFT, and nuclear many-body methods all played a role in the multi-scale problem. We discuss recent developments in each of these areas.
EFT. Before the Topical Collaboration, the connection of nucleon operators with fundamental sources of leptonnumber violation tended to be ad hoc, with BSM models analyzed individually, and unsystematically. The first application of the framework of chiral EFT to the problem, for 0⌫ decay induced by heavy-particle exchange, appeared in Ref. [39]. In the last few years, work of this kind has grown much more systematic. References [35,40,41] systematized the work of Ref. [39], showing how the parameters that determine the rates of very heavy-particle leptonnumber violating physics work their way down into nucleon-level operators. At around the same time, Ref.
[42] treated light-neutrino exchange, showing that working to N 2 LO requires "non-factorizable" diagrams (those that cannot be broken in two by cutting the line representing the exchanged neutrino) that had never been considered before. Shortly after that, researchers made the surprising discovery [36] that a contact interaction, representing the e↵ects of high virtual-neutrino momenta that are integrated out of the chiral EFT, occurs at leading order. Though the coe cient of the contact operator was initially unknown, it was later determined approximately through a resonancemodel-based interpolation between perturbative QCD and low-energy pion and nucleon dynamics. [43,44]. For the first time, nuclear many-body computations of M 0⌫ in the nuclei used in experiments are taking the contact term into account. So far it has caused all ab initio matrix elements to increase.
LQCD. The hope is that LQCD will soon be able to directly supply the coe cient of the aforementioned contact terms, as well all other relevant LECs. In the last few years, the field has made significant progress toward that goal. A contribution to 0⌫ decay with TeV-scale LNV is produced by the exchange of BSM heavy particles between two pions, each of which are then absorbed by protons as they turn into neutrons. The exchange between these virtual pions is easier to compute with LQCD than the direct exchange between nucleons, and in recent work the dependence of the resulting 0⌫ nucleon-level matrix elements on parameters that specify BSM models has been calculated [45,46]. Pionic matrix elements in the light-neutrino exchange scenario have also been computed in LQCD, and the corresponding LECs in chiral perturbation theory have been constrained [47,48]. We anticipate progress toward direct calculations of nn ! pp matrix elements and are developing the formalism for constraining contact LECs from future LQCD calculations [49,50].
Nuclear Structure. At the nuclear-structure scale, recent progress has been mostly in applying newly developed non-perturbative ab initio many-body methods to decay. Such methods start with interactions and operators determined from QCD and/or fit to data in very light nuclei (A = 2 , 3, or 4), and then produce (approximate) solutions to the Schrödinger equation in heavier nuclei. Three distinct ab initio methods have been applied together with chiral-EFT interactions to the heavy open-shell nuclei of interest for 0⌫ experiments. The first two, the In-Medium Generator Coordinate Method (IM-GCM) and the Valence Space IMSRG (VS-IMSRG) are variants of the In-Medium Similarity Renormalization Group (IMSRG), an approach in which one uses renormalization-group flow equations to decouple a predefined "reference" state, ensemble, or subspace from the bulk of the many-body Hilbert space. The third method is Coupled Cluster (CC) Theory; it uses an ansatz for the ground state in which FIG. 2. The light-neutrino-exchange M0⌫ for the transition 48 Ca ! 48 Ti, computed in various approaches. The four rightmost values, in green, all result from the the same chiral-EFT interaction. References: EDF [51,52], IBM [53], QRPA [54], SM-pf [55,56], SM-sdpf [57], SM-MBPT [58], RSM [59], QMC+SM [60], IM-GCM [61], VS-IMSRG [62], CCSD,CCSD-T1 [63].
particle-hole excitation operators are exponentiated before being applied to a Slater determinant.
All three of these methods, along with many older and more phenomenological schemes, have been applied to the computation of M 0⌫ for light-neutrino exchange in 48 Ca, the lightest isotope that can be used in an experiment. Figure 2 displays the compiled results. Those produced by the methods just described-the "most ab initio"-are shown in green on the right of the figure. The uncertainty range is more significant for these than for other methods, but still omits most systematic error.
The next section reviews the state of the art in uncertainty quantification. The rest of this document then discusses both the ways in which physics methods can be improved in accuracy, and the ways in which remaining uncertainty in their predictions for M 0⌫ can be reliably estimated.
In 2011 Physical Review A published an Editorial that stated ". . . there is a broad class of papers where estimates of theoretical uncertainties can and should be made. Papers presenting the results of theoretical calculations are expected to include uncertainty estimates for the calculations whenever practicable." [64]. Uncertainty quantification is crucial for calculations of 0⌫ -decay in nuclei. The planning and-eventually, we hope-the interpretation of 0⌫decay measurements requires that theorists deliver not just a expectation value for M 0⌫ , but also an uncertainty that represents the range of probable values that matrix element can take and does so in a statistically meaningful way. The goal of the uncertainty quantification (UQ) is not a precise evaluation of whatever is missing from the calculation. Quoting again from Ref. [64]: "The aim is to estimate the uncertainty, not to state the exact amount of the error or provide a rigorous bound."
UQ in nuclear-physics calculations pre-dated that Editorial but standard regression analysis was prevalent for many years [65]. Since then, nuclear-theory UQ has become much more sophisticated. This progress has taken place on several fronts.
The first, and most straightforward kind of UQ, is the estimation of error bars on the parameters ✓ in the nuclearphysics model being employed, then the propagation of those uncertainties-including their correlations-to model predictions. An early example of such an e↵ort is the estimation of the parameters in a sophisticated nuclear energydensity functional [66]. There are many recent examples of nuclear-structure calculations that do this, but a particularly striking one from the ab initio world constrained the parameters of nuclear forces using data from light nuclei and propagated the resulting uncertainties to predictions for properties of 208 Pb [67].
Three pieces of theory technology are commonly employed in such studies:
• Bayes' theorem, which relates the multi-dimensional posterior probability density p of the model parameters ✓ to the data, y, used to constrain those parameters, according to:
where p(✓) is the a priori distribution of the parameters ✓.
• A method by which a representative set of samples of the posterior probability distribution p(✓|y) can be obtained. Markov Chain Monte Carlo (MCMC) sampling is commonly employed, and was combined with a technique called "history matching" in Ref. [67]. We note that with such a set of samples in hand, it is conceptually straightforward to obtain a predictive probability distribution for, say, M 0⌫ . That distribution is found by repeated forward evaluation of the model for M 0⌫ at the di↵erent parameter values ✓ in the set of samples.
• Emulators that allow rapid evaluation of the (approximate) model at di↵erent values of the parameters ✓. This makes practical the computation of the likelihood p(y|✓) within whatever sampling framework is chosen.
Nuclear theorists have also become more attuned to the imperfections in their models. The inclusion of a "model discrepancy" term in the analysis of data is known to be crucial for accurate parameter estimation [68,69]. This means that one must admit that not just data, but also calculations, have imperfections that may cause a mismatch between theory and experiment. This idea can be formalized as,
where the last two terms encode, respectively, the experimental error (often taken to be independent at the x's corresponding to di↵erent data y) and the theory uncertainty (which is almost certainly not independent, i.e., we expect to be correlated across di↵erent x's). Significant e↵ort has gone into building models of y th for EFT calculations [32,70]; since EFT methods are characterized by a systematic expansion in a small parameter, one can predict how they will fail and so write down candidate functional forms for y th . But, even when such control is not available, model defects can still be productively introduced, e.g., Gaussian processes can be used to model the discrepancy between density-functional-theory calculations of masses and experimental data thereon [71]. Ultimately, though, the complex dynamics of nuclei means that di↵erent theoretical models will be employed to describe them. This diversity of models is advantageous because the methods have complementary strengths but also di↵erent systematic model discrepancies. This becomes a virtue by exploiting the third area of progress, which has been in the use of forms of Bayesian Model Averaging (BMA) or Bayesian Model Mixing (BMM) to incorporate insights from di↵erent models into a unified prediction in a statistically rigorous way. BMM can only be done reliably if individual models M k have had their uncertainties quantified in the ways described in the previous two paragraphs. Once that has taken place, the predictions of those models for the observable of interest y ⇤ can be weighted according to "scoring criteria":
Here, p(y ⇤ |y, M k ) is the posterior for the observable y ⇤ , given the data y, in a particular model M k , and w k (y ev ) is a weight that is determined by the model's performance on a target (or evidence) dataset, y ev . While in the BMA expression (3) the weights are constant across the domain, in the more advanced BMM they can also depend on x. We pause here to make two crucial points:
• The data, y ev that are used to assess the aspects of model performance that are pertinent for predicting y ⇤ need not be the same as the data set(s) used to calibrate the models. Ideally, the data y ev will be chosen because they are understood to be, or analyzed to be, a good proxy for the quantity of interest, y ⇤ , i.e., models' ability to predict y ⇤ is highly correlated with their ability to predict whatever observables are selected to be part of y ev .
• We note that that performance will almost certainly be addressed within the context of the model discrepancy y th of each model, and hence an understanding of those model discrepancies plays a critical role in betweenmodel UQ.
Early nuclear-physics applications of model averaging can be found in Refs. [72,73].
In 2020, the Bayesian Analysis of Nuclear Dynamics (BAND) collaboration [74] began its e↵ort to lower the barrier for nuclear theorists to perform all three of these types of uncertainty quantification. A particular interest within BAND is methodological work on BMM. The main product the collaboration seeks to deliver is software packages and use cases that facilitate emulation, model calibration, and model mixing. The goals of BAND are described in Ref. [75].
Within the context of the 0⌫ Topical Collaboration, some e↵orts were proposed to quantify the uncertainties in calculations of ab initio M 0⌫ . However, these e↵orts were limited by the ability to rapidly evaluate these matrix elements for the nuclei of interest in 0⌫
experiments. This made it di cult to even accomplish the first, parametric, kind of UQ. The determination of model discrepancy for the di↵erent many-body methods employed in 0⌫ studies remains a topic of forefront research, see Sec. IV.
Future e↵orts in the community to deliver reliable 0⌫ nuclear matrix elements will probably focus on advancing LQCD and EFT calculations of the underlying matrix elements in the few-nucleon sector, ab initio nuclear many-body calculations that use the LQCD and EFT input in experimentally-relevant isotopes. We discuss these subjects in this section and lay out a path for rigorous UQ in the next section.
The goal of a combined EFT and LQCD e↵ort in the 0⌫
program will be the identification and computation of the LECs multiplying nn ! pp transition operators in chiral EFT. Because BSM LNV interactions originally involve leptons and quarks, one has to evaluate matrix elements of quark operators in hadronic states in order to link LNV parameters such as m to the LECs. The program of constructing consistent and predictive nuclear EFTs has a long history and a recent review summarizing its status and prospects can be found in Ref. [30].
As mentioned in Sec. II A, analyses of the nn ! pp 0⌫ amplitude revealed that new nn ! pp contact interactions are needed at leading order in chiral EFT, even for light Majorana-neutrino exchange [36,42]. The associated LEC, called g NN ⌫ , is not determined by symmetry considerations or experiment (at least not in a straightforward way) and so must be obtained theoretically. The LEC g NN ⌫ has been studied so far by applying both large-N c and dispersive methods, while LQCD studies require methods that are still under development. Large-N c QCD [76] relates g NN ⌫ to LECs that can be extracted from the charge-independence-breaking combination of nucleon-nucleon scattering lengths in the 1 S 0 channel. Meanwhile, the dispersion-theory approach, inspired by the Cottingham formula for electromagnetic hadron masses [77,78], leads to a prediction for the nn ! pp amplitude near threshold, from which g NN ⌫ can be extracted in any EFT regularization and renormalization scheme, including those used in nuclear manybody theory. (See Ref. [79] for an early use of this nn ! pp input in ab initio nuclear many-body calculations.) The main uncertainty in this approach comes from inelastic intermediate states that can appear between the two insertions of the weak current (e.g. NN⇡ states). The existing estimates of g NN ⌫ can be improved by analyzing suitable I = 2 observables, thus anchoring g NN ⌫ to data. Finally, as discussed below, LQCD can play a major role in a first-principles determination of g NN ⌫ [50,80,81]. The leading-order LECs associated with TeV-scale LNV operators are currently completely unknown. Their determination will be possible through the use of dispersion-theory techniques similar to those developed in Ref. [77,78], as well as by through a direct calculation in LQCD. In fact, direct LQCD calculations can in principle determine the entire nn ! pp amplitude. (See Refs. [82,83] for recent reviews of the role of LQCD in constraining nuclear observables). The interplay between LQCD and EFT is symbiotic: On one hand, matching EFT and LQCD will enable an assessment of the theoretical foundation of nuclear EFT and a calibration of its truncation scheme. On the other hand, EFT descriptions allow better quantification of the systematic uncertainties in LQCD calculations, providing reasonable extrapolation forms for taking continuum and infinite-volume limits. Furthermore, in order to play a role in the 0⌫ program, LQCD calculations need to be performed at quark masses that are su ciently close to the physical values to allow reliable extrapolations to the physical point. Such extrapolations rely on EFTs, which in turn rely upon LQCD input to determine the relevant LECs. Thus, an interplay between LQCD calculations of two-nucleon (NN) observables and EFT will be necessary to determine at which quark masses one may trust results for 0⌫ -decay observables.
Before calculating the nn ! pp amplitudes with LQCD, however, the low-energy spectra and scattering amplitudes in the NN system need to be calculated with precision. Doing so allows one to determine which operators couple su ciently well to the ground states of interest, to understand the systematic uncertainties inherent to NN systems, and to match finite-volume Euclidean matrix elements to infinite-volume transition amplitudes. The NN studies to date have been largely carried out at very large quark masses, where extrapolation to the physical point cannot be controlled. Fully understanding the systematic uncertainties will become even more crucial as the quark masses are lowered toward their physical values because of a signal-to-noise problem for nucleons, in which statistical noise grows exponentially with the pion mass, atomic number, and Euclidean time [84][85][86]. Furthermore, in calculations at lighter quark masses (which require larger lattice volumes), the energy gaps that dictate the exponential decay of excited states with Euclidean time become very small, causing a slow approach to the ground state that may be obscured by the growth in noise. Thus, improved operators, analysis, and understanding of excited-state contamination are of critical importance.
These complications mean that we still do not know whether two nucleons form bound states, even at large quark masses that make precise calculations easier. Recent work within the LQCD community has highlighted the importance of fully-controlled calculations in the NN sector. First, the use of improved interpolating-operator sets and analysis techniques based on the variational principle of quantum mechanics has led to results [87][88][89] that cast doubt on earlier spectroscopy calculations at similar quark masses. Second, a preliminary study of the discretization e↵ects of two-baryon calculations has shown large shifts in the binding energies away from the continuum limit [90]. The latter finding, in particular, needs to be verified by di↵erent groups with di↵erent lattice actions, and may indicate that NN calculations must be performed at multiple fine lattice spacings. This would significantly increase the cost of calculations.
To use LQCD to access the 0⌫ amplitude, one needs to develop indirect mapping relations. This is because the notion of asymptotic states is absent in the finite Euclidean spacetime that is that used in the LQCD setting. A general mapping exists to obtain matrix elements of local (short-range) operators such as those appearing in the highscale models of 0⌫ decay within two-nucleon states [49]. As an input, this mapping requires two-nucleon spectra and the energy dependence of elastic scattering amplitudes. The existing mapping for the matrix element associated with light neutrino-exchange involves a matching to the leading-order nucleonic EFT [50], and requires as input the two-nucleon spectra and scattering amplitudes. The mappings for such long-range matrix elements are in general more complex than those for local matrix elements because a straightforward analytic continuation in the presence of on-shell intermediate states is not possible [91][92][93][94]. With properly infrared-regulated neutrino propagators, however, analytic continuation will be straightforward in future calculations of the 0⌫ amplitude [47,48]. Techniques for computing both the short-and long-distance contributions to 0⌫ processes have already been developed and applied in studies of the ⇡ ! ⇡ + e e and ⇡ ⇡ ! e e processes [45,47,48], which also constrain pionic contributions within nuclear 0⌫ decays. Calculations of the nn ! ppe e process will be significantly more involved for the reasons discussed above, but will be a crucial next step.
With the broad goal of achieving a systematic quantification of nuclear uncertainties, we must face the challenges of extending the analysis of the 0⌫ transition operator beyond leading order in chiral EFT, especially in the case of light-neutrino exchange. We must also go beyond two-nucleon observables to reliably determine the role of multinucleon e↵ects in double-beta decay. Regarding the nn ! pp amplitude, both in Weinberg's power counting (WPC) and in renormalized chiral EFT, the first corrections arise at next-to-next-to-leading order (N 2 LO) [95]. In the twobody sector, the transition operator includes contributions from the nucleon vector, axial, and induced pseudoscalar form factors (which are customarily included in nuclear calculations), from pion-neutrino loops [42], and from new contact interactions required to absorb the divergences in these loops. These include the couplings of two electrons to two pions (g ⇡⇡ ⌫ ), to two nucleons and one pion (g ⇡N ⌫ ), and to four nucleons (a correction to g NN ⌫ ). g ⇡⇡ ⌫ is well determined by LQCD [47,48], while extracting the correction to g NN ⌫ will require the matching of LQCD and chiral-EFT amplitudes at higher orders. The short-range structure of the two-body 0⌫ operator at N 2 LO is at the moment unknown beyond WPC. Reference [95] pointed out that the promotion of g NN ⌫ to LO implies that certain derivative operators in the spin-singlet channel are also more important than in WPC, but a full analysis of the LNV scattering amplitudes to N 2 LO does not yet exist, and needs to be developed to interpret anticipated LQCD results. A deeper question is whether the chiral and momentum expansions of chiral EFT converge (and converge to what is observed in Nature). This question is open even for single baryons [83,96] and relatively light systems [97][98][99], and needs to be answered as the community moves beyond purely phenomenological approaches. LQCD input for the unknown LECs at successively higher orders can help resolve power-counting questions for LNV processes.
Moving to the question of multi-nucleon corrections, we note that two-body currents, which are important in the g A -quenching problem in decays [100,101], first contribute to 0⌫ decay at N 2 LO, by generating three-body corrections to the operator. These corrections were considered in Refs. [102,103] and found to be compatible with power-counting estimates. Furthermore, in the three-body sector, a goal for the chiral-EFT community is to validate or falsify WPC's expectations, by studying suitable few-body amplitudes. Once calculations of two-body transitions have been achieved with systematic control, LQCD studies of 0⌫ decay of A 2 {4, 6} systems, if they can be carried out, will provide valuable additional information. Such calculations can reduce systematic uncertainties in the process of matching 0⌫ amplitudes to the chiral EFT used in nuclear many-body calculations. In particular, constraining the same LECs from LQCD calculations of di↵erent processes will not only reduce statistical uncertainty, but will also, through benchmarking, reduce the uncertainties in nuclear EFTs that arise from choices of scheme or regulator. Useful transitions will probably include the A = 4 processes 4 H ! 4 Li e e and 4 n ! 4 He e e , and the A = 6 transitions 6 He ! 6 Be e e (for which nuclear-structure calculations have been performed [104]), and 6 H ! 6 Li e e (which introduces additional challenges for many-body approaches because 6 H is unstable). To reduce the cost of extrapolating such LQCD calculations to infinite volume, directly matching matrix elements to finite-volume EFT calculations [105][106][107][108] may be a valuable approach to precisely determining the LECs.
In summary, while significant outstanding challenges must be overcome to reliably determine the nn ! pp amplitude, for both the light-neutrino exchange and the short-distance I = 2 4-quark operators, there exists a clear road map for addressing them. Following it will require a concerted e↵ort in LQCD, EFT, and the coupling of these theories, as well as computing resources at the exascale and beyond, both to quantify the uncertainties in the relevant two-body process and to build an understanding of multi-nucleon corrections.
candidate nuclei, with the exception of 48 Ca, are open-shell and at least of medium mass. Consequently, only a subset of the currently-available ab initio many-body methods can be used to compute the nuclear matrix elements that govern their decay. First computations of the nuclear matrix elements have been performed in coupled cluster theory [63], the IM-GCM [61], and the VS-IMSRG [62]. We describe each of these methods and prospects for improving them next.
In the coupled cluster method [109][110][111], the exact wave function | i is parameterized by the exponential ansatz | i = e T | 0 i, where the reference | 0 i is a product state, and the cluster operator T generates particle-hole excitations. One expresses T in terms of single, double, triple etc. particle-hole excitations and (usually) truncates it at the so-called doubles or triples level. This is the main approximation. The calculation of transition matrix elements in coupled cluster theory is complicated by the fact that T is purely an excitation operator, i.e. the fact that the similarity transform e T Ôe T of a Hermitian operator Ô is not Hermitian. This implies that the bra version of a state needs to be parameterized through a de-excitation operator rather than an excitation operator. An additional complication for 0⌫ decay is that the initial and final states are the ground states of di↵erent nuclei, with each in principle requiring its own T operator and reference state. In practice, one expresses the final state as a generalized excitation of the initial state through the equation-of-motion method as | F i = e T R| 0 i, where R is a double-charge-changing excitation operator [63]. Alternatively, one can express the initial state as an excitation of the final state. In the absence of any truncation, these two choices should yield identical results, so the di↵erence between the two is an indication of the truncation error.
Like T , the excitation operator R is expanded in terms of charge-changing few-nucleon "excitations" and truncated at a doubles or triples level. This approximation may not be accurate when the initial and final nuclei are very di↵erent in structure, because, for example, they di↵er in their intrinsic deformation. Indeed, the spherically-symmetric coupled cluster method works well for computing properties of closed-shell nuclei such as 48 Ca. However, the ground state of 48 Ti (the final state in the decay of 48 Ca) is open-shell and is better treated with an intrinsically deformed, (though axially symmetric) reference state, which is computationally more expensive. In benchmarks performed so far [63], it appears that taking | 0 i to be the deformed 48 Ti state yields more accurate results, though the reason is not entirely known.
We can expect this approach to be applied to more nuclei and with more accuracy in the next few years. With enough computation time, it can be generalized to allow triaxial deformation of the reference state (and thus a good calculation, e.g., in 76 Ge) and the restoration of rotational symmetry through projection onto states with good total angular momentum [112]. These developments will turn the method into a much more versatile tool for the computation of 0⌫ nuclear matrix elements.
The IM-GCM [61,113] is a combination of the GCM [114] and the Multi-Reference In-Medium Similarity Renormalization Group (MR-IMSRG) [115,116]. ("Multireference" refers to a generalization of the renormalization-group flow equations to work with a reference state that is more complex than a Slater determinant.) The GCM e ciently captures the collective long-range correlations which are important in deformed nuclei, while the MR-IMSRG captures short-range correlations associated with the repulsive core of realistic nuclear interactions.
On can view the IMSRG as a way to generate a unitary transformation U of the Hamiltonian that brings it to a form more amenable to solution. The transformation is parameterized by a flow parameter s, yielding a di↵erential equation for the transformed Hamiltonian H(s) = U (s)HU † (s) and for other consistently-transformed operators O(s) = U (s)OU † (s). The unitary transformation, which is conveniently expressed in the Magnus formulation as U (s) = e ⌦(s) [117], is designed so that with increasing s, a reference state | 0 i increasingly approximates an eigenstate of H(s).
In the IM-GCM, the approach is to take | 0 i to be the ground state of a GCM calculation. The GCM ground state is expressed as a linear combination of configurations | (q)i labeled by a set of generator coordinates q (e.g. quadrupole deformation), so that | 0 i = P q f (q)| (q)i. The amplitudes f (q) are obtained by minimizing the energy via the Hill-Wheeler-Gri n equation, which amounts to a diagonalization the space of GCM states | (q)i. As we noted in the context of coupled cluster theory, the initial and final states in any 0⌫ decay are di↵erent, a fact that complicates most computations. The transformations e ⌦ I (s) and e ⌦ F (s) that decouple the initial and final states are not equivalent. In Ref. [113], this complication was addressed by combining the IMSRG transformation and GCM calculations for the initial and final states in di↵erent ways, again with the understanding that, without approximations, all these combinations should give the same results. In more recent work, a powerful alternative was presented in the form of an ensemble composed of reference states in both the initial and final nuclei that allows one to use a single transformation rather than two separate ones [61].
The main approximation in the MR-IMSRG flow equations is that all operators are truncated at the normal ordered two-body (NO2B) level. In the next few years, with enough computational capacity, we will be able to go beyond this approximation by either exactly or approximately including the e↵ects of three-body operators that are induced by the flow equations. A first step in this direction indicated that the correction due to induced three body operators is sub-leading (on the order of 10% of the NO2B correction) [61,118]. The result is encouraging, but the corrections are large enough that they should be included.
Another approximation is in the selection of generator coordinates. In principle, one can continue to add more coordinates that are believed to be relevant (for example, proton-neutron pairing gaps) and confirm that the answer does not change, but it is di cult to establish that all important degrees of freedom have been explored. Historically, this has been a significant issue for the GCM calculations based on phenomenological interactions. In the IM-GCM, this issue can be overcome because dependence of the transformation on the flow parameter s o↵ers a powerful diagnostic tool: If su cient degrees of freedom are included in the MR-IMSRG flow and the GCM basis, the unitarity of the transformation will not be spoiled by truncation errors, and all observables should be independent of s [119,120].
In the VS-IMSRG [121], as in the IM-GCM, the strategy is to perform a unitary transformation to bring the Hamiltonian into a form more amenable to solution. In this case, the transformation block-diagonalizes the Hamiltonian such that an additional diagonalization in a valence shell-model space yields exact results (assuming no truncations are made in the transformation).
The VS-IMSRG calculations carried out thus far have generally performed the normal ordering with respect to a closed-shell reference state or an uncorrelated "ensemble" reference that has the correct number of protons and neutrons on average. As with coupled cluster theory, one needs to choose the initial or final state as the reference and, in the absence of truncation, this should not a↵ect the answer. The simpler reference used in the VS-IMSRG (compared to, e.g., the IM-GCM) is somewhat compensated for by the subsequent exact diagonalization in the valence space, resulting in a complementary approximation scheme. Like IM-GCM, the VS-IMSRG as currently practiced truncates operators at the two-body level after normal ordering, and the clear next step is the approximate inclusion of the e↵ects of induced three-body operators. Again, with su cient computational resources and personpower, this can be done.
Quasi-exact ab initio methods, namely quantum Monte-Carlo (QMC) and the no-core shell model (NCSM), are generally limited to light systems, which are not directly relevant for 0⌫ experimental searches. They play an important role, however, in benchmarking the methods we've discussed, which can reach the relevant heavier systems. The three methods described above, coupled cluster theory, IM-GCM, and VS-IMSRG, have all been benchmarked against the NCSM in light systems up to 14 C (and up to 22 O with the importance-truncated NCSM) [61,63,118,122]. The benchmarks showed that coupled cluster calculations that use a deformed reference state are usually more accurate than those that use a spherical reference.
In contrast to the NCSM, IM-GCM, VS-IMSR, and CC theory, QMC approaches do not rely on a single-particle basis expansion. Variational Monte Carlo (VMC) approximates the solution of the many-body problem by an accurate trial wave function T , obtained by applying two-and three-body correlation operators to a Slater determinant of A single-particle wave functions [123,124]. The optimal set of variational parameters defining the trial wave function is obtained by minimizing the energy expectation value h T |H| T i with dedicated optimization algorithms [125]. The limitations of the variational ansatz are overcome by the Green's function Monte Carlo (GFMC) method that propagates the trial wave function in imaginary-time to extract the ground-state of the system
QMC methods have no di culty in using "sti↵" forces that can generate wave functions with high-momentum components, but they are limited to local (or nearly local) Hamiltonians because nonlocalities exacerbate the fermion-sign problem [126]. There have been QMC studies of the 0⌫ -decay nuclear matrix elements for light nuclei (see, e.g., [127,128]), but the (nearly) local Hamiltonians [129][130][131][132][133] used in these studies pose a substantial hurdle for direct benchmarks against the configuration-space methods that we discussed above. More recently developed local chiral interactions with typical cuto↵s around ⇤ = 500 MeV lead to much slower convergence than their nonlocal counterparts with the same scales. Renormalization group transformations may help to mitigate this problem, but the uncertainties due to the omission of induced contributions to the interaction and transition operators may also be more substantial than in a nonlocal regularization scheme. Nevertheless, once RG and EFT truncation errors have been propagated to the 0⌫ matrix element, a comparison of QMC and configuration-space methods will represent an important check.
We have focused on the several ab initio methods that have already been applied to experimentally relevant transitions, but there are others that may also be able to tackle these nuclei soon. Applications of QMC have been limited almost entirely to the p-shell and below because the number of spin/isospin states scale exponentially with particle number A (see, e.g., [134] and references therein). However, within the auxiliary field di↵usion Monte Carlo (AFDMC) approach [135] the spin-isospin degrees of freedom are described by single-particle spinors, the amplitudes of which are sampled with Monte Carlo techniques based on the Hubbard-Stratonovich transformation. The transformation reduces the computational scaling from exponential to polynomial in A. AFDMC calculations for 16 O have been reported [131], and calculations of 48 Ca are conceivable in the near future.
A recently proposed alternative is to use QMC to compute M 0⌫ for light nuclei, and match an e↵ective shell-model operator to these calculations [60], using the generalized contact formalism (GCF). The e↵ective operator is then employed in shell-model calculations of heavier nuclei, where QMC is not feasible. This approach can be viewed as the QMC providing synthetic data to which a shell-model e↵ective operator can be fit. It is justified by the factorization of physics at the scale of nucleon-nucleon interactions from the nuclear environment. This factorization is seen in the application of renormalization group (RG) methods to QMC wave functions, where short-distance physics in those wave functions evolves into e↵ective operators at the lower resolution appropriate to the shell model [136]. The GCF implements the leading-order consequences of factorization. One challenge for the future will be quantifying the long-range correlations missed by the shell model; these will in general depend on the valence space (see e↵ective charges for E2 transitions as an example). Such quantification is one aspect of a more general question about what the sub-leading corrections to the GCF calculation carried out in Ref. [60] are.
The RG approach to this problem makes it clear that, in any of the approaches to the nuclear many-body problem described here, the 0⌫ operator must be evolved consistently with the methods used to reduce the e↵ective size of the space in which the many-body wave functions are computed. It follows that the 0⌫ contact operator will not necessarily be the same as the one computed in Ref. [78], or obtained in the future from LQCD. Instead that short-distance operator must absorb the physics between the hadronic scale of LQCD/sum-rule calculations and the low-energy nuclear-structure scale; i.e., it will account not just for hadronic excitations that have been integrated out of the Hilbert space, but for high-energy nuclear correlations that are integrated out too.
The NCSM may also be applied to heavier istopes in the future. Although the method in its original form is typically limited to A . 16, the importance-truncated NCSM (IT-NCSM) [137] can significantly reduce the dimensions of the Hamiltonian and thus reach higher in mass, conceivably up to 48 Ca. One challenge will be to obtain a better understanding of the extrapolation of M 0⌫ in the importance truncation parameter min .
The symmetry-adapted no-core shell model (SA-NCSM) is a version of the NCSM that uses irreducible representations of the symplectic symmetry group Sp(3, R) rather than particle-hole energy to truncate its basis [138][139][140][141]. The alternative truncation scheme allows it to e ciently capture deformation, which is important in either the mother or daughter nucleus in all experimentally relevant 0⌫ -decay candidates. Applications of the SA-NCSM have mostly focused on p-and sd shell nuclei so far, but first results for 48 Ti have been reported in Ref. [140]. These particular results serve as a demonstration that convergence of a SA-NCSM calculation is mainly a↵ected by the strength of the mixing between irreps in a particular nucleus rather than the mass number: the model space dimension for 48 Ti is more than an order of magnitude smaller than the dimension for 20 Ne. At present three-nucleon forces have not yet been included in the SA-NCSM, but once this challenge is overcome, the SA-NCSM will be a valuable complementary approach to the coupled cluster and VS-IMSRG methods that employ particle-hole based truncations. It is also complementary to the IM-GCM because the Sp(3, R) irreps o↵er a more systematic approach to basis construction than the selection of relevant generator coordinates.
Both the conventional NCSM and SA-NCSM can also be combined with an (MR-)IMSRG preprocessing of the Hamiltonian and transition operators to accelerate convergence [119]. A combination of MR-IMSRG evolution and SA-NCSM, in particular, would embrace a similar philosophy as does the IM-GCM.
Besides the ab initio methods described in the previous subsections, a variety of other methods have been used to compute M 0⌫ in nuclei of interest to experimentalists. These others can be broadly grouped into categories: the interacting shell model [142][143][144][145][146][147][148][149][150], energy-density-functional (EDF) methods [151], the quasi-particle random phase approximation (QRPA) [54,[152][153][154][155], and the interacting boson model (IBM) [156,157].
In this report the emphasis is on methods that can, in principle, quantify the theoretical uncertainties of the underlying strong-interaction Hamiltonian and of the transition operators. However, methods such as the "phenomenological" shell model still have a valuable role to play, because they preserve underlying nuclear many-body symmetries and thus capture the most relevant degrees of freedom. In addition, the heavy work of finding optimized e↵ective shellmodel Hamiltonians and e↵ective transition operators that capture the landscape of realistic nuclear spectra and of the experimentally accessible nuclear observables has already been done in these approaches. We envision that semi-phenomenological methods, such as the shell model, can be used to explore correlations between observables, helping us identify the quantities that best reflect the accuracy of an ab initio calculation. For example, an ensemble of shell model Hamiltonians can be generated by adding random contributions to the two-body matrix elements of some "seed" Hamiltonians [158]. These Hamiltonians can then be used to obtain M 0⌫ , as well as excitation spectra and electroweak transitions and moments (for which data exist or could be obtained). Any observables which are significantly correlated with M 0⌫ would then be explored in the more expensive ab initio calculations, producing input for subsequent model-mixing analysis. An initial study along these lines for the 0⌫ decay of 48 Ca-48 Sc- 48 Ti system can be found in Ref. [159].
The preceding sections outlined a variety of many-body methods that can be used to perform ab initio calculations of 0⌫ nuclear matrix elements. It might be supposed that the goal of a UQ analysis should be to determine the "best" of these methods and that whichever method turns out to be "best" should then be used exclusively. In fact, these methods have complementary strengths and deficiencies, so the goal instead is to use all of them to optimize the overall predictions. Therefore in this section we outline a procedure, depicted in Fig. 3, by which the results for M 0⌫ obtained in those di↵erent methods-as well as their uncertainties-can be combined into a single, unified prediction for M 0⌫ . We also explain how that procedure will naturally suggest alternative strategies for calibration of the ab initio calculations which should, in turn, lead to refined predictions for M 0⌫ .
Throughout this section we have in mind that we are considering ab initio predictions for M 0⌫ that are obtained with chiral-EFT forces and decay operators. Di↵erences between di↵erent implementations of the chiral-EFT force should therefore be encompassed within the uncertainty assigned due to truncation of the EFT expansion, cf. Sec. II A above. The source of uncertainty that is hardest to assess is therefore that due to the use of di↵erent methods for solving the A-body problem: these are associated with di↵erent ways of truncating the A-body Hilbert space. In what follows we denote the di↵erent many-body methods that have been, or may in the future be, adopted to solve this problem as M k . We treat these as di↵erent "models" in the statistical sense of the term "model" and seek to combine their predictions into a single prediction that accurately assesses uncertainties in the evaluation of M 0⌫ . We assume that uncertainties due to the truncation of the chiral-EFT expansion are reflected in the posterior distribution that must be provided by each many-body method, M k .
Method M k 's prediction also has an inherent parametric uncertainty, coming both from the parameters of the Hamiltonian used to obtain the wave function of the initial and final state in the double-beta-decay process, and from the contact piece of the 0⌫
operator. In what follows we denote the low-energy constant that multiplies the contact piece by ⌘ and the parameters of the Hamiltonian as ✓.
Recently, ⌘ has been determined [161] by reproducing the synthetic datum, y synth provided in Refs. [43,44]. Meanwhile, for most of this section we will assume that the parameters ✓ are calibrated to a dataset y (cf. Sec. II B) that does not have to include observables that we expect are correlated with 0⌫ decay. This is, after all, the stated orientation of most ab initio approaches, which calibrate the parameters of NN and three-nucleon forces (3NFs) to NN scattering data and a few observables in light nuclei. The posterior probability distribution for the parameters ✓ that is obtained from such an analysis is denoted p(✓|y).
From a Bayesian perspective, each many-body method's prediction of M 0⌫ also comes with a systematic error that depends on parameters of the approach employed, e.g., Hilbert-space size, accuracy of treatment of 3NFs, etc. The statistical modeling of this systematic error is referred to as discrepancy learning (cf. Eq. ( 2)). Simultaneous learning of discrepancy and parameters is a complicated practical and theoretical exercise. Moreover, with no information to leverage near the quantity of interest M 0⌫ , verification of discrepancy can be di cult. Nonetheless, grounded, informed priors on the discrepancy can improve prediction-especially when we seek to leverage the predictions made FIG. 3. The road to calculations of M0⌫ with UQ that accounts for all limitations of the nuclear-physics calculation: truncation errors in chiral EFT, uncertainties in the theory's parameters, and deficiencies of the many-body methods used. Emulator development is the first step, as it is key to facilitating subsequent calculations. Model parameters ✓ can be calibrated against a dataset y. Weights for the di↵erent many-body methods will be obtained by assessing methods' performance on a set of observables {yev}. The weights w k (yev) are to be computed via scoring rules that gauge predictive accuracy [160] and will also take account of the extent to which di↵erent members of {yev} are correlated with M0⌫ .
across several many-body methods. We therefore write:
where M M 0⌫ (✓, ⌘; ) is the prediction obtained in method M with method hyperparameters (and at specific Hamiltonian and operator parameter values) and M M 0⌫ ( ) is the corresponding model uncertainty. If we, for the moment, ignore the issue of the discrepancy function, then the method-M prediction of M 0⌫ is formed by marginalizing over ✓ and ⌘ using the distributions established for them from the data {y, y synth }:
Here it should be noted that we have allowed for the possibility that the probability density obtained for both the Hamiltonian parameters and ⌘ is di↵erent for di↵erent methods, i.e., depends on the method M. We have, however, assumed that all all methods are calibrated using a common data set y.
But the problem of model discrepancy is critical in predictions for neutrinoless double-beta decay: di↵erent approaches to the nuclear many-body problem are based on di↵erent physics assumptions, and so have di↵erent model discrepancies. In order to get a handle on the model discrepancy we propose to assess that method's ability to predict observables that may share similar physics features to the 0⌫ decay in the nucleus of interest. Candidate processes include:
• Single -decay rates in neighboring nuclei, e.g., in the intermediate nucleus in 0⌫ decay;
• -strength distributions;
• Known 2⌫ decay rates;
• Magnetic moments and B(M 1) rates in three nuclei involved in 0⌫ decay;
• Energies of the lowest J ⇡ = 2 + states and B(E2, 2 + ! 0 + ) rates in initial and final nuclei; • Charge radii;
• Observables probing a 100 MeV momentum-transfer scale, e.g., in muon capture.
The idea is then that a method that performs well on these observables, which we denote collectively as y ev , should be a more accurate predictor of M 0⌫ than one that does not.
But which of these observables are most important for constraining the 0⌫ decay rates? Until now discussion on this point has been largely driven by qualitiative arguments. We propose that, by using properly calibrated Hamiltonians, this question can be answered by analyzing the correlations between the observables on the list above and M 0⌫ . Those correlations can be well approximated by drawing a finite number of samples from p(✓|y) (say ⇡ 100), using model M to compute each observable in the set y ev and M 0⌫ , and extracting the empirical correlation coe cient of M 0⌫ and each quantity in y ev for that model M. Examples of such sensitivity studies can be found in, e.g., Refs. [162][163][164].
A few supplementary points regarding this correlation analysis need to be made:
• This correlation need not be the same in every method. Di↵erent methods have di↵erent discrepancy functions, because di↵erent methods truncate the nuclear many-body problem in di↵erent ways. A particular example of this is that methods with di↵erent resolution scales may di↵er in whether their discrepancy function reflects errors in the long-distance physics or errors in the short-distance physics. Indeed, the balance between these two types of errors could shift within a particular method as the value of the hyperparameters changes. It follows that the correlation found for method M at particular values of may depend on either or M. But, analysis of these correlations, when combined with data on the observables in the set y ev , will help us pin down the discrepancy functions, or at least minimize their impact on the M 0⌫ prediction.
• The prediction of the observables y ev may depend on additional parameters , that are not part of the set ✓, and are not a priori needed to predict M 0⌫ . The posterior prediction for y ev is then formed by marginalizing over , using a probability distribution function that can be thought of as a prior for our purposes, but may be informed by studies of the pertinent observable(s) in nuclei that are some distance from the 0⌫ candidate p(y ev |y, M) = Z p(y ev , ✓, , M)p(✓|y)p( |M)d✓d ,
Marginalizing over ⌘ will certainly be necessary for the prediction of M 0⌫ . Such marginalization (over ⌘ and ) may a↵ect the correlation.
• As discussed in Sec. III B 6, the correlation analysis does not need to be initiated within computationallyexpensive ab initio calculations. In the first instance, it can be carried out using lower resolution approaches that can be viewed as computationally less expensive emulators of ab initio methods. An example of ongoing work along these lines can be found in Ref. [159].
Such a correlation analysis is useful in its own right. But, with these correlation coe cients in hand we can form an "improved Bayesian Model Average" of the results from the di↵erent many-body methods M. In BMA [75], a set of candidate methods M 1 , . . . , M p are combined to form a predictive distribution via Eq. (3):
where w k (y ev ) represents the underlying weight given to each method. The weight w k is proportional to p(y ev |M k ) ⇥ p(M k ), where p(y ev |M k ) represents the evidence for method M k present in data y ev and p(M k ) is a prior probability that method k is correct. This prior is usually taken to be flat across the methods M k , i.e., the methods are all taken to be equally plausible. The formula ( 7) is aspirational in that there are complications in its deployment that are both practical and theoretical.
The selection of the weights requires careful documentation of the source of the systematic errors between method predictions and observables. Extraneous observables, i.e., observables that are no more than weakly related to 0⌫ , are less dangerous to the resulting inference if all error sources were only experimental and independent throughout the observable space as no method has a specific advantageous bias. However, in the context of 0⌫ one expects that a significant portion of the error can be attributed to systematic method deficiencies y th . The result is significantly related error that exists across the observable space.
It is therefore critical to carefully select observables in y ev . Observables that are closely related to the target observable, y ⇤ = M 0⌫ should receive higher weights-something that classical BMA does not do. If observables in y ev that are not singificantly related to M 0⌫ can influence the weight of a method in the BMA formula (7), they are likely to dilute-or even bias-the prediction for M 0⌫ . Parsimony is also important for a practical reason: the Bayesian model averaging formula requires the same y ev should be used across all many-body methods, meaning all of them need to be able to produce predictions of these quantities. A parsimonious choice of observables makes it more feasible that y ev can be predicted in all candidate approaches to the nuclear many-body problem.
The program for UQ that we have laid out up to this point in the section could be carried out using ab initio methods and already calibrated chiral-EFT forces and operators. We now discuss a longer-term strategy for refining the prediction of M 0⌫ . Once it has been established which observables in y ev are strongly correlated with M 0⌫ the predictive power of the methods M k can be improved by including those members of the dataset y ev in the dataset used to calibrate the nuclear Hamiltonian and decay operators. The parameter vector ✓ would then be readjusted within each calculation M k . This would open a window for a more refined combined prediction of the di↵erent many-body methods.
One critical challenge to overcome is providing predictions of quantities for various values of ✓, ⌘, and for all methods. Part of this is needed for the integration to compute p(y|M k ), where MCMC must be leveraged for both the core Hamiltonian parameters ✓ as well as the ancillary parameters and the 0⌫ parameter ⌘. This problem is amplified in the case of high-dimensional parameter space.
There are a few statistical and computational tools that can be deployed to resolve this problem. Firstly, reducing the space of parameters through screening will be critical for each method, thereby including only the parameters that are critical for predicting y ev and M 0⌫ . Secondly, emulators (or surrogates) can play a vital role to conduct Bayesian inference from only a few full ab initio evaluations of the matrix element. Emulators can take on a variety of forms but the function is shared: to provide a computationally inexpensive approximation of y ev and M 0⌫ for any value of ✓, ⌘, and . Gaussian process-type emulators exploit smoothness in the computation's response to the parameters. Other reduced-basis emulators instead build modified, cheaper alternatives to the full ab initio result when it is subject to specific structure [165][166][167][168][169][170]. Variational machine learning methods [171][172][173] form another path that produce e cient emulators that intrinsically learn the optimal latent parameter space needed for robust interpolation and prediction. These methods have an advantage over many forms of emulation as they can learn highly non-linear manifolds while still generating a notion of the emulators' internal uncertainty.
One common theme across all methods is they rely on using specific computations at designated parameter combinations to build predictions on other ones. The amount of computations needed to build an adequate emulator varies across methodologies. Classical literature on Gaussian process-type emulators suggest computing at ten times the dimension of the parameters, but that suggestion has recently been reconsidered and higher-dimensional parameter spaces perhaps need more computation. Good emulators will naturally come with their own uncertainty quantification, which is critical for producing valid approximations of p(y ev |M k ) as well as predictions of p(M 0⌫ |M k ).
As mentioned above, lower resolution approaches, such as DFT, can be used to construct emulators of higher resolution ab initio methods. This construction can follow Eq. ( 2) in which y exp (x; ✓) is replaced by ab initio predictions and the model discrepancy y th would model the di↵erence between predictions of ab initio and lower resolution models.
Accurate calculation of the nuclear matrix elements governing neutrinoless double-beta decay, with quantified uncertainty, is essential for the success of the impressive experimental and theoretical worldwide e↵ort in this area [37,174]. The purpose of this report is to lay out the challenges to the nuclear-theory community-in regard to both nuclear-physics calculations and uncertainty quantification therein-and map out a path for near-and long-term progress.
It will require a concerted e↵ort in both LQCD and EFT, as well as the coupling of these theories, to fully quantify the theoretical uncertainties related to the 0⌫ -decay operator and associated matrix elements at the hadronic level. Systematic improvements in nuclear many-body methods are underway, and should be ready to produce a new generation of M 0⌫ matrix elements in the next few years. The complexity of the problem makes the use of multiple many-body methods well worthwhile. Their complementary strengths and deficiencies can not only be exploited for validation but also combined through Bayesian methods to yield better overall predictions.
The implementation and application of both QCD/chiral EFT and nuclear many-body approaches will require exascale computing resources and beyond. It will also require an investment in personnel to advance EFT calculations of 0⌫ operators and to develop and maintain the codes that implement new theoretical methods. Without such an investment computing time cannot be used e ciently. One focus of future work in this area will be ensuring optimal use of the heterogeneous architectures that characterize leadership-class computers.
In addition, the cohesion within and among research groups needs to be strengthened, e.g., through the establishment of joint project resources and inter-institutional collaborations in both pure and computational theory. The first wave of ab initio nuclear matrix elements came from the combined e↵orts of multiple groups, each with its own computing resources, and were often obtained at the expense of other groups. Wait times were long, not only because of a lack of computing resources but also because certain implementation steps required the work of a single specific (and busy) person, or because large amounts of data needed to be moved between computing systems. The community needs structures that reduce the severity of these kinds of bottlenecks.
A crucial feature of this report is the emphasis it places on UQ. Without principled UQ, the usefulness of predicted M 0⌫ values for guiding experimental e↵orts, interpreting measurements, and assessing new physics will be limited. At present, few physicists working on the problem of M 0⌫ make informed choices about UQ, understand the modern UQ glossary, or consider UQ to be an essential part of "the answer". This situation can be improved through coherent inter-disciplinary collaboration of nuclear physicists with applied mathematicians, statisticians, and computer-science experts.
Such a collaboration could carry out the concrete, multi-staged, and interwoven program of nuclear-physics and UQ methodological improvements and computations laid out in this report. In concert with continued strong support for the e↵orts of PIs and research groups working on 0⌫ decay, this will make the ultimate goal of accurate and precise M 0⌫ predictions achievable.
•