Chiral Effective Field Theory Calculations of Weak Transitions in Light Nuclei

In this work, we present Quantum Monte Carlo (QMC) calculations, including both variational Monte Carlo (VMC) and Green’s function Monte Carlo (GFMC) calculations, of Gamow-Teller (GT) matrix elements entering beta-decay and electron-capture rates in $A\,$=$\,3$–$10$ nuclei. These observables are experimentally known (in most cases) at the sub-percent level, and are used here primarily to validate our microscopic theoretical modeling of the nucleus as a system of nucleons interacting with each other via effective interactions, and with electroweak probes via effective currents. Specifically, this modeling is based on local two- and three-nucleon interactions formulated in configuration space, and derived from a chiral effective field theory ($\chi$EFT) that retains, in addition to nucleons and pions, $\Delta$-isobars as explicit degrees of freedom Piarulli et al. (2015, 2016); Baroni et al. (2016a, 2018). They are referred to below as the Norfolk interactions and are denoted as NV2+3. Accompanying these interactions are one- and two-body axial currents—local in configuration space—derived within the same $\chi$EFT formulation Baroni et al. (2016b, a, 2018). In the calculations to follow, we include up to tree-level contributions at next-to-next-to-next-to-leading order (N3LO) in the chiral expansion, and disregard subleading corrections involving loops and higher order contact terms.

An analogous study was recently reported by some of the present authors in Ref. Pastore et al. (2018a). There, the QMC calculations were based on the Argonne-$v_{18}$ (AV18) two-nucleon Wiringa et al. (1995) and Illinois-7 (IL-7) three-nucleon Pieper (2008) interactions, in combination with the axial currents of Ref. Baroni et al. (2016b). We found agreement with the experimental Gamow-Teller matrix elements at the $\approx 2\%$ level for the $A\,$=$\,6$ and $7$ systems, and at the $\approx 10\%$ level in ¹⁰C. The large uncertainty in the $A\,$=$\,10$ transition is primarily systematic and results from the narrow energy separation between the first two $J^{\pi}\,$=$\,1^{+}$ states in ¹⁰B, which makes it hard to precisely disentangle them Pastore et al. (2018a). The study of Ref. Pastore et al. (2018a) found that two-body currents generate an additive contribution of the order of $\approx 3\%$ and concluded that the agreement with the data is mainly attributable to the use of fully correlated nuclear wave functions, rather than two-body effects in the currents.

In the meantime, no-core shell-model calculations of weak matrix elements based on chiral interactions and currents Gysbers et al. (2019) found the sign of the overall correction generated by two-body currents to be opposite to that obtained in Ref. Pastore et al. (2018a) for the same systems (but in agreement with a hybrid calculation of the $A\,$=$\,6$ decay reported in Ref. Gazit et al. (2009a)). This discrepancy was attributed to the hybrid nature of the calculation of Ref. Pastore et al. (2018a), i.e., to the mismatch between the two- and three-body correlations implemented to construct the nuclear wave functions—and induced by the AV18 and IL7—and those entering the axial currents which were instead derived from $\chi$EFT.

In this work, by reexamining the evaluation of these weak matrix elements with the NV2+3 chiral interactions Piarulli et al. (2015, 2016); Baroni et al. (2016a, 2018), in combination with consistent chiral axial currents at tree-level Baroni et al. (2016b), we aim to address and explore the aforementioned claim. We investigate the sensitivity of the calculated matrix elements with respect to different choices of regulators and to different strategies adopted to constrain the three-body Norfolk interactions (NV3). This latter aspect is important in order to understand the interplay between these interactions and the axial currents, since the strength of the contact current and that of the three-body interaction of one-pion-range are rigorously related to each other by the symmetries imposed in the $\chi$EFT formulation Gardestig and Phillips (2006); Gazit et al. (2009b); Schiavilla (2017).

This study has several merits. First, we report new GFMC results of the energy spectra of $A\leq 10$ nuclei based on two classes of NV2+3 interactions. In addition to the systems studied in Ref. Pastore et al. (2018a), we study weak transitions in $A\,$=$\,8$ nuclei where we find that two-body axial currents provide a large correction to the one-body results. Finally, we provide the first calculations of two-body weak transition densities which shed light on the role of short-range physics in these observables.

