The isotensor axial polarisability and lattice QCD input for nuclear double-$\beta$ decay phenomenology

TL;DR

This paper uses lattice QCD and effective field theory to evaluate short-distance nuclear effects in double-beta decay, highlighting the importance of the isotensor axial polarisability for accurate nuclear modeling.

Contribution

It introduces the first lattice QCD calculation of the isotensor axial polarisability relevant for double-beta decay, emphasizing its significance in nuclear models.

Findings

Short-distance effects are comparable to axial charge quenching.

Nuclear models should include isotensor axial polarisability effects.

Lattice QCD can help constrain nuclear double-beta decay parameters.

Abstract

The potential importance of short-distance nuclear effects in double-$\beta$ decay is assessed using a lattice QCD calculation of the $nn\rightarrow pp$ transition and effective field theory methods. At the unphysical quark masses used in the numerical computation, these effects, encoded in the isotensor axial polarisability, are found to be of similar magnitude to the nuclear modification of the single axial current, which phenomenologically is the quenching of the axial charge used in nuclear many-body calculations. This finding suggests that nuclear models for neutrinoful and neutrinoless double-$\beta$ decays should incorporate this previously neglected contribution if they are to provide reliable guidance for next-generation neutrinoless double-$\beta$ decay searches. The prospects of constraining the isotensor axial polarisabilities of nuclei using lattice QCD input into nuclear many-body calculations are discussed.

Figures

Equations

$$[T_{1/2}^{2\nu }{]}^{-1}$$

$$M_{GT}^{2\nu }$$

$$⟨pp|J_3^{+}(x)J_3^{+}(y)|nn⟩$$

$${\mathcal{L}}^{(1)}$$

$${\mathcal{L}}^{(2)}$$

$$M_{nn\to pp}$$

$$C_{{\lambda }_u;{\lambda }_d}^{(h)}(t)$$

$$C_{{\lambda }_u;{\lambda }_d=0}^{(h)}(t)$$

$$C(t)$$

$$\mathcal{R}(t)$$

$$a^2\hat{\mathcal{R}}(t)$$

$${\mathcal{R}}^{\text{(lin)}}(t)=\frac{(e^{a\Delta }+1)\hat{\mathcal{R}}(t+a)-\hat{\mathcal{R}}(t+2a)-e^{a\Delta }\hat{\mathcal{R}}(t)}{e^{a\Delta }-1}$$

$$\displaystyle{\mathcal{R}}^{\text{(full)}}(t)={\mathcal{R}}^{\text{(lin)}}(t)-\frac{|\langle pp|\tilde{J}_{3}^{+}|d\rangle|^{2}}{a\Delta}\stackrel{{\scriptstyle t\to\infty}}{{\longrightarrow}}{M^{2\nu}_{GT}\over 6\,aZ_{A}^{2}}.\ \ \$$

$$\frac{\Delta }{g_A^2}\sum _{{\mathfrak{l}}^{\prime }\ne d}\frac{⟨pp|{\tilde{J}}_3^{+}|{\mathfrak{l}}^{\prime }⟩⟨{\mathfrak{l}}^{\prime }|{\tilde{J}}_3^{+}|nn⟩}{E_{{\mathfrak{l}}^{\prime }}-E_{nn}}$$

$$\frac{1}{6}\frac{\Delta }{g_A^2}{M}_{GT}^{2\nu }$$

Full text

††preprint: INT-PUB-16-056, MIT-CTP-4867

Double-$\beta$ ($\beta\beta$) decays of nuclei are of significant phenomenological interest; they probe fundamental symmetries of nature and admit both tests of the Standard Model (SM) and investigations of physics beyond it Mohapatra et al. (2007). Consequently, these decays are the subject of intense experimental study, and next-generation $\beta\beta$-decay experiments are currently being planned Avignone et al. (2008); Bilenky and Giunti (2012); Dell’Oro et al. (2016). At present, both the robust prediction of the efficacy of different detector materials, necessary for optimal design sensitivity, and the robust interpretation of the highly sought-after neutrinoless $\beta\beta$-decay ($0\nu\beta\beta$) mode are impeded by the lack of knowledge of second-order weak-interaction nuclear matrix elements. These quantities bear uncertainties from nuclear modelling that are both significant and difficult to quantify Engel and Mene ndez (2016). Controlling the nuclear uncertainties in $\beta\beta$-decay matrix elements by connecting the nuclear many-body methods to the underlying parameters of the SM is a critical task for nuclear theory.

