Beta-decay formulas revisited (I): Gamow--Teller and spin-dipole contributions to allowed and first-forbidden transitions

We propose formulas of the nuclear beta-decay rate that are useful in a practical calculation. The decay rate is determined by the product of the lepton and hadron current densities. A widely used formula relies upon the fact that the low-energy lepton wave functions in a nucleus can be well approximated by a constant and linear to the radius for the $s$-wave and $p$-wave wave functions, respectively. We find, however, the deviation from such a simple approximation is evident for heavy nuclei with large $Z$ by numerically solving the Dirac equation. In our proposed formulas, the neutrino wave function is treated exactly as a plane wave, while the electron wave function is obtained by iteratively solving the integral equation, thus we can control the uncertainty of the approximate wave function. The leading-order approximation gives a formula equivalent to the conventional one and overestimates the decay rate. We demonstrate that the next-to-leading-order formula reproduces well the exact result for a schematic transition density as well as a microscopic one obtained by a nuclear energy-density functional method.


Introduction
The physics of exotic nuclei away from the stability line has been a major subject in nuclear physics. The lifetime of neutron-rich nuclei is governed by beta decay. Since the beta decay determines the time scale of the rapid-neutron-capture process (r-process) and the production of heavy elements together with the beta-delayed neutron(s) emission, the beta-decay rates of exotic nuclei are an important microscopic input for the simulation of nucleosynthesis [1]. The multi-messenger observations from a binary neutron star merger [2,3] imply that heavy neutron-rich nuclei that are even close to the drip line are involved in the r-process. Thus, the Coulomb effect on the beta particle (emitted electron) should be carefully examined under the extreme environment where the Q value for the beta decay, Q β , is high and the nuclear charge Z is large.
A careful analysis of the Coulomb effect is also useful for a precision test of the standard model to find a signal of new physics. For example, the effect on spectra of the beta particle and angular correlation as well as beta-decay rates has been studied to test the unitarity of the Cabibbo-Kobayashi-Maskawa (CKM) matrix, the scalar and tensor interactions, and the effect of neutrino mass in the allowed and first-forbidden transitions [4][5][6][7].
The formulation of nuclear beta decay within the distorted-wave impulse approximation of the electron Coulomb interaction has been matured [8][9][10][11][12][13][14][15]. The crucial part is how to handle the electron Coulomb wave function with a potential of the finite-size nuclear-charge distribution. Using the Maclaurin expansion of the nuclear radius r, the exact electron wave function was included in Ref. [16]. An iterative solution of the integral equation was found to have a better convergence by Behrens and Bühling [10,17]. The formula is arranged in the order of O(r a V b C E c e m d e ), where V C , E e and m e represent the Coulomb potential, the energy and the mass of an electron, respectively. It has been widely used in the calculations such as in Refs. [18][19][20] and in the recent application to the r-process nuclei [21][22][23][24][25][26][27][28][29][30][31][32]. In most of those works, however, the leading-order approximation of the formula in Refs. [10,12] is adopted. Instead of expanding the lepton wave functions, one can incorporate the numerical solution of the charged lepton wave functions thanks to the advance of the computational ability. The muon capture [33] and the beta decay [15] are formulated suitable for this purpose. In this formulation, the nuclear matrix element is defined in a transparent way and appears and S KLJ (κ ′ , κ) = 2(2j κ + 1)(2j κ ′ + 1)(2l κ + 1)(2l κ ′ + 1)(2K + 1)(l κ , 0, l κ ′ , 0|L, 0) Using the above form of the effective Hamiltonian, the beta-decay rate is given by integrating the scattering angles of the neutrino and electron as J,L,κe,κν where J i is the angular momentum of the initial state. Neglecting the mass of a neutrino, the maximum energy of an electron E 0 is the Q value of the nuclear transition, Q β . See the Appendix A for the derivation.
3 Electron Coulomb Wave function