Searches for physics beyond the Standard Model (BSM) via beta decay are the focus of current and planned experimental programs carried out at the Facility for Rare Isotope Beams (FRIB), the the University of Washington, and Argonne National Laboratory (see, e.g., Ref. Gonzalez-Alonso et al. (2019) and references therein). Among the targets under consideration are ⁶He, ⁸Li, ⁸B and ¹⁰C. A systematic study of axial-current matrix elements in these systems is a prerequisite for all further investigations and BSM searches in beta decay. Beta decays are ideal processes to assess the validity of the dynamical inputs of ab inito calculations, namely many-body correlations and weak currents. The latter also impact calculations of neutrinoless double beta decay matrix elements, whose knowledge is critical to the neutrinoless double beta decay experimental program Engel and Menendez (2017).

This paper is structured as follows: a brief review of the QMC computational method and Norfolk interactions is given in Secs. II and III. The many-body axial currents used in this work are reported in Sec. IV. Results and conclusions are provided in Secs. V and VI.

The Quantum Monte Carlo methods used in this study have been described in detail in several review articles, the most recent of which being Refs. Carlson et al. (2015); Lynn et al. (2019). Here, we sketch the computational procedure and refer the interested reader to Refs. Carlson et al. (2015); Lynn et al. (2019) and references therein.

We seek accurate solutions of the many-nucleon Schrödinger equation

where $J^{\pi}$ are the total angular momentum and parity of the state, and $T$ and $T_{z}$ are the total isospin and its projection, respectively. We use the Hamiltonian

where $K_{i}$ is the non-relativistic kinetic energy, and $v_{ij}$ and $V_{ijk}$ are the NV2 and NV3 local chiral interactions Piarulli et al. (2015, 2016); Baroni et al. (2016a, 2018), collectively denoted as NV2+3.

The VMC trial function $\Psi_{V}(J^{\pi};T,T_{z})$ for a given nucleus is constructed from products of two- and three-body correlation operators acting on an antisymmetric single-particle state of the appropriate quantum numbers. The correlation operators are designed to reflect the influence of the interactions at short distances, while appropriate boundary conditions are imposed at long range Wiringa (1991); Pudliner et al. (1997); Nollett et al. (2001); Nollett (2001). The $\Psi_{V}(J^{\pi};T,T_{z})$ contains variational parameters that are adjusted to minimize the expectation value

which is evaluated by Metropolis Monte Carlo integration Metropolis et al. (1953). The lowest value for $E_{V}$ is then taken as the approximate ground-state energy of the exact lowest eigenvalue of $H$, $E_{0}$, for the specified quantum numbers.

A good trial wave function is given by

The Jastrow wave function $\Psi_{J}$ is fully antisymmetric and has the $(J^{\pi};T,T_{z})$ quantum numbers of the state of interest, while $U_{ij}$ and $\widetilde{U}^{TNI}_{ijk}$ are two- and three-body correlation operators.

The GFMC method Carlson et al. (2015) improves on the VMC wave functions by acting on $\Psi_{V}$ with the operator $\exp\left[-\left(H-E_{0}\right)\tau\right]$. The operator is applied in a sequence of small imaginary-time steps $\Delta\tau$ to produce a propagated wave function

Obviously $\Psi(\tau\!=\!0)\,$=$\,\Psi_{V}$ and $\Psi(\tau\rightarrow\infty)\,$=$\,\Psi_{0}$. Quantities of interest are evaluated in terms of a “mixed” expectation value between $\Psi_{V}$ and $\Psi(\tau)$:

where the operator $O$ acts on the trial function $\Psi_{V}$. The desired expectation values would, of course, have $\Psi(\tau)$ on both sides; by writing $\Psi(\tau)=\Psi_{V}+\delta\Psi(\tau)$ and neglecting terms of order $[\delta\Psi(\tau)]^{2}$, we obtain the approximate expression

where $\langle O\rangle_{\rm V}$ is the variational expectation value.

For off-diagonal matrix elements required by the transitions we are interested here, the generalized mixed estimate is given by the expression

where

and $\langle O(\tau)\rangle_{M_{i}}$ is defined similarly. For more details see Eqs. (19)–(24) and the accompanying discussions in Ref. Pervin et al. (2007).