In this Letter, lattice QCD and pionless effective field theory (EFT(${\pi\hskip-5.50003pt/}$)) are used to investigate the strong-interaction uncertainties in the second-order weak transition of the two-nucleon system in the SM by determining the threshold transition matrix element for $nn\to pp$. This matrix element receives long-distance contributions from the deuteron intermediate state whose size is governed by the squared magnitude of the $\langle pp|\tilde{J}_{\mu}^{+}|d\rangle$ matrix element of the axial current that has been recently calculated using lattice quantum chromodynamics (LQCD) Savage et al. (2017). In that work, the two-body contribution to the matrix element (i.e., that beyond the coupling of the axial current to a single nucleon) was constrained, quantifying the effective modification (quenching) of the axial charge of the nucleon from two-body effects. Here, it is highlighted that the $nn\to pp$ matrix element receives additional short-distance contributions beyond those in $|\langle pp|\tilde{J}_{\mu}^{+}|d\rangle|^{2}$ arising from the two axial currents being separated by $r<\Lambda^{-1}\sim m_{\pi}^{-1}$ (where $\Lambda$ is the cutoff scale of EFT(${\pi\hskip-5.50003pt/}$)), referred to herein as the isotensor axial polarisability. Using EFT(${\pi\hskip-5.50003pt/}$) to analyse both the second-order weak transition calculated here for the first time and the first-order $\langle pp|\tilde{J}_{\mu}^{+}|d\rangle$ transition Savage et al. (2017), the short- and long-distance contributions to the $nn\to pp$ matrix element are separately determined. Interestingly, the short-distance contribution to the total matrix element from the axial polarisability is found to be of comparable size (within the uncertainties of the present calculation) to the two-nucleon current contribution to $|\langle pp|\tilde{J}_{\mu}^{+}|d\rangle|^{2}$. Described in phenomenological terms, the polarisability thus appears to be as important as the effective quenching of the axial charge in the two-nucleon system.

The numerical calculations in this work are performed at unphysical values of the quark masses and for a disallowed decay. While there is no immediate phenomenological impact of the numerical values of the matrix elements that are extracted, the observed behavior does provide an important lesson for many-body calculations. In typical calculations of two-neutrino ($2\nu\beta\beta$) decay, the nuclear matrix elements are calculated using two insertions of the axial current in a truncated model space, with a quenched value of $g_{A}$ tuned to reproduce experiment. If the findings presented here persist at the physical values of the quark masses, they would imply that a potentially significant contribution has been ignored in standard $2\nu\beta\beta$ calculations, resulting in a source of uncertainty in the nuclear matrix elements that remains to be quantified. Importantly, this uncertainty can only be constrained using $\beta\beta$-decay measurements or numerical calculations. In $0\nu\beta\beta$ decays, the situation becomes even less certain, in part due to dependence on possible scenarios of physics beyond the SM. With a light Majorana neutrino, generalisations of the axial polarisability will also likely be relevant.

In what follows, the lattice QCD and EFT(${\pi\hskip-5.50003pt/}$) calculations and the analysis of the axial polarisability are summarised, with complete details presented in a subsequent paper Tiburzi et al. (2017). The potential for future lattice QCD calculations to provide the necessary input to constrain many-body calculations of $2\nu\beta\beta$ and $0\nu\beta\beta$ matrix elements, and thereby reduce the uncertainties in calculated $\beta\beta$-decay rates, is also discussed.

Two-neutrino $\beta\beta$-Decay: The focus of this Letter is on $2\nu\beta\beta$ decay of the dinucleon system. The decay width is given by

[TABLE]