Parametrization of Lepton Wave Function
The general formula given in Eq. (12) is ready for the use of any allowed and forbidden transition rates by evaluating the nuclear transition density. However, an explicit formula for the allowed and first-forbidden transitions helps extract nuclear structure information from the beta-decay observables. Since the electron Coulomb wave function is rather involved in evaluating the beta-decay rate, we briefly describe the derivation of the expression of charged lepton wave functions by iterating the integral equation following Refs. [10,17].
A Dirac wave function of an electron is given as The electron wave functions G κ , F κ satisfy the coupled first-order differential equation with the Coulomb potential V C : Throughout this paper we keep the electron mass explicit so that in future we can use the formula for the muon neutrino reactions. Electron wave functions are parametrized by taking 6 into account the behavior of the wave function at the origin r ∼ 0 [10] as Here k > 0 and H k (0) = 1 and h k (0) = 0. The normalization of the electron wave functions are determined by constants α κ . This parametrization of G κ , F κ incorporates the boundary condition of the wave function at the origin. Then the following set of coupled integral equations is obtained: At this stage R is just a parameter of dimension length. We take R as the nuclear radius though the final formulas are independent of the choice of R.

Iterative Solution of Integral Equation
Taking into account the boundary condition, H k , h k , D k , and d k are expanded according to the number of iteration as 7 The first iteration of the integral equation gives and further iterations give and for n = 1, 2, . . . . The exact electron wave functions in terms of E e , m e , and V C are obtained from the iterative solution of the above equations.

LO and NLO electron wave functions
We denote the leading order (LO) electron wave function as Adding the next-to-leading order (NLO), the NLO wave function is given as The explicit expressions of H The beta-decay rate is usually expressed in terms of the Fermi function F (Z, E e ) and the shape correction factor C(E e ) as [12,14] Using the matrix element of the effective operator Ξ JLM , we obtain and F (Z, E e ) = α 2 −1 + α 2 1 .

Allowed and first-forbidden transitions of axial vector current
We focus on the transition rate due to the space component of the axial vector current.
In the impulse approximation, the axial vector current is given as with the nucleon field operators ψ, ψ † , at position r, spin σ, and isospin τ . The transition density ρ JL (r) represented in the radial coordinate is defined as The reduced matrix element of the effective operator Ξ JLM is given in terms of the radial integral of the transition density ρ JL (r) multiplied by combination of the electron and neutrino wave functions with the coefficients c g and c f given in the Appendix C: The leading-order formula by Behrens-Bühring (LOB) of Ref. [10] conventionally used in the nuclear structure calculations can be derived by using approximate lepton wave functions in Eqs. (56) and (59). We take the LO electron wave function and the leading-order approximation of the neutrino wave function. For the allowed transition with ∆J π = 1 + , we approximate the s-wave wave functions as a constant number: and neglect all other partial waves. For the spin-dipole transition with ∆J π = 0 − , 1 − , and 2 − , in addition to the above approximation to the s-wave wave function, we use the following leading-order approximation for the p-wave wave functions 10 For the allowed ∆J π = 1 + transition, two partial waves of leptons (κ e , κ ν ) = (−1, −1) and (1, 1) contribute within the LOB, In the last step, we use the approximation for the lepton wave functions. As a result the shape correction factor is given as The second example is the first-forbidden transition ∆J π = 0 − . The leading-order partial waves are (κ e , κ ν ) = (−1, 1) and (1, −1). We then obtain In order to compare our formula with LOB, for example, Eq. (10.56) of Ref. [11], introducing nuclear matrix elements we obtain where Here α is the fine structure constant. A similar comparison can be done for the transitions to 1 − and 2 − states, and we can confirm the use of the approximate lepton wave function within our formalism leads to the 'conventional' formula of the decay rate.