We base our calculations of weak transitions in $A$=6–10 on the local NV2 and NV3 interactions developed in Refs. Piarulli et al. (2015, 2016); Baroni et al. (2016a, 2018). The NV2 model has been derived from a $\chi$EFT that uses pions, nucleons and $\Delta$’s as fundamental degrees of freedom. It consists of a long-range part, $v_{ij}^{\rm L}$, mediated by one- and two-pion exchanges, and a short-range part, $v_{ij}^{\rm S}$, described in terms of contact interactions with strengths specified by unknown low-energy constants (LECs). The strength of the long-range part is fully determined by the nucleon and nucleon-to-$\Delta$ axial coupling constants $g_{A}$ and $h_{A}$, the pion decay amplitude $F_{\pi}$, and the LECs $c_{1}$, $c_{2}$, $c_{3}$, $c_{4}$, and $b_{3}+b_{8}$, constrained by reproducing $\pi N$ scattering data Krebs et al. (2007). The LECs entering the contact interactions are fixed by fitting nucleon-nucleon scattering data from the most recent and up-to-date database collected by the Granada group Navarro Pérez et al. (2013, 2014); Navarro Perez et al. (2014). The value for $h_{A}$ is taken from the large $N_{c}$ expansion or strong-coupling model Green (1976). The value of the nucleon axial coupling constant used to construct the nucleon-nucleon interaction accounts for the Goldberger-Treiman discrepancy Arndt et al. (1994); Stoks et al. (1993) and, to distinguish it from the experimental value of $g_{A}\,$=$\,1.2723\,(23)$ Tanabashi et al. (2018) entering the axial currents, we denote it with $\tilde{g}_{A}$. For completeness, Tables 1 and 2 report the values of these constants, along with the pion and nucleon masses, the $\Delta$-nucleon mass difference, the electron mass, and the fine structure constant $\alpha$ used in the NV2 interactions (these last two characterize the electromagnetic part of the NV2s Piarulli et al. (2015)).

The contact terms are implemented using a Gaussian representation of the three dimensional delta function, with $R_{S}$ denoting the Gaussian parameter Piarulli et al. (2016, 2015); Baroni et al. (2016a, 2018). The pion-range operators are strongly singular at short range in configuration space, and are regularized by a radial function characterized by a cutoff $R_{L}$ Piarulli et al. (2016, 2015); Baroni et al. (2016a, 2018). There are two classes (I and II) of NV2s, differing only in the range of energy over which they are fitted to the database—class I up to 125 MeV, and class II up to 200 MeV. For each class, two combinations of short- and long-range regulators have been used, namely ($R_{S}$, $R_{L}$)=(0.8, 1.2) fm (models NV2-Ia and NV2-IIa) and ($R_{S}$, $R_{L}$)=(0.7, 1.0) fm (models NV2-Ib and NV2-IIb). Class I (II) fits about 2700 (3700) data points with a $\chi^{2}$/datum $\lesssim 1.1$ ($\lesssim 1.4$) Piarulli et al. (2016, 2015).

The NV2 models were found to provide insufficient attraction in GFMC calculations of the binding energies of light nuclei Piarulli et al. (2016). To remedy this shortcoming, a consistent three-body interaction was constructed up to N2LO in the chiral expansion. It consists of a long-range part mediated by two-pion exchange and a short-range part parametrized in terms of two contact interactions van Kolck (1994); Epelbaum et al. (2002) proportional to the LECs $c_{D}$ and $c_{E}$. These LECs have been obtained by fitting either observables that involve exclusively strong interactions Lynn et al. (2016); Tews et al. (2016); Lynn et al. (2017); Piarulli et al. (2018) or a combination of observables that involve both strong and weak interactions Gazit et al. (2009b); Marcucci et al. (2012); Baroni et al. (2018). This last strategy is feasible because of the relation established in $\chi$EFT Gardestig and Phillips (2006) that links $c_{D}$ with the LECs entering the contact axial current at N3LO Gazit et al. (2009b); Marcucci et al. (2012); Schiavilla (2017) (see next section for details).