where $Q=E_{nn}-E_{pp}$, $G_{2\nu}(Q)$ is a known phase-space factor Kotila and Iachello (2012); Stoica and Mirea (2013), and the Gamow-Teller matrix element in the two-nucleon system is

[TABLE]

Here, $J_{3}^{+}=(J_{3}^{1}+iJ_{3}^{2})/\sqrt{2}$ is the $3$rd-component of the $\Delta I_{3}=1$ axial-vector current, $J_{\mu}^{a}(x)=\overline{q}(x)\gamma_{\mu}\gamma_{5}\frac{\tau^{a}}{2}q(x)$, and ${\mathfrak{l}^{\prime}}$ indexes a complete set of zero-momentum hadronic states with the quantum numbers of the deuteron. The factors of 6 in Eq. (2) are due to rotational symmetry and our normalization of the currents. We employ $\tilde{J}^{+}_{3}=\int d\bm{x}\,J_{3}^{+}(\bm{x},t=0)$ to denote the zero-momentum current at $t=0$.

As with forward Compton scattering, the amplitude can be written in terms of a Born term, corresponding to an intermediate deuteron state, and the isotensor axial polarisability which absorbs the contributions from the remaining states in the above summation. By isospin symmetry, this polarisability is most cleanly identified as the forward matrix element of the $I=2$, $I_{3}=0$ component of the time-ordered product of two axial-vector currents in the ${}^{1}\hskip-2.168ptS_{0}$ $np$ ground-state with the deuteron pole (the Born term) omitted. For use below, isospin relations allow this matrix element to be written as

[TABLE]

where $J_{3}^{(q)}(x)=\overline{q}(x)\gamma_{3}\gamma_{5}q(x)$.

Pionless effective field theory: EFT(${\pi\hskip-5.50003pt/}$) Kaplan (1997); Kaplan et al. (1998a, b); van Kolck (1999); Chen et al. (1999); Beane and Savage (2001) efficiently describes two-nucleon systems in the regime where momenta are small compared to the pion mass. This is an appropriate tool with which to address $2\nu\beta\beta$ decays at heavier quark masses, but the inclusion of explicit pion degrees of freedom will likely be required at the physical quark masses ($0\nu\beta\beta$ decay probes higher momenta, $k\sim 100$ MeV, in large nuclei and likely also requires an EFT with explicit pion degrees of freedom). In what follows, the dibaryon formalism of EFT(${\pi\hskip-5.50003pt/}$) is utilised, using the conventions for the strong-interaction sector described in Ref. Beane and Savage (2001). The nucleon degrees of freedom are encoded in the field $N$, and the two-nucleon degrees of freedom enter as the isosinglet, $t_{i}$, and isotriplet, $s_{a}$, dibaryon fields while $y_{t}$ and $y_{s}$ describe the couplings between two nucleons and the corresponding dibaryon fields. In this formalism, the single axial-current interactions enter through the Lagrangian Butler and Chen (2000); Butler et al. (2001); Kong and Ravndal (2001); Butler and Chen (2001)

[TABLE]

where $r_{s(t)}$ is the effective range in the ${}^{1}\hskip-2.168ptS_{0}(^{3}\hskip-1.8063ptS_{1})$ two-nucleon channel, $\sigma_{i}(\tau^{a})$ are Pauli matrices in spin(flavour) space, $g_{A}$ and $\tilde{l}_{1,A}$ are the one- and two-nucleon axial couplings, and $W_{3}^{a}$ is an axial isovector field aligned in the $j=3$ spatial direction. The second term is constructed so that $\tilde{l}_{1,A}$ corresponds to a purely two-body current effect. The second-order isotensor axial interaction in the ${}^{1}\hskip-2.168ptS_{0}$ channel enters as

[TABLE]

where ${\mathcal{W}}^{ab}=W_{3}^{\{a}W_{3}^{b\}}$ is the traceless symmetric combination of two axial fields at the same location, $\tilde{h}_{2,S}$ is the scalar isotensor weak two-nucleon coupling and $\gamma_{s}=\sqrt{MB_{nn}}$ with the binding energy of the ${}^{1}\hskip-2.168ptS_{0}$ system being $B_{nn}$ (at the unphysical masses used herein, the ${}^{1}\hskip-2.168ptS_{0}$ system is bound Beane et al. (2013a)).