Analysis with a schematic model
In the following, we examine the validity of the approximation for the electron wave function proposed in this work by using a schematic model of transition density. Three sets of treatment of the lepton wave function (i) exact, (ii) LO, and (iii) NLO are defined.
By (i) 'exact', we use the electron wave function obtained by a numerical solution of the Dirac equation and the spherical Bessel function for the neutrino wave function. In (ii) LO and (iii) NLO, we approximate the electron wave function by the LO and NLO wave functions described in the previous section. Notice that we do not expand the neutrino wave function. We use the uniform charge distribution for the nuclear charge with a charge radius R A = 1.2 × A 1/3 fm, and a transition density given by a sum of two Gaussians. The analytic expressions for the LO and NLO terms of the electron wave functions are summarized in the Appendix B. We found that the numerical results of LO are very close to those of 'conventional' formula LOB.
The explicit forms of the LO and NLO approximation of the electron wave functions of κ e = −1 with the s-wave large component (G −1 ) and the p-wave small component (F −1 ) are given as where x = r/R A and ξ = αZ/(2R A ). Figure 1 shows the electron wave functions G −1 and F −1 at E e = 10 MeV for Z = 82 and A = 208. The 'exact' and the 'LO' wave functions are shown by the solid and short-dashed curves, respectively. The deviation of the LO wave function from the exact one grows as r increases. One can see that the deviation is larger for an s-wave than a p-wave wave function. By taking into account the NLO correction, the wave functions are greatly improved, but a slight deviation from the 'exact' wave function still remains for a larger r region. For the uniform charge distribution, the Coulomb potential for r > R A agrees with the point Coulomb potential. Therefore, by connecting the NLO wave function with the combination of the analytic form of the regular and irregular point Coulomb wave functions, we can obtain the electron wave functions for r > R A with improved accuracy (NLO * ). This is indeed the case, as shown in the blue dashed curves in Fig. 1. Figure 2 is the same as Fig. 1 but for Z = 28, A = 80. One sees that the effects of the NLO correction are smaller than in the Z = 82 case, though the deviation of the LO wave function is distinct for the s-wave.
For β + decay, the sign changes for the odd power terms of ξ. The Coulomb potential enters in the Dirac equation in the form of E e − V C as given in Eqs. (18) and (19). The Coulomb effect is constructive to E e for an electron, while it is destructive for a positron. For E e > |V C |, the deviation from the LO wave function becomes smaller for a positron than for an electron.
The E e and Z dependence of the NLO correction is parametrized essentially by two non-dimensional parameters R A E e and R A ξ = αZ/2. Figure 3 shows the deviation of the approximate LO and NLO electron wave functions G −1 at the nuclear surface r = R A , , as a function of R A E e and Z. One sees a considerable deviation for a larger nuclear charge Z and a higher R A E e value. The LO approximation overestimates the amplitude of the wave function at the nuclear surface. By including the NLO correction, an error is notably reduced. For E e = 10 MeV, the case of Z = 82, A = 208 and Z = 28, The difference between the LO and 'exact' lepton wave functions observed above certainly affects the beta-decay rate. The magnitude of the effect depends on the transition density of nuclear weak currents. To examine the effects on the beta-decay rate, we take the following simple form of the transition density for the Gamow-Teller and spin-dipole transitions: Here we take r 1 = 0. transition density multiplied by r 2 , (r/R A ) 2 ρ tr (r), is shown in Fig. 4 with a as a parameter. For a = −0.2, the contribution of the r < R A region is predominant for the transition matrix element, while for a = −1, the r ∼ R A region gives a prevailing contribution to the matrix element. For a = −0.6, a strong suppression of the matrix element would take place.
In what follows, we investigate the decay rate using the transition density (81). First, we use the exact electron wave function and set Z = 82 ± 1 for β ∓ decay and A = 208. The decay rate shown in obtain a similar a dependence of the β + decay rate for the normalized rate as shown by filled symbols.
We then study the validity of the approximation of the electron wave function thoroughly. At the end of the study with the schematic model, we investigate the Z dependence of the NLO correction. Figure 7 shows the Z dependence of the β − decay rate for E 0 = 10 MeV for the transition involving heavy neutron-rich nuclei. However, it is apparent that our NLO approximation works well for a wide range of the nuclear charge.

