Beta Decay in Medium-Mass Nuclei with the In-Medium Similarity Renormalization Group

We review the status of ab initio calculations of allowed beta decays (both Fermi and Gamow-Teller), within the framework of the valence-space in-medium similarity renormalization group approach.


Introduction
Beta decays in atomic nuclei have long been a source of fundamental discoveries in physics [1][2][3], and precise measurements of beta decays continue to be a promising path to search for physics beyond the Standard Model (BSM) [4][5][6][7]. A major challenge in the search for a signal of new physics is understanding the Standard Model "background", especially the effects of low-energy quantum chromodynamics which manifest as nuclear structure. The situation is aggravated by the fact that nuclei which are preferred experimentally are often difficult to treat theoretically in a framework that allows quantified uncertainties.
Nevertheless, progress has been made over the past few decades so that the internucleon interaction can be systematically constructed within an effective field theory framework [8][9][10]. Simultaneously, advances in many-body theory and computational resources have enabled ab initio treatment of the medium-mass nuclei which are often relevant for BSM searches [11][12][13][14][15][16][17]. Of course more work remains to be done, both on the effective field theory side [18][19][20], and on understanding how approximation schemes in ab initio calculations impact the observables in question.
In this paper, I will focus on one particular many-body method-the valence-space inmedium similarity renormalization group (VS-IMSRG)-and consider two topics in allowed beta decay, the quenching in Gamow-Teller decays and correction factors for superallowed 0 + → 0 + Fermi decays.