In Ref. Piarulli et al. (2018), $c_{D}$ and $c_{E}$ were determined by simultaneously reproducing the experimental trinucleon ground-state energies and $nd$ doublet scattering length. These first-generation NV2+3 interactions, denoted with NV2+3-Ia/b and NV2+3-IIa/b, have been implemented in both VMC and GFMC codes and used to study static properties of light nuclei Piarulli et al. (2016, 2018); Lynn et al. (2019); Cirigliano et al. (2019); Piarulli and Tews (2020); Gandolfi et al. (2020), and in auxiliary-field diffusion Monte Carlo (AFDMC) Schmidt and Fantoni (1999), Brueckner-Bethe-Goldstone (BBG) Day (1967); Baldo and Burgio (2012) and Fermi hypernetted chain/single-operator chain (FHNC/SOC) Fantoni and Rosati (1975); Pandharipande and Wiringa (1979) approaches to investigate the equation of state of neutron matter Piarulli et al. (2019); Bombaci and Logoteta (2018).

In more recent work Baroni et al. (2018), $c_{D}$ and $c_{E}$ were constrained by fitting, in addition to the trinucleon energies, the empirical value of the GT matrix element in tritium $\beta$ decay. These second-generation NV2+3 interactions were designated as NV2+3-Ia^∗/b^∗ and NV2+3-IIa^∗/b^∗. These two different procedures for fixing $c_{D}$ and $c_{E}$ produced rather different values for these LECs. They are reported in Table 3.

Many-body axial currents have been first examined within $\chi$EFT by Park and collaborators in Ref. Park et al. (1993). In that work, the authors retained pions and nucleons in their effective theory and calculated the two-body axial currents up to one-loop terms. The derivation neglected, for example, pion-pole contributions. More recently, two-body axial currents with pions and nucleons have been derived by the Bonn group Krebs et al. (2017) using the unitary transformation method, and by the JLab-Pisa group using time-ordered perturbation theory Baroni et al. (2016b, a). The two derivations differ in the treatment of reducible diagrams. When calculating box diagrams entering the electromagnetic charge and current operators Pastore et al. (2008, 2009, 2011); Kolling et al. (2009, 2011), the two methods lead to results that are in agreement. However, as discussed at length in Refs. Baroni et al. (2016a); Krebs et al. (2020), the two groups find different results for the box diagrams in the two-body axial current operator at N4LO. The numerical impact of this difference has been investigated in Refs. Baroni et al. (2016a, 2018), where both the JLab-Pisa and Bonn versions of the N4LO current operators have been implemented to calculate the GT matrix element in triton beta decay. In those studies, it was found that the corrections generated by the JLab-Pisa and Bonn N4LO operators are qualitatively in agreement (they both quench the GT matrix element at leading order), and provide, respectively, a $\approx 6\%$ and $\approx 4\%$ contribution to the total GT matrix element.

Here, we consider two-body axial currents derived within the same $\chi$EFT used to construct the NV2+3 interactions Baroni et al. (2018). Moreover, we base our calculations on tree-level corrections only, and disregard the (problematic) N4LO loop contributions discussed above. This choice is advantageous also because it allows for a clearer comparison with the no-core shell model and coupled-cluster calculations of Ref. Gysbers et al. (2019), which are also based on tree-level axial currents alone (albeit derived in a $\Delta$-less $\chi$EFT). Corrections at N4LO in the present formulation are, in practice, subsumed in the LECs of the theory, which have been determined by fits to experimental data.

Before moving on to a (brief) discussion of these axial currents, it is worthwhile pointing out that many-body corrections to leading one-body transition operators have been shown to be crucial for providing a quantitatively successful description of many nuclear electroweak observables Bacca and Pastore (2014), such as nuclear electromagnetic form factors Piarulli et al. (2013); Schiavilla et al. (2019); Nevo Dinur et al. (2019); Marcucci et al. (2016), low-energy electroweak transitions Pastore et al. (2009); Girlanda et al. (2010); Pastore et al. (2011, 2013); Datar et al. (2013); Pastore et al. (2014); Pastore et al. (2018a), and electroweak scattering Pastore et al. (2019). They have also been used in studies of double beta decay matrix elements Pastore et al. (2018b); Cirigliano et al. (2018, 2019); Wang et al. (2019).

The N3LO axial currents used in this work are represented diagrammatically Fig. 1. We refer the interested reader to Refs. Baroni et al. (2016b, 2018) for additional details and explicit expressions of the operators; here, we only note that we do not show diagrams that lead to vanishing contributions as well as pion-pole terms which give negligible corrections to the matrix elements under study.