Calculation of the $nn\to pp$ amplitude, presented in detail in Ref. Tiburzi et al. (2017), shows that

[TABLE]

where $\Delta=E_{nn}-E_{d}$ is the difference of the ground-state ${}^{1}\hskip-2.168ptS_{0}$ and ${}^{3}\hskip-1.8063ptS_{1}$ energies. The parameter $\mathbb{H}_{2,S}$ is the short-distance two-nucleon, two-axial current coupling $\tilde{h}_{2,S}$ of Eq. (5), redefined and rescaled to capture the effects beyond the deuteron pole and two-nucleon states at energies below $\Lambda$ Tiburzi et al. (2017). This expression depends on both the long-distance contribution from the deuteron pole (the first term) and the short-distance contributions encapsulated in the second and third terms (where “short-distance” here means the “nonadiabatic” contribution; that is, every process other than those proceeding via the bound ${}^{3}\hskip-1.8063ptS_{1}$ ground state). The deuteron-pole contribution includes the effective quenching of the axial charge through $\tilde{l}_{1,A}$. A determination of the $nn\to pp$ transition matrix element, along with an extraction of the $pp\to d$ amplitude, allows for the isolation of the unknown short-distance contribution, $\mathbb{H}_{2,S}$. Once this counterterm is determined, few-body methods based on EFT(${\pi\hskip-5.50003pt/}$) or matched to them (for example, see Ref. Barnea et al. (2015)) can incorporate the axial polarisability in computations of decay rates of larger nuclei.

Lattice QCD calculations: The present lattice QCD calculations extend those of the $pp$-fusion cross section and tritium $\beta$-decay in Ref. Savage et al. (2017). The same hadronic correlators calculated in the presence of external axial fields are analysed further to access second-order weak responses to the external field. Recent calculations by the RBC/UKQCD collaboration in the kaon sector Christ et al. (2013, 2015, 2016a, 2016b) have demonstrated that long-distance second-order weak effects can be constrained using lattice QCD. These methods are extended to determine the second-order weak matrix elements of the two-nucleon system.

As discussed in Ref. Savage et al. (2017), calculations are performed on one ensemble of gauge-field configurations generated using a Lüscher-Weisz gauge action Lüscher and Weisz (1985) and a clover-improved fermion action Sheikholeslami and Wohlert (1985) with $N_{f}=3$ degenerate flavours of quarks. The quark masses are tuned to the physical strange-quark mass, producing a pion of mass $m_{\pi}\approx 806~{}{\rm MeV}$. The ensemble has a spacetime volume of $L^{3}\times T=32^{3}\times 48$ and a gauge coupling that corresponds to a lattice spacing of $a\sim 0.145~{}{\rm fm}$. For these calculations, 437 configurations spaced by 10 hybrid Monte Carlo trajectories are used and seven different sets of compound propagators are generated from sixteen smeared sources on each configuration with both smeared (SS) and point (SP) sinks. The compound propagators are produced with a single insertion of $J^{(u,d)}_{3}$ with couplings $\lambda_{u,d}=\{0,\pm 0.05,\pm 0.1,\pm 0.2\}$ . These propagators are used to produce correlation functions

[TABLE]

for all the allowed spin states of the one- and two-nucleon systems, $h\in\{p$, $np(^{3}\hskip-1.8063ptS_{1})$, $nn$, $np(^{1}\hskip-2.168ptS_{0})$, $pp\}$. Results for all source locations on each configuration are averaged before subsequent analysis.

The calculations use a lattice axial current with the finite renormalisation factor $Z_{A}=0.867(43)$ Savage et al. (2017). Because of the isotensor nature of the bilinear insertions, mixing with other operator structures is highly suppressed. As the dinucleon and deuteron states are both compact bound states at this value of the quark masses Beane et al. (2013a), only exponentially small volume effects are anticipated in the extracted matrix element. This will become a more subtle issue for future calculations with quark masses near the physical values, as discussed in Ref. Tiburzi et al. (2017).