IMSRG formalism
There are several review articles detailing both the free-space SRG [21,22] and the inmedium SRG [15,[23][24][25][26], and so here I will review only what is needed for our present purposes.

Similarity renormalization group
The basic idea of the SRG is to perform a unitary transformation U on the Hamiltonian H (and all other operators) in such a way that the resulting nuclear wave function is simpler. This is achieved by performing a sequence of infinitessimal unitary transformations, labeled by a flow parameter s, so that H(s) = U(s)H(0)U † (s) (1) with U(0) = 1. The way in which U changes with s is determined by the generator η, which we are free to specify, dU(s) ds = η(s)U(s) (2) arXiv:2109.13462v1 [nucl-th] 28 Sep 2021 as long as η is anti-hermitian η † = −η. Combining (1) and (2) The flow equation for any other operator O is obtained by replacing H → O in (3) [27][28][29][30]. It remains to specify η(s). In the free-space SRG, we choose where T is the kinetic energy and V(s) is the potential so H(s) = T + V(s). This generator drives V(s) towards a band-diagonal form in momentum space, with a width λ SRG ≡ s −1/4 . When the SRG flow equation (3) is formulated in Fock-space (i.e. in terms of creation and annihilation operators), many-body forces are inevitably induced, and these must be truncated in order to make the calculation tractable. For this reason, the "free-space" SRG evolution is typically performed out to λ SRG 2fm −1 .

In-medium SRG
The truncation of many-body forces is rendered less severe if all operators are normalordered with respect to a reference state |Φ , which should be a reasonable first approximation of the exact wave function |Ψ . This approach is called the in-medium SRG (IMSRG). If we choose the generator to suppress the parts of H which lead to excitations out of the reference |Φ , then for s → ∞ |Φ becomes an eigenstate of H(s) with an eigenvalue corresponding to the energy of the exact wave function |Ψ , up to approximation errors in solving (3). In all calculations presented here, I neglect three-body operators after the initial normal-ordering step, resulting in the IMSRG(2) approximation.
One possible choice for the generator which achieves the desired suppression was proposed by White [31] where the "off-diagonal" part of the Hamiltonian, denoted H od , is any part of H which connects |Φ to a different state. The energy denominator ∆ is typically defined with Epstein-Nesbet or Møller-Plesset partitioning. In this work, I use Epstein-Nesbet denominators, and a modification of (5)-also suggested by White-called the arctangent generator (see Ref. [26] for more details). As a further generalization, we may define a valence space (e.g. the sd-shell above an 16 O core), and define H od such that any part of H which connects a valence configuration to a non-valence configuration is suppressed 1 . The Hamiltonian is then driven to a blockdiagonal form and we may diagonalize in the (typically much smaller) sub-space of valence configurations. Such a diagonalization directly corresponds to a standard large-scale shell model calculation with an effective interaction defined by H(∞). This approach is referred to as the valence-space IMSRG (VS-IMSRG), and it is used for all calculations presented.
Generally, the states we wish to target in a valence space approach are not well-described by a single closed-shell configuration, and the choice of |Φ becomes less clear. In this work I use the ensemble normal-ordering (ENO) approach [32], which amounts to taking fractional occupation numbers such that the reference has spherical symmetry and the correct number of particles on average. In addition, I use the Magnus formulation of the IMSRG [33], in which we write U(s) ≡ e Ω(s) (6) where Ω = −Ω † is the Magnus operator. Equations (6) and (2) In (7), each commutator is truncated at the normal-ordered two-body level, and the series is computed iteratively until the size of a term falls below a numerical threshold.

Aspects relevant to beta decay
Two additional details of the calculation are relevant for β decays. The first is the choice of single-particle basis in which we express the operators at s = 0. In this work I use a Hartree-Fock basis with Coulomb and isospin-breaking strong forces included, so that for a given set of single-particle quantum numbers {n, , j}, the proton and neutron radial wave functions are not identical. The second is the choice of reference. Because we wish to compute the initial state, with N neutrons and Z protons, consistently with the final state with N ± 1 neutrons and Z ∓ 1 protons, there is some ambiguity about which reference |Φ should be used. If we retain all induced many-body terms during the SRG evolution, the choice of reference is irrelevant. However, the accuracy of the IMSRG(2) approximation depends on the choice of reference. The two natural choices for β decay are to use the N, Z of the initial state or the final state. I will discuss this in more detail in Section 4.

Gamow-Teller decays
In a Gamow-Teller decay, the leptons carry one unit of angular momentum and leave the parity of the nucleus unchanged. The relevant nuclear transition operator is obtained from the space-like part of the hadronic axial-vector current. The leading term in the non-relativistic reduction is g A στ, where g A ≈ 1.27 is the axial coupling constant, and σ and τ are the spin and isospin Pauli matrices.
Historically, when the leading operator was combined with shell model wave functions, a systematic "quenching" of the decay strength was observed, i.e. experimental matrix elements were smaller than the predicted ones, with a similar effect in isovector M1 observables [34][35][36][37][38][39]. It was quickly surmised that the source of the discrepancy should be the some combination of inadequate wave functions (missing correlations) and an inadequate transition operator (missing currents), that neither of these obviously dominated, and that the two effects were not independent [35,40,41]. It was also suggested that pions ought to have something to do with the renormalization of the axial current in the nuclear medium [42].
These physics arguments survive in the modern EFT point of view, which organizes the nuclear interaction and coupling to external fields in powers of a ratio of scales. The distinction between short and long-distance physics is made by a cutoff, and the arbitrariness of the cutoff is reflected in the requirement that observables be independent of its value. The relationship between pions and the axial current arises as a consequence of broken chiral symmetry [43]. Importantly, chiral EFT enables a systematic and consistent construction of three-nucleon forces and two-body currents [44][45][46][47]. The result in the limit that the momentum carried by the leptons vanishes, up to order Q 0 (leading order is Q −3 ) is [44,48,49] where p i , p i are the incoming and outgoing momenta of the ith nucleon, and f π is the pion decay constant. The low-energy constants c 3 , c 4 and c D also enter into the NN and 3N forces, and so are not additional free parameters. Equations (9), (10), and (11) correspond to diagrams (a), (b), and (c) in Figure 1, respectively. Note that there are also corrections to the one-body operator of order p 2 i /m 2 N . Depending on how the nucleon mass is counted, these corrections will enter at different orders. In the counting of e.g. Park et al [44], these corrections are also Q 0 , while in the counting used by other authors [47,49,50], including the calculations in this paper, these corrections are Q 1 .
In Refs. [48] and [49], these currents were normal ordered with respect to uniform nuclear matter to obtain an in-medium quenching factor for the one-body operator. In Ref. [51], the full two-body current was constructed, consistently 2 with the NN+3N force, and Gamow-Teller decays of 14 C, 22 O and 24 O were computed using the coupled cluster method 3 . In all three of these cases, a quenching of about the right size was obtained. In Ref [55], axial currents up to N 4 LO were used in quantum Monte Carlo calculations of A=6-10 nuclei, where it was found that correlations beyond the shell model accounted for most of the quenching, with subleading currents playing a minor role. In Ref. [50], the full two-body current up to N 3 LO was constructed consistently with the NN+3N force, consistently SRG evolved, and evaluated in a range of nuclei in the p, sd, and p f shells, as well as 100 Sn, using no-core shell model, coupled cluster, or VS-IMSRG to solve the many-body problem. Here, I will provide some additional calculations not presented in [50], and some further discussion.
The experimental Gamow-Teller matrix elements are obtained from the f t values by with K ≡ (2π 3h7 ln 2)/(m 5 e c 4 ), and K/G 2 V ≈ 6140s. The Gamow-Teller matrix element is defined as Note different definitions have been used in the literature, e.g. one may divide the right hand side by g A , as was done in [50]. I consider Gamow-Teller transitions in nuclei in the p, sd and p f shells, with experimental data taken from refs. [37][38][39]. I have selected transitions with large transition matrix elements, with the goal of reducing sensitivity to fine-tuned cancellations. I also consider the decay of 100 Sn, which was treated with equations-of-motion coupled cluster in [50], and for which the experimental picture is still somewhat conflicted [56][57][58][59]. I adopt the average value presented in [59]. In the VS-IMSRG calculation of 100 Sn, I use valence space consisting of the 0 f 5 , 1p 3 , 1p 3 , 0g 9 orbits for protons and 0 f 7 , 1d 5 , 1d 3 , 2s 1 , 0h 11 for neutrons.
For the theoretical calculations, I use the NN+3N(lnl) interaction developed by Navrátil [60]. The interaction and current are consistently SRG-evolved to a scale λ SRG = 2.0 fm −1 and evaluated in an oscillator space defined by 2n + ≤ e max = 12 andhω = 16. The 3N matrix elements are further truncated with e 1 + e 2 + e 3 ≤ E 3max = 14. All operators are transformed to the Hartree-Fock basis, and then the residual 3N operators are truncated (the NO2B approximation). Next, a VS-IMSRG calculation is performed using the code imsrg++ [61], yielding and effective valence space interaction and operator. The valence space diagonalization is carried out either using NuShellX@MSU [62] with operators evaluating using the code nutbar [63], or with KSHELL [64,65]. The results are listed in Table A1 and plotted in Fig. 2. Table A1 in the appendix contains the numerical results. The column labeled M(GT) exp lists the experimental Gamow-Teller matrix elements defined by (13) (experimental uncertainties are not listed). The column labeled στ bare is the obtained by evaluating the operator στ (assuming identical radial wave functions for protons and neutrons) between valence space wave functions obtained using the VS-IMSRG evolved interaction. The column labeled στ IMSRG is obtained by consistently SRG and VS-IMSRG evolving the στ operator (including the radial mismatch due to the Hartree-Fock basis). Finally, M(GT) th also includes the twobody currents, consistently SRG and VS-IMSRG evolved. In a few cases, the listed strength is summed over multiple final states with the same spin and parity.
In Fig. 2, panel (a) shows a scatter plot of M(GT) exp vs M(GT) bare , while panel (b) shows MGT) th vs M(GT) exp . The solid line shows y = x corresponding to perfect agreement between theory and experiment. The dashed line shows a best-fit slope, which is indicated as a quenching factor at the top of the figure. For this quenching factor, I only include sd and p f shell nuclei because the p shell nuclei have a large scatter due to nuclear structure details. The quantity in parenthesis indicates the standard deviation about the best-fit line. If I include p shell nuclei in the fit, the full theory quenching factor changes to q = 0.99, but the standard deviation increases to 0.21.
It is evident from Table A1 that both the correlations included in στ IMSRG and the twobody currents lead to a reduction of the Gamow-Teller matrix element. As discussed in Ref. [50], the detailed breakdown of the quenching into correlations and currents is schemeand scale-dependent; some of the effects attributed to correlations when using a hard interaction get shuffled into currents, when using a soft interaction.
It is also evident from Figure 2 that the systematic quenching effect, observed when using the bare στ operator, essentially vanishes when using the consistently evolved operator including two-body currents.
In the right panel of Figure 2 I highlight the transition 12 N→ 12 C 0 + as an illustration of the cancellation effects in the p shell which wash out the quenching signal. When evaluating the bare στ operator 4 (including g A ), there are four terms that contribute, corresponding to proton-to-neutron transitions p 3 → p 3 , p 3 → p 1 , p 1 → p 3 , and p 1 → p 1 . The contributions are +0.366, -0.955, +1.592, and -0.173, respectively, totalling to 0.830. Evidently there is significant cancellation so that a relatively small change of the individual terms can lead to a relatively large change on the final matrix element. For example, if I use valence-space wave functions obtained with the phenomenological Cohen-Kurath interaction [67], the bare operator yields a matrix element of 1.219. This difference, due to configuration mixing within the valence space, is larger than the systematic quenching effect of interest. We may view this "noise" in the quenching from another perspective. Fig. 3 shows the theory-to-theory quenching factor M(GT) th /M(GT) bare as a function of mass number. This indicates the quenching factor needed if we wanted to approximately account for the correlations and currents contained in M(GT) th . The main point here is to emphasize that the quenching is not a smooth function of A, but in fact has considerable state dependence.
Moving forward, including two-body currents consistently with the NN+3N interaction helps to remove the ambiguity of emprical quenching factors, especially in the context of double beta decay [68]. Of course, sensitivity to details of the configuration mixing must be checked.