The LO term, which scales as $Q^{-3}$ in the power counting ($Q$ denotes generically a low-momentum scale), is shown in panel (a) and reads

where $g_{A}$ is the nucleon axial coupling constant ($g_{A}\,$=$\,1.2723$ Tanabashi et al. (2018)), $f_{\pi}$ is the pion-decay amplitude ($f_{\pi}\,$=$\,92.4$ MeV), ${\bm{\sigma}}_{i}$ and ${\bm{\tau}}_{i}$ are the spin and isospin Pauli matrices of nucleon $i$, ${\bf q}$ is the external field momentum, ${\bf r}_{i}$ is the position of nucleon $i$, and the subscript $a$ specifies the isospin component ($a\,$=$\,x,y,z$).

At N2LO there are two contributions (scaling as $Q^{-1}$). The first one is a relativistic correction to the single-nucleon operator at LO and is diagrammatically illustrated in panel (b), while the second involves the excitation of a nucleon into a $\Delta$ by pion exchange, as illustrated in panel (c). In the tables and figures below, we will denote these two contributions with N2LO-RC and N2LO-$\Delta$, respectively. We use the same notation introduced in Ref. Baroni et al. (2018) and write the cumulative N2LO contribution as

At N3LO (or $Q^{0}$ in the chiral expansion), there is a term of one-pion range illustrated in panel (d), and a contact term shown in panel (e), which together give the following N3LO correction

We will denote the individual terms with N3LO-OPE and N3LO-CT, respectively.

The configuration-space expressions of these currents are given in Eqs. (2.7)–(2.10) of Ref. Baroni et al. (2018). Here, we limit ourselves to report the expression of the N3LO contact term used in this work to explicitly show the relation between the LECs entering this axial current and the LEC $c_{D}$ in the three-nucleon interaction. In $r$-space the N3LO-CT current reads

where

and $r_{ij}$ is the interparticle distance. The $\delta$-function in the contact axial current has been smeared by replacing it with a Gaussian cutoff of range $R_{S}$,

as previously done for the contact-like terms of the NV2 interactions. The adimensional LEC $z_{0}$ (reported in Table 3) is given by

where $c_{D}$ is the LEC multiplying one of the contact terms in the three-nucleon interaction Epelbaum et al. (2002) given in Table 3, $\Lambda_{\chi}=1$ GeV is the chiral symmetry breaking scale, $c_{3}$ and $c_{4}$ are given in Table 1, and $m$ and $m_{\pi}$ are the average nucleon and pion masses. It has recently been realized Schiavilla (2017) that the relation between $z_{0}$ and $c_{D}$ had been given erroneously in the original reference Gazit et al. (2009b), a – sign and a factor 1/4 were missing in the term proportional to $c_{D}$.

The GFMC energies of the nuclei of interest calculated using the NV2+3-Ia and NV2+3-Ia* models are listed in Table 4 along with the dominant spatial symmetry (s.s.) of the variational wave functions Wiringa (2006). The energies are obtained using $\approx 80,000$ walkers, and are all well converged by 30 unconstrained steps Wiringa et al. (2000a). All the GFMC results presented in this article (but for the two cases discussed below) are averages over the imaginary time $\tau$ from 0.2 to 0.82 MeV^-1. Results obtained with the NV2+3-Ia interaction are in statistical agreement with those published in Ref. Piarulli et al. (2018) based on the same nuclear Hamiltonian. Model NV2+3-Ia leads to predictions that are in excellent agreement with the data. We also report for the first time results based on the second generation of NV2+3 interactions, specifically model NV2+3-Ia*, whose three-nucleon interaction has been constrained by fits to the experimental trinucleon binding energies and tritium GT matrix element. Results obtained with the NV2+3-Ia^∗ Hamiltonian display a somewhat less satisfactory agreement with the experimental data, but still less than 4% away from them.