The second-order axial responses of the dinucleon system are the primary focus of the current work. For an up-quark axial current, the relevant background-field correlators have the form

[TABLE]

from which the second-order term in the field strength, $\lambda_{u}$, can be extracted from determinations at multiple values of $\lambda_{u}$. Combining these correlators, and those for the down-quark axial current, as specified in Eq. (3), leads to the isotensor matrix element

[TABLE]

where each term on the right-hand-side is the component of the correlation functions that is second-order in $\lambda_{u}$ or $\lambda_{d}$, as denoted. Using Eq. (3) and isospin symmetry, Eq. (LABEL:eq:Cnnpp) can be written as

[TABLE]

Up to discretisation effects, insertion of appropriate complete sets of states allows this expression to be written as

[TABLE]

where $|\mathfrak{n}\rangle$, $|\mathfrak{m}\rangle$ and $|\mathfrak{l}^{\prime}\rangle$ are zero-momentum energy eigenstates with the quantum numbers of the $pp$, $nn$ and deuteron systems, respectively. Here $Z_{\mathfrak{n}}=\sqrt{V}\langle 0|\chi_{pp}|\mathfrak{n}\rangle$ and $Z_{\mathfrak{m}}=\sqrt{V}\langle 0|\chi_{nn}|\mathfrak{m}\rangle$ are overlap factors, and $E_{\mathfrak{l}^{\prime}}=E_{nn}+\delta_{\mathfrak{l}^{\prime}}$ and $E_{\mathfrak{n}}=E_{nn}+\delta_{\mathfrak{n}}$ are the energies of the $\mathfrak{l}^{\prime}$th and $\mathfrak{n}$th excited states in the ${}^{3}\hskip-1.8063ptS_{1}$ and ${}^{1}\hskip-2.168ptS_{0}$ channels, respectively.

Forming a ratio of Eq. (Isotensor axial polarisability and lattice QCD input for nuclear double-$\beta$ decay phenomenology) to the zero-field two-point function,

[TABLE]

it is straightforward, utilizing the isospin symmetry of the calculation, to show that Tiburzi et al. (2017)

[TABLE]

where $c$ and $d$ involve complicated combinations of ground- and excited-state transition amplitudes, and $\hat{\delta}$ is the minimum energy gap to the first excited state in either channel; and, for these calculations, $\hat{\delta}\gg\Delta$. Importantly, the coefficient of the linear term determines the axial polarisability and can be extracted from

[TABLE]

at late times. Finally, this result can be combined with the deuteron-pole contribution to give a quantity that asymptotes to the bare Gamow-Teller matrix element at late times,

[TABLE]

The four ratios used to determine $M_{GT}^{2\nu}$ are shown in Fig. 1 for both SS and SP source–sink combinations.

Fits are performed to the statistically more precise SP correlators and the values of the total matrix element and the short-distance contribution, normalised by the naive deuteron-pole matrix element $g_{A}^{2}/\Delta$”, are given by

[TABLE]

In these expressions, the first uncertainties arise from statistical sampling and from systematic effects from fitting choices and deviations from Wigner symmetry Tiburzi et al. (2017). The second uncertainties encompass differences between analysis methods. The leading discretisation effects, which are potentially large on the numerically smaller polarisability term, are removed by normalising to the square of the proton axial charge computed using the same lattice axial current on the same ensemble.

Discussion: The computed value of $M_{GT}^{2\nu}$ above can be used to determine the unknown EFT(${\pi\hskip-5.50003pt/}$) low-energy constant $\mathbb{H}_{2,S}$. Taking the values of $g_{A}$ and the two-body single-current matrix element from Ref. Savage et al. (2017), and using the calculated binding energies and effective ranges of the two-nucleon systems Beane et al. (2013a, b), the result is $\mathbb{H}_{2,S}=4.7(1.3)(1.8)\ {\tt fm}$. The dominant contribution to $M_{GT}^{2\nu}$ comes from the deuteron pole with coupling $g_{A}^{2}$. This is modified by two-body effects in the axial current ($\tilde{l}_{1,A}$); this contribution shifts the leading-order result from $g_{A}^{2}/\Delta$ to $|\langle pp|\tilde{J}_{3}^{+}|d\rangle|^{2}/\Delta$, which is approximately a 5% shift. Interestingly, the calculated result suggests that the additional axial polarisability contribution is of similar size. The existence of this short-distance contribution precludes accurately predicting $\beta\beta$-decay matrix elements in a nuclear many-body calculation by simply rescaling (quenching) $g_{A}$.