EDF transition density and NLO electron wave function
To investigate the validity of our formalism in realistic cases, we use the transition densities microscopically calculated by a nuclear energy-density functional (EDF) method. Since the details of the formalism can be found in Ref. [44], here we recapitulate the basic equations relevant to the present study. In the framework of the nuclear EDF method we employ, the ground state of a mother nucleus is described by solving the Kohn-Sham-Bogoliubov (KSB) equation [45] h where the KS potentials h andh are given by the EDF. An explicit expression of the potentials can be found for example in the Appendix of Ref. [46]. The chemical potential λ is determined so as to give the desired nucleon number as an average value. The superscript q denotes n (neutron, τ z = 1) or p (proton, τ z = −1). The excited states |f ; J π in a daughter nucleus are described as one-phonon excitations built on the ground state |i of the mother nucleus as where a † n (a † p ) and a n (a p ) are the neutron (proton) quasiparticle (labeled by α and β) creation and annihilation operators that are defined in terms of the solutions of the KSB equation (82) with the Bogoliubov transformation. The phonon states, the amplitudes X f , Y f and the vibrational frequency ω f , are obtained in the proton-neutron quasiparticle-random-phase approximation (pnQRPA). The residual interactions entering into the pnQRPA equation are given by the EDF self-consistently. With the solutions of the pnQRPA equation, the transition density is given as in a standard quasi-boson approximation. One obtains the transition density in the radial coordinate as which is independent of K in the present case for spherical systems. Thus, the input transition density is obtained by ρ JL (r) = √ 2J + 1ρ JL0 (r).
We apply our formula for the medium-heavy Ni and Sn isotopes. Since a considerable contribution of the first-forbidden transition is predicted in the Sn isotopes [30], we take 160 Sn as an example in the present study. Furthermore, an interplay between the allowed and first-forbidden transitions has been discussed around 78 Ni [47], and we thus take 80 Ni as a target of the present study as well and employ the same Skyrme and pairing EDF as in Ref. [47]. Within the pnQRPA, the maximum electron energy is given as The transition densities ρ JL of the J π = 1 + (E 0 = 12.1 MeV) and 0 − (16.1 MeV) states in 160 Sn are shown in Fig. 8. Those states give the largest contribution to the transition rate for each J π . One sees there are nodes similarly to the transition densities of the schematic model. In such a case, the contribution around the nuclear surface r ∼ R A is important.
Using the transition densities microscopically calculated by the EDF method, we evaluate the half-life of β − decay of the allowed Gamow-Teller and the first-forbidden spin-dipole transitions of 80 Ni and 160 Sn. The β decay rates are calculated within the impulse approximation for the space component of the axial vector current only. Here we use the effective axial vector coupling constant g A = 1. The half-life is calculated using the 'exact' formula without approximation for the lepton wave functions. The contribution of all the states up to E 0 ∼ 14 MeV (16 MeV) for 80 Ni ( 160 Sn) are included. Shown in Tab. 1 is the half-life thus calculated for each J π . We show the ratios of the half-life t 1/2 (LO)/t 1/2 (exact) and t 1/2 (NLO)/t 1/2 (exact) in the table as well. As suspected the LO overestimates the transition rate by about 5 to 15% depending on the type of the transition and nuclide. Therefore, the half-lives are underestimated. The deviation from the 'exact' calculation is larger for Sn than for Ni. Introducing the NLO correction, those errors are nicely restored, as shown in the third column of Tab. 1. We can therefore argue that our simple NLO formula is very effective in realistic calculations.  Table 1 Half life t 1/2 of the β − decay of 80 Ni and 160 Sn to the daughter nucleus with the states J π . The ratios of the half life t 1/2 (LO)/t 1/2 (exact) and t 1/2 (NLO)/t 1/2 (exact) are shown in the second and the third column. We have investigated the Coulomb effects on the beta-decay rate. The decay rate is determined by the product of the lepton and hadron current densities. A widely used formula relies on the fact that the low-energy lepton wave functions in a nucleus can be well approximated by a constant and linear to the radius for the s-wave and p-wave wave functions, respectively. We found, however, the Coulomb wave function is conspicuously different from such a simple approximation for heavy nuclei with large Z by numerically solving the Dirac equation. We then have proposed formulas of the nuclear beta-decay rate that are useful in a practical calculation.

20
In our proposed formulas, the neutrino wave function is treated exactly as a plane wave, while the electron wave function is obtained by iteratively solving the integral equation; thus, we can control the uncertainty of the approximate electron wave function order by order. The leading-order approximation gives a formula that is almost equivalent to the widely used one and overestimates the decay rate by about 50-100% for heavy nuclei with Z ∼ 80. We demonstrated that the next-to-leading-order formula reproduces well the exact result for a schematic transition density as well as a microscopic one obtained by a nuclear energydensity functional method. For the beta decay involving heavy neutron-rich nuclei, the NLO will be needed for the determination of the Gamow-Teller strength from the beta-decay rate.
We considered only the space component of the axial vector currents and kept only the lowest multipoles. The time components as well as the vector currents can have a comparable contribution to the decay rate, and we plan to present these improvements in a sequel to the present article. The beta decay provides a unique spectroscopic tool of exotic nuclei, that is, the angular correlation contains rich information of nuclear structure. Furthermore, the electron/muon capture is an important process in the application to astrophysics and fundamental physics. It is straightforward to extend our formalism towards these directions.