Superallowed 0 + → 0 + Fermi decays
For transitions between J = 0 states, B(GT) = 0 by conservation of angular momentum. Furthermore, in the limit in which isospin is a perfect symmetry, a "superallowed" transition between T = 1 isobaric analogue states yields B(F) = 2, and so (12) reduces to 5 (where here f = f V ). This implies all superallowed 0 + → 0 + should have the same f t value, and that from this one may measure the coupling constant for semileptonic decay G V , which is in turn related to the constant G F obtained from muon decay by G V = V ud G F , where V ud is the up-down element of the Cabibbo-Kobayashi-Maskawa (CKM) quark mixing matrix.
Consequently, precise f t measurements of superallowed 0 + → 0 + decays provide a sensitive test of the Standard Model: non-universality of superallowed f t values, or non-unitarity of the CKM matrix would be signs of new physics. Of course, isospin is not a perfect symmetry of the Standard Model. It is broken by the quark electric charges, and the up-down mass difference. This is manifested at the nuclear level as the Coulomb force between protons and isospin-violating strong interactions. The Standard Model corrections to (14) have been parameterized by Towner and Hardy [6] as In (15) ∆ V R is a process-independent radiative correction [69], δ R is a radiative correction only depending on the electron energy and the charge of the daughter nucleus, and δ NS is a radiative correction depending on the detailed nuclear structure. The isospin-symmetrybreaking correction δ C accounts for the fact that the final state is not exactly an isospin rotation of the initial state.
Consequently, only δ NS and δ C are the purview of nuclear structure theory. To draw an analogy with the situation for Gamow-Teller decays, δ C corresponds to including the effects of correlations for the leading operator τ, while the radiative corrections correspond to sub-leading corrections to the operator, with δ NS corresponding to two-body currents. The difference here is that the corrections are sub-leading in the fine structure constant α ≈ 1/137 (or Zα), as opposed to the chiral EFT expansion parameter Q ∼ 1/4. The various corrections are illustrated in Fig. 4.
In this paper we focus on the δ C correction, for no better reason than the operator is the simplest to implement. Towner and Hardy decompose δ C into a correction due to isospinbreaking configuration mixing effects, and a correction due the mismatch in single-particle  wave functions between protons and neutrons. As we will be treating both within a consistent calculation, such a decomposition is not necessary (and ambiguous) and we will simply use where M F = Ψ f τ Ψ i is the result of the ab initio calculation. Nevertheless, it is useful to keep the two mechanisms (configuration mixing and wave function mismatch) in mind when considering the impact of various approximations. The wave function mismatch effect is taken into account primarily by the fact that we use a Hartree-Fock single-particle basis, with Coulomb and nuclear ISB effects included in the potential (see also [70]).
To get an idea of what the VS-IMSRG framework produces for the isospin-breaking correction δ C , I consider three transitions spanning the p, sd and f p shells: 14 O→ 14 N, 34 Ar→ 34 Cl, and 46 Cr→ 46 V. I take the 1.8/2.0 (EM) interaction [71] with oscillator frequencyhω = 16 MeV and E 3max = 16. The resulting δ C values are plotted in Fig. 5 as a function of the e max truncation. I show results where the normal-ordering reference |Φ is taking to be either the inital or the final nucleus, and also the results from including only the one-body part of the evolved operator. For reference, I also indicate the δ C values adopted by Towner and Hardy [6].
If the calculation were under control, we should observe the following: convergence with respect to e max ; independence of the choice of reference; and a relatively small correction from including induced two-body terms, indicating a converging hierarchy of the cluster expansion. For the 46 Cr→ 46 V transition, we observe reference independence and a small two-body correction, but only a hint of convergence in e max . For the lighter nuclei the situation is worse, especially for 14 O→ 14 N. In all cases, it appears that the large e max behavior will need to be incorporated in some manner, possibly by utilizing natural orbitals [72,73], or by obtaining an extrapolation formula [74,75]. However, I leave this for future work. It appears that, at least for the near term, such calculations will tell us more about the IMSRG than about physics beyond the Standard Model.
The transition 14 O→ 14 N warrants a closer inspection, because it is light enough that it can be benchmarked against the no core shell model [77], although the observed e max dependence suggests converged results may be challenging. Moreover, the dramatic reference-dependence and contribution of two-body terms make this a good system for studying such effects, which are also important 6 in neutrinoless double beta decay [68,[78][79][80]. The two references used are shown schematically in Fig. 6.
Because 14 O is a closed shell, it may also be treated by the coupled cluster method [12] with the transition handled by the isospin-breaking equation-of-motion approach [51]. The strong e max dependence suggests sensitivity to infrared physics, and this may be probed by including the continuum in the coupled cluster calculation [81]. The results at e max = 12 are  [6]. For the 14 O decay, we also include coupled cluster points [76].  indicated with yellow diamonds in Fig. 5. The large δ C correction obtained with coupled cluster, as well as the substantial effect of including the continuum reinforces the notion that the pathology is not unique to the IMSRG. The reference dependence persists even down to e max = 3, which is sufficiently small that we may directly perform truncated configuration interaction (CI) calculations and extrapolate to the full CI result. At e max = 3, the VS-IMSRG yields δ C values of 0.839% and 0.056% with 14 O and 14 N references, respectively. The full CI result is 0.081%, which is considerably closer to the 14 N reference value. This reinforces our suspicion that the calculation with the 14 O reference is misbehaving.
As I mentioned above, results would be independent of the reference if all induced many-body terms were retained during the IMSRG evolution. In practice we only retain up to two-body operators, so any reference dependence indicates the impact of discarded three-body or (if we are especially unlucky) higher-body terms. If I truncate the expansion (7) for the transformed operator at two nested commutators, the reference dependence remains, while truncating at one nested commutator eliminates the effect. This suggests that three-body operators-which first show up at two nested commutators if O(0) is purely one-body-are the culprit.
A more complete investigation would benefit from calculations at the IMSRG(3) level, which are becoming available [82], but which will not be pursued here. For the moment, I will speculate. It appears that using the 14 O reference produces the largest error on the δ C value out of all the calculations presented in Fig. 5. This is naively counterintuitive because 14 O is a closed shell and should be best approximated by the single-configuration reference. But we are interested in the extent to which isospin symmetry is violated. Choosing a reference that approximates 14 O well but does a poor job for 14 N artificially breaks this symmetry, leading to an overestimate of the correction δ C . On the other hand, the ensemble 14 N reference used is not a great approximation of the state of interest in 14 N, which is actually open shell. Using the 14 N ensemble reference, we make an error in our description of 14 N, and we make a similar error in the description of 14 O, so the artificial breaking of isospin symmetry is reduced. In support of this, using a 12 C reference, which is a poor approximation of both 14 O and 14 N, results in δ C values in very close agreement to those obtained with the 14 N reference. This is in line with the behavior observed with other relative quantities, namely excitation energies and separation energies. It is a robust finding that the VS-IMSRG predicts too-high 2 + excitation energies for closed-shell nuclei [83,84]. This can be understood by considering that the reference is a good approximation of the 0 + ground state, and a worse approximation for the 2 + excited state. The truncation of three-body terms has a more severe impact on the 2 + state, leading to missed correlation energy, and consequently an excitation energy which is too high. For open-shell nuclei, the reference is a mediocre approximation of both the ground and excited states, and so they are treated on more equal footing, leading to a more accurate excitation energy. Likewise, it was observed that separation energies are more accurately obtained when using the same valence space for both the mass A and mass A − 1 nuclei, even if a different valence space might produce a more accurate ground state for one of the nuclei [85].
Preliminary studies have found that including three-body terms ameliorates the issue with 2 + states, and an analogous effect is observed in coupled cluster [86]. Presumably retaining three-body terms would also reduce the dependence of separation energies on the choice of valence space, but this has not yet been explored. One hopes that then including three-body terms will also help with the δ C calculation.

Conclusion and outlook
We are in an exciting time in nuclear structure theory, in which it is becoming possible to address issues which were long plagued by ambiguities arising from inconsistent modelling. On the question of Gamow-Teller decays, considerable progress has been made and it appears the issue is understood. However, in order to fully put the question to rest one should properly assess theoretical errors from the EFT truncation and the many-body solution, and demonstrate systematic improvement in both.
In addition, a similar quenching is observed in strong-interaction charge exchange reactions [87], where it is assumed that the transition operator for the target nucleus is proportional to στ. The quenching is often ascribed to missing strength at higher energies (the equivalent of "correlations" in beta decay). However, three-body forces should play an analogous role as two-body currents do for beta decay (replace the axial current in the diagrams in Figure 1 with an additional nucleon). It will be interesting to quantitatively investigate this parallel.
For superallowed Fermi decays, the work of Towner and Hardy [6] has laid a clear path, but there is more work to do on the many-body side. A more detailed understanding of how errors creep into the calculations will be essential.
Funding: This work was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, contract no. DEAC02-06CH11357 Table A1. Matrix elements for the Gamow-Teller transitions plotted in Figures 2 and 3. In the last column, † indicates the lowest 3 states of the listed J π in the final nucleus are summed in the matrix element, while * indicates that 4 states are summed.