The present results are obtained at an unphysical quark mass without the inclusion of electromagnetism and isospin breaking effects, and at a single lattice spacing and volume. All these caveats may be important, particularly given that the short-distance two-nucleon effects are only few-percent contributions to the matrix elements. Such effects require further investigation but are not expected to qualitatively alter the conclusions of this work. Despite these qualifications, the potential for a relatively large contribution from the isotensor axial polarisability is important, as terms of this form are not included in current phenomenological analyses of $\beta\beta$-decay. This observation, supported and motivated by the numerical calculations, is the central point of this Letter. In order to accurately predict $2\nu\beta\beta$ decay rates with fully quantified uncertainties, the isotensor axial polarisabilities of nuclei must be determined. Future lattice QCD calculations in few-nucleon systems and light nuclei, matched to nuclear many-body methods such as EFT(${\pi\hskip-5.50003pt/}$), offer the possibility of determining these contributions which are difficult to access experimentally. However, in order to undertake such calculations at the physical quark masses, a number of difficulties related to the bi-local nature of such weak processes must be overcome Tiburzi et al. (2017). Additional complications arise for $0\nu\beta\beta$ where these polarisability contributions are also likely to be important in the case of light Majorana neutrinos. Furthermore, the short-distance strong interaction contributions encapsulated in the isotensor axial polarisability provide an inherent background to extracting contributions from short-distance lepton-number violating operators that are possible beyond the light Majorana scenario Savage (1999); Prezeau et al. (2003); Graesser (2016); Nicholson et al. (2016); Cirigliano et al. (2017).

Acknowledgments: This research was supported in part by the National Science Foundation under grant number NSF PHY11-25915 and ZD, WD, MJS, PES, BCT and MLW acknowledge the Kavli Institute for Theoretical Physics for hospitality during completion of this work. Calculations were performed using computational resources provided by NERSC (supported by U.S. Department of Energy grant number DE-AC02-05CH11231), and by the USQCD collaboration. This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract number DE-AC05-00OR22725. The PRACE Research Infrastructure resources at the Très Grand Centre de Calcul and Barcelona Supercomputing Center were also used. Parts of the calculations used the chroma software suite Edwards and Joó (2005) and the quda library Clark et al. (2010); Babich et al. (2011). EC was supported in part by the USQCD SciDAC project, the U.S. Department of Energy through grant number DE-SC00-10337, and by U.S. Department of Energy grant number DE-FG02-00ER41132. ZD, WD and PES were partly supported by U.S. Department of Energy Early Career Research Award DE-SC0010495 and grant number DE-SC0011090. The work of WD is supported in part by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, within the framework of the TMD Topical Collaboration. KO was partially supported by the U.S. Department of Energy through grant number DE-FG02-04ER41302 and through contract number DE-AC05-06OR23177 under which JSA operates the Thomas Jefferson National Accelerator Facility. MJS was supported by DOE grant number DE-FG02-00ER41132, and in part by the USQCD SciDAC project, the U.S. Department of Energy through grant number DE-SC00-10337. BCT was supported in part by a joint City College of New York-RIKEN/Brookhaven Research Center fellowship, and by the U.S. National Science Foundation, under grant number PHY15-15738. MLW was supported in part by DOE grant number DE-FG02-00ER41132. FW was partially supported through the USQCD Scientific Discovery through Advanced Computing (SciDAC) project funded by U.S. Department of Energy, Office of Science, Offices of Advanced Scientific Computing Research, Nuclear Physics and High Energy Physics and by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under contract DE-AC05-06OR23177.