Acknowledgment
We would like to thank Prof. K. Koshigiri for useful discussions. This work was in part supported by the JSPS KAKENHI Grants Nos. JP18H01210, JP18H04569, JP18K03635, JP19H05104, JP19H05140, and JP19K03824, the Collaborative Research Program 2019-2020, Information Initiative Center, Hokkaido University, and the JSPS/NRF/NSFC A3 Foresight Program "Nuclear Physics in the 21st Century." The nuclear EDF calculation was performed on CRAY XC40 at the Yukawa Institute for Theoretical Physics, Kyoto University.

A Derivation of the decay rate
In this appendix, we show the derivation of the decay rate Eq. (16) for β − decay. In the present case, it is useful to expand the effective Hamiltonian in terms of the angular momentum. With the partial wave expansion of the neutrino wave function, we get where we use the abbreviated notation (upper sign) defined in Eq. (14). Since the neutrino mass is negligible, we used here g κ (r) = −S κ f −κ (r) and f κ (r) = S κ g −κ (r). We also have an 21 alternative expression: Thus, we obtain two equivalent expressions: where we used Eq. (13). As in Eq. (2), the hadron current is composed of the vector and axial vector components.
We can thus write The products of the two-component spinors are given as where ǫ is a unit vector. Using these relations, we obtain the effective Hamiltonian Eqs. (11) and (12). According to the Wigner-Eckart theorem, the M-dependence of the spherical tensor Ξ JLM , Eq. (12), is known as where the reduced matrix element f Ξ JL (κ e , κ ν ) i is independent of M, and J f is the angular momentum of the final nuclear state.
With the obtained Hamiltonian H eff , the decay rate is given by where J i is the angular momentum of the initial nuclear state. For an arbitrary function X (κ, µ), we have

B Explicit formula for the uniform charge distribution
For the uniform charge distribution of nuclei with radius R A , the Coulomb potential for an electron is given as where α is the fine structure constant and x = r/R A . The electron wave functions D (i) , d (i) , H (2) , and h (2) are given as where ξ = αZ/(2R A ) for an electron and ξ = −αZ/(2R A ) for a positron. The functions s a , t a , w a , y a , and z a for k = 1 are given as Similar formulas for k = 2 are given as  Tables C1 and C2 list the explicit numbers of the coefficients of each (κ e , κ ν ) for the Gamow-Teller and spin-dipole transitions, respectively.

D Tables of electron and positron wave functions
We provide numerical tables of the four constants α 1 , α −1 , α 2 , and α −2 needed to construct the electron and positron wave functions in this paper. Since these constants are strongly dependent on the electron momentum p e and charge number Z of a nucleus, we rewrite them to L 0 , λ 2 , µ 1 , and µ 2 according to Ref. [12]: where F 0 (Z, E) = 4(2p e R A ) −2(1−γ 1 ) e πν Γ(γ 1 + iν) Γ(2γ 1 + 1) for the uniform nuclear charge distribution with radius R A . Actually, these variables have a milder momentum and charge dependence than α ±1,2 . First, we calculate α ±1,2 by numerically solving the Dirac equation and convert them into L 0 , λ 2 , µ 1 , and µ 2 at various p e /m e and Z. These generated tables are respectively interpolated by assuming the following polynomial function at three regions, p e /m e =0.01-1, 1-10, and 10-100: The weights of the polynomial W st are determined by the least-square method. Finally, we reconstruct α ±1,2 from L 0 , λ 2 , µ 1 , and µ 2 . Since all α ±1,2 values are positive, the reconstruction can be made easily. We confirm that the resulting numerical tables are accurate more than 3-4 digits with (m, n) = (3,4). For the convenience of a user, we provide a FORTRAN program code to generate α ±1,2 with a given p e /m e (=0.01-100) and Z(= 1-90) for β ∓ decay as supplemental material. The nuclear charge radius R A is set to be 1.2A 1/3 fm with the nuclear mass number A. To cover stable and neutron-rich unstable nuclei for the β − decay, a variation of the charge radius can be considered among five options: A = 2Z, 2.5Z, 3Z, 3.5Z, and 4Z.