A typical imaginary-time evolution of the GFMC transition matrix elements is shown in Fig. 2. As can be seen, there is a rapid drop of 3% from the initial VMC estimate at $\tau\,$=$\,0$ that reaches a stable value around 0.2 MeV^-1. The results for all transitions presented in this article are averages over $\tau$ from 0.2 to 0.82 MeV^-1, as indicated by the dashed lines, with statistical errors denoted by the solid lines. The calculations of weak transitions involving the $(J^{\pi},T)=(2^{+},0)$ state of ⁸Be and the ground state of ⁸B are treated differently. For these two states, we observe that the binding energy, magnitude of the quadrupole moment, and point-proton radius all increase monotonically as the imaginary time increases. This can be appreciated in Fig. 3 where we show the point-proton radii of the two nuclear states. We interpret this behavior as an indication that the resonant excited state of ⁸Be is dissolving into two separated $\alpha$’s, while ⁸B is breaking into $p+^{7}$Be. In the case of ⁸Be, this issue has been addressed already in Refs. Wiringa et al. (2000b); Pastore et al. (2014); Datar et al. (2013). Here, we use similar techniques to treat these systems and extract matrix elements from the GFMC data. In particular, we note that the ground-state energy of ⁸Be drops very quickly as the imaginary time increases and reaches stability around $\tau\sim 0.1$ MeV^-1. The transition matrix elements involving the two dissolving states have been determined assuming that, also for these states, $\tau\approx 0.1$ MeV^-1 is the point at which spurious contaminations in the nuclear wave functions have been removed by the GFMC propagation. We then average in a small interval around this point, typically between 0.06 and 0.14 MeV^-1. Calculations involving these two systems are clearly affected by a systematic error. To have a rather rough estimate of this error, we study the sensitivity of the extracted matrix elements with respect to variations in the imaginary time interval selected for the averaging. We find that such a procedure generates an additional uncertainty of $\approx 5\%$ which is added in quadrature to the statistical one, and quoted in Table 6 below of GFMC results.

In this work, we reported on a detailed study of weak matrix elements in $A\,$=$\,3$–10 systems based on chiral (two- and three-nucleon) interactions and associated (one- and two-body) axial currents at high orders in the chiral expansion. A summary of our results is displayed in Fig. 6. Agreement with the experimental data is obtained when correlated nuclear wave functions are adopted. For these transitions the contribution of corrections beyond LO in the axial current is typically at the $\approx 2\%$ level of the value calculated with the LO Gamow-Teller operator. These findings are in line with those reported in the hybrid study of Ref. Pastore et al. (2018a) for the same transitions. Here, we also present calculations of matrix elements entering the rates of the ⁸Li, ⁸B, and ⁸He beta decays. These matrix elements are found to be suppressed at LO, and N2LO and N3LO currents provide a large correction ($\approx 20$–30%) which is, however, insufficient to explain the experimental data. We attribute the large suppression at LO to the fact that the Gamow-Teller operator is, in these transitions, connecting large to small components of the initial and final wave functions. Improving on these calculations will require the development of more sophisticated wave functions with better constrained small components.

Finally, we also reported on a careful analysis of one- and two-body transition densities shown in Figs. 4–8. The latter are especially interesting because they allow us to understand the spatial distributions of the various two-body current operators, that is their behavior as functions of interparticle distance. We have shown that, for each set of interactions and consistent currents (either NV2+3-Ia or NV2+3-Ia*), the two-body transition densities exhibit a universal behavior at short distance across all nuclei we have considered in the present study.

We would like to thank V. Cirigliano, R.  Charity, W. Dekens, J. Engel, A. Hayes, E. Mereghetti, and L.  Sobotka for useful discussions at various stages of this work. S. P. and G. K. would like to thank Grigor Sargsyan for useful discussions and for sharing information on the experimental log$(ft)$ value of ⁸Li. The many-body calculations were performed on the parallel computers of the Laboratory Computing Resource Center, Argonne National Laboratory, the computers of the Argonne Leadership Computing Facility (ALCF) via the 2019/2020 ALCC grant “Low energy neutrino-nucleus interactions” for the project NNInteractions. The work of J.C., S.G., and R.B.W. has been supported by the NUclear Computational Low-Energy Initiative (NUCLEI) SciDAC project. This research is supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under contracts DE-AC05-06OR23177 (R.S.) and DE-AC02-06CH11357 (R.B.W.), DE-AC52-06NA25396 (S.G. and J.C.) and U.S. Department of Energy funds through the FRIB Theory Alliance award DE-SC0013617 (M.P. and S.P.).