Electron wave functions in 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 similarly in Refs. [34,35] for the semi-leptonic nuclear processes and electron scattering [36]. It is thus straightforward to apply it to the charged-current neutrino reaction and lepton capture reaction. Developing an analytic formula of beta decay based on Ref. [15] would also contribute to a precise understanding of the neutrino-nucleus reactions to extract neutrino properties from neutrino experiments as discussed in Ref. [37,38].
The high-energy forbidden transitions occur under the exotic environment with high Q β value [39]. Therefore, in this work we revisit the formulation of beta decay for not only the allowed but the first-forbidden transitions induced by the Gamow-Teller and spin-dipole type operators. We provide a simple way to improve the widely used formula in the nuclear beta-decay study to apply for nuclei with large Z and away from the stability line. We start from the formulation of Koshigiri et al. [15] and use iterative solutions of the integral equation [10,17]. In the previous formalism, one often expands the electron and neutrino wave functions in the long wave-length approximation and collect terms in a systematic way.
Here we avoid this expansion of the neutrino wave function. We use an analytic form of the LO and NLO electron wave functions combined with the numerical table of the electron wave function at the origin. This makes the formula of the beta-decay rate simple and easy to use. This paper is organized in the following way. We start from the formulation of the beta decay with the partial wave expansion for the lepton wave functions in Sect. 2. We provide an explicit expression of the first (LO) and the second (NLO) iteration of the integral equation for an electron wave function in Sect. 3. Formulas of the beta-decay rate are given and compared with the widely used one in Sect. 4. The formulas of LO and NLO are examined in Sect. 5, using a schematic transition density that is given by a sum of two Gaussians. We then in Sect. 6 apply the formulas to the neutron-rich Ni and Sn isotopes where the transition densities are microscopically obtained by a nuclear energy-density functional (EDF) method.
Summary and perspectives are given finally in Sect. 7.

Effective Hamiltonian
An effective Hamiltonian for a low-energy charged-current reaction is given as where (x) represents either the electron, muon, or tau field andψ = ψ † γ 0 . The hadron current J µ (x) is given by the vector and axial vector currents where G F = 1.166 × 10 −5 GeV −2 is the Fermi coupling constant and V ud = 0.9737 is the CKM matrix [40]. Here we take natural units = c = 1.
The effective Hamiltonian describes semi-leptonic nuclear weak processes such as lepton capture, neutrino reaction, and β ± decay. For β − decay, i → e − (p e ) +ν e (p ν ) + f , where i and f respectively denote the initial and final nuclear states, and p is the lepton momentum, the transition matrix element is given as and for β + decay where u and v are the Dirac spinors of the neutrino and antineutrino, respectively. The electron scattering wave functions ψ with the superscript (−) and (+) satisfy the incoming and outgoing boundary conditions, respectively.
The normalization of the scattering wave function in the plane wave expansion ψ e (x) → u(p e ) exp(ip e · x) is given as It is noticed the electron wave functions (G κ , F κ ) in Ref. [15] are defined by multiplying e i∆κ e to ours, while those of Refs. [10,16] are given by multiplying √ 2p e to ours.
The neutrino and antineutrino wave functions are respectively expanded as where j l (x) is the spherical Bessel function of order l, S κ = sgn(κ) is the sign of κ, and With the partial wave expansion of the electron and neutrino wave functions, one obtains the following form [15]: and Ξ JLM (κ e , κ ν ) = S κe dr where [O k 1 ⊗ O k 2 ] k 3 m 3 denotes the tensor product. Here we adopt the following simplified notation for el and neu : el = κe,µe i −lκ e e i∆κ e (l κe , m e , 1/2, s e |j κe , µ e )Y lκ e me (p e ), 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 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 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 where The explicit expressions of H k , D 4 Decay rate and comparison with the conventional formula

Decay rate
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 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 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 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.9R A , b = R A /4 and r 2 = 3r 1 /4. By varying −1 ≤ a ≤ −0.2, we investigate the validity of the approximation for the electron wave function on the decay rate. The 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 examine the validity of the LO and NLO approximation of the electron wave function. The decay rate is studied for Z = 82, A = 208 with E 0 = 10 MeV using the transition density (81). In Fig. 5, the decay rate evaluated using the 'exact' electron wave function is shown by the solid curve with ∆a = 0.025 in arbitrary unit normalized to unity at  where the contribution at r ∼ R A is more important than r < R A . The suppression takes place at larger a for LO than the exact calculation. The use of the augmented NLO electron wave function significantly improves for a < −0.7, and works reasonably well even when a severe cancellation between the inner and the outer contribution of the integration takes place around a = −0.6. We obtain a similar a dependence of the β + decay rate shown in Fig. 6.
For β + decay, LO approximation gives reasonable description for Z = 82 and E 0 = 10MeV. 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 = suggests that the transition probabolity extracted from the beta-decay rate using the LO can be underestimated 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 q (rσ) − λ qhq (rσ) 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 evalu- 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.

Summary
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.
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 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 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 r R d 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

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 Table C2 Same as Tab. C1 but for the spin-dipole transition (L = 1, ∆J π = 0 − , 1 − , 2 − ). 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.