Three-Dimensional Approach Applied to Quasi-stationary States of Deformed $\alpha$-Emitters

The quantum tunneling problem of $\alpha$-cluster decay provides insight into not only the quasi-stationary states of unstable nuclei but also the nature of $\alpha$-clustering in surrounding nuclear medium DLW ; Delion09 ; Buck ; cxu06 ; Mohr06 ; review1 ; review2 ; Ismail16 ; Perez19 ; Denisov09 ; Royer ; Poenaru . Much of our understanding of $\alpha$-cluster decay comes from the well-established quantum picture where an $\alpha$-cluster tunnels through the Coulomb barrier after its formation on the surface of nuclei. The $\alpha$-clustering process is rather complex to handle, which in principle involves a quantum four-body problem embedded in nuclear medium THSR ; quartet2017 ; Po ; Yasu21 ; Eiji20 ; SY20 ; SY21 ; Xu2017 ; Xu2016 . In contrast, the tunneling probability can be quantitatively estimated from the WKB approximation by assuming an $\alpha$-cluster interacting with a spherical daughter nucleus. This assumption is appropriate for $\alpha$-emitters with small deformations Buck ; cxu06 ; Mohr06 . However, large deformations could be involved in the $\alpha$-cluster decay of heavy nuclei and exotic nuclei far away from the $\beta$-stability line Oganessian ; Batchelder ; Te2019 ; Zhang . The exact 3D decay width of well-deformed $\alpha$-emitters is difficult to obtain by treating self-consistently both deformation and mixing of angular momenta. An empirical way to include effect of large deformations is to average the tunneling probabilities of all angles, which yields an enhanced $\alpha$-cluster tunneling probability. For $\alpha$-transitions with measured core excitation energies, the total wave function of the system can be expanded in terms of the ground and excited states channels and the coupled Schrödinger equations can be solved Barmore ; channel1 ; Delion06 .

The exact solution of quasi-stationary state with an anisotropic potential is quite interesting in many quantum systems QT ; Spins ; Electrons ; QD ; Superconductor ; TPA ; MTPA . In this work, we apply a 3D approach to investigate the quasi-stationary state of $\alpha$-cluster decay problem, namely the 3D two-potential approach (3D-TPA) MDTPA . We emphasize that the quasi-stationary state of $\alpha$-cluster decay highly resembles a bound state more than a scattering state. The anisotropic $\alpha$-core potential is properly divided into inner and outer potentials on a separation surface inside the Coulomb barrier. The formed $\alpha$-cluster is considered to stay initially in the bound state generated by the inner potential, and then transforms to a quasi-stationary state by switching on the outer potential. Using the time-dependent perturbation theory, both the decay width and the energy shift of $\alpha$-cluster decay can be well defined by inner 3D bound state wave function and outer scattering state wave function on the separation surface. The multi-dimensional quantum tunneling problem for well-deformed $\alpha$-emitters is reduced to a problem of true bound state plus non-resonance scattering state MDTPA . It is worth noting that the choice of separation surface does not affect the final results as long as it is inside the classical forbidden region. The challenge here is the numerical solution of the inner 3D bound state wave function on the separation surface, whose value is approximately of the order of 10^-15. Although several numerical methods can be applied to solve the 3D Schrödinger equation to obtain the bound state wave function such as the grid-based approach, the imaginary time propagation method, and the basis expansion method, however, it is still an open question to obtain accurately the inner 3D bound state wave function at large distances due to the limitation of matrix size or the number of bases Imaginary ; Gaussian . In this sense, the anisotropic harmonic oscillator potential with exact solutions is a good choice to simulate the internal $\alpha$-core 3D potential. The 1D and 3D decay widths are compared for several typical non-spherical $\alpha$-emitters with large quadrupole and hexadecapole deformations. Note that the formation process of the $\alpha$-cluster on nuclear surface is not touched here, which in principle does not affect the comparison between 1D and 3D decay widths.

The paper is organized as follows: In Sec.II, we give the formalism of 3D two-potential approach for the multi-dimensional quantum tunneling problem. The anisotropic $\alpha$-core effective potential and the choice of separation surface are introduced in Sec.III. In Sec.IV, both the inner 3D bound state wave function and outer scattering state wave function are explicitly given. In Sec.V, the 3D approach is applied to several deformed $\alpha$-emitters ${}^{236}\textrm{U}$, ${}^{244}\textrm{Pu}$, ${}^{246}\textrm{Cm}$, ${}^{250}\textrm{Cf}$, and ${}^{254}\textrm{Fm}$. The possible theoretical uncertainties are also discussed. Conclusions are presented in Sec.VI.

The time-dependent perturbation theory is applied to solve the multi-dimensional tunneling problems in 3D-TPA MDTPA . The main integrant is the separation of the 3D $\alpha$-core potential $V({\mbox{\boldmath$r$}})$ into the inner $U({\mbox{\boldmath$r$}})$ and outer $W({\mbox{\boldmath$r$}})$ potentials on the separation surface $\mathcal{S}$

The inner potential is

where $U_{0}$ is the minimal value of $V({\mbox{\boldmath$r$}})$ on the separation surface $\mathcal{S}$. The outer potential is

The quasi-stationary state of $\alpha$-cluster decay can be described by the time-dependent Schrödinger equation,

where $H_{0}=-\frac{\hbar^{2}}{2\mu}{\mbox{\boldmath$\nabla$}}^{2}+U({\mbox{\boldmath$r$}})$ is the Hamiltonian of bound state $|\Phi_{i}\rangle$ confined in the inner potential $U({\mbox{\boldmath$r$}})$, and the corresponding time-independent Schrödinger equation is $[-\frac{\hbar^{2}}{2\mu}{\mbox{\boldmath$\nabla$}}^{2}+U({\mbox{\boldmath$r$}})]|\Phi_{i}\rangle=E_{0}|\Phi_{i}\rangle$. The separation surface $\mathcal{S}$ is not necessarily spherical, as long as it is taken between the equipotential surfaces $S_{1}$ and $S_{2}$. At $t>0$, the “unperturbed” bound state $|\Phi_{i}\rangle$ is no longer an eigenstate of the total Hamiltonian $H=H_{0}+W({\mbox{\boldmath$r$}})$, but a wave packet spreading in time due to the perturbation $W({\mbox{\boldmath$r$}})$

where $b_{0}(t)$ and $b_{\textbf{k}}(t)$ are the probability amplitudes of finding the system in the eigenstates $|\Phi_{i}\rangle$ and $|\varphi_{\textbf{k}}\rangle$, respectively. The amplitudes $b_{0}(t)$ and $b_{\textbf{k}}(t)$ can be found from Eq.(4) with the initial condition: $b_{0}(t)=1,b_{\textbf{k}}(t)=0$. The energy shift $Re(\epsilon_{0})$ and the width $\Gamma=-2Im(\epsilon_{0})$ of the quasi-stationary state are directly related to the pole in the complex E-plane using the Green’s function technique TPA

The Green’s function $\tilde{G}$ is given by

where $\tilde{W}=W+U_{0}$ is used instead of $W$, to ensure the potential vanishes for $r\rightarrow\infty$. $G_{0}(E)$ is given by

Note that the above derivations are general, but the numerical solution of $\tilde{G}$ is rather difficult and converges very slowly. To make it feasible within the capacity of computer calculation, $\tilde{G}$ is expanded in powers of $G_{\tilde{W}}$, namely the Green function corresponding to $\tilde{W}({\mbox{\boldmath$r$}})$

Substitute Eq.(9) into Eq.(6), then retain the first order and assume that the energy shift is small compared to $E_{0}$, one can replace $G_{\tilde{W}}(E)$ by $G_{\tilde{W}}(E_{0})$

The Schrödinger equation of the Green’s function is

and the spectral representation for the Green’s function is

where $E_{\textbf{k}}=\hbar^{2}\textbf{k}^{2}/2\mu$. The $|\varphi_{\textbf{k}}\rangle$ is the eigenstate corresponding to non-resonance scattering states, which satisfies the following Schrödinger equation

One can obtain the total width $\Gamma$ as an integral over the partial width $\Gamma_{\textbf{k}}$ MDTPA

where $k_{0}=\sqrt{2mE_{0}/\hbar^{2}}$. The partial width $\Gamma_{\textbf{k}}$ is

where the symbol $\stackrel{{\scriptstyle\leftrightarrow}}{{{\mbox{\boldmath$\nabla$}}}}_{\textbf{r}}$ means the gradient on the right minus the gradient on the left. ${\mbox{\boldmath$k$}}=\{k_{1},k_{2},k_{3}\}$ is the momentum vector. We note that Eq.(15) is similar to the Bardeen formula for the tunneling coupling between adjoining wells from a many-particle point of view, which is widely used in the solid state and atomic physics Bardeen ; Gurvitz . Finally, the decay half-life is given by

where P_α is the formation probability of $\alpha$-cluster on nuclear surface. Recent microscopic calculation of the $\alpha$-cluster formation probability for the ideal Po isotopes has been performed by the quartetting wave function approach (QWFA), in which the intrinsic motion of the four nucleons forming the $\alpha$-cluster and the center of mass motion between the $\alpha$-cluster and the core are correctly treated SY20 ; SY21 ; Xu2016 ; Xu2017 . Empirically, the $\alpha$-cluster formation probability $P_{\alpha}$ is known to change abruptly across the major shell closures and differs for even-even, odd-A and odd-odd nuclei, as indicated by the experimental systematics. The decay width $\Gamma$ can be directly obtained by using Eq.(14) and Eq.(15), however, it is more convenient to re-formulate $\Gamma$ in the following way

where $R_{s}(\theta,\phi)$ is the radius of the separation surface $\mathcal{S}$. For spherical emitters, Eq.(LABEL:gam3) reduces exactly to the 1D decay width $\Gamma={\hbar^{2}\over\mu k_{0}}|\phi_{i}(R)\chi_{k}^{\prime}(R)-\chi_{k}(R)\phi_{i}^{\prime}(R)|^{2}$ TPA .

The 3D-TPA requires the information on 3D bound-state wave function deep inside the classical forbidden region where the wave function decreases exponentially and is extremely small (of the order of  10^-15). This poses a big challenge for numerical calculations, even with high-performance parallel computing. Here we simulate the 3D $\alpha$-core potential with an anisotropic harmonic oscillator potential that is of rotational symmetry. The relevant Hamiltonian is of Nilsson-form, which has a $L^{2}$ correction term Nilsson ; SPM

with

where $\mathcal{R_{s}}(\theta)$ describes the boundary condition $V_{I}(r,\theta)-\kappa\omega_{0}\hbar\big{[}L^{2}-\langle L^{2}\rangle_{N}\big{]}=V_{C}(r,\theta)$. The anisotropic harmonic oscillator potential with depth parameter $D$ is

where $\omega$, $\omega_{0}$ and $\gamma$ read

The deformed Coulomb potential is

The half-density radius $R(\theta)$ is given by

where the parameter $R_{0}=1.07A_{d}^{1/3}$rrr , and $\beta_{2}$ and $\beta_{4}$ are quadrupole and hexadecapole deformations, respectively. The deformed Coulomb potential $V_{C}$ can be regarded as an isotropic potential at very large distances. The $L^{2}$ correction term is introduced to eliminate the angular momentum degeneracy, resulting in different depths of potential for different angular momenta.

The details of the $\alpha$-core 3D potential $V_{0}(r,\theta)$ are shown in Fig.1, in which the inner potential joins with the outer Coulomb potential at $\mathcal{R_{s}}(\theta)$. The total $\alpha$-core potential is divided into two parts by the separation surface $\mathcal{S}$, which can be chosen between the equipotential surfaces $S_{1}$ and $S_{2}$ with the same energy $E_{0}$. Here, we take the separation surface $\mathcal{S}$ as coinciding with $\mathcal{R_{s}}(\theta)$ (see Fig.2). In this way, one can minimize the numerical errors associated with $\Phi_{i}(r,\theta,\phi)$ and $\varphi_{\textbf{k}}(r,\theta,\phi)$. Note that the boundary $\mathcal{R_{s}}(\theta)$ is azimuth-independent ($\phi$-independent). Thus, the separation surface $\mathcal{S}$ is also azimuth-independent in our calculations.

We use the analytical 3D wave function of anisotropic harmonic oscillator potential as an approximation of $\Phi_{i}(r,\theta,\phi)$, which is the eigenstate of $E_{0}$ corresponding to the inner potential. As a matter of fact, we only need the information of $\Phi_{i}(r,\theta,\phi)$ on the surface $\mathcal{S}$ where the 3D wave function of anisotropic harmonic oscillator is considered to be a good approximation for well-deformed $\alpha$-emitters

where

and

with parameters $\beta_{\perp}$, $\beta_{z}$, $N_{n_{\rho}}^{m}$ and $N_{n_{z}}$

The expectation value $E_{0}$ is given by SPM

where the depth $D$ can be obtained by matching the value of experimental decay energy $Q_{\alpha}$ with $E_{0}$. The choice of $n_{\rho}$ and $n_{z}$ should fulfill the so-called Wildermuth-Tang rule WT

where $G$ is the global quantum number and $g_{i}$ are the oscillator quantum numbers of the nucleons forming the $\alpha$-cluster. The $G$ value is restricted by the Pauli principle. For instance, $G$ is usually taken as $G=22$ for heavy nuclei with neutron numbers $N>126$. In principle, $n_{z}$ can be taken as an even number from 0 to $G$. Here, $n_{z}$ is taken as $n_{z}=10$ as an example.

The scattering wave function $\varphi_{k}(r,\theta,\phi)$ is solved from the scattering Schrödinger equation Eq.(13) with the outer potential (see Fig.3(b)), which has the form scattering ; QM

where $\chi(r,\theta)$ can be represented by the linear combination of the regular $F_{l}(r)$ and irregular $G_{l}(r)$ solutions of Coulomb potential. The corresponding coefficients in the linear combination of $F_{l}(r)$ and $G_{l}(r)$ are directly related to the non-resonant scattering phase shift for the outer potential $\tilde{W}(r)$ MTPA

where the relative phase $\delta_{l}$ obeys

and

The inner wave function $\Phi_{i}(r,\theta,\phi)$ and the outer Coulomb function $\chi(r,\theta)$ are shown in Fig.3 and Fig.4, respectively. For spherical emitters, the inner wave function shown in Fig.3(a) is the spherical harmonic oscillator wave function. For well-deformed emitters, the inner wave function is approximated by the anisotropic harmonic oscillator wave function as shown in Fig.3(b). As for the scattering state, the wave function $\chi(r,\theta)$ becomes almost isotropic at large distances (see Fig.4).

As shown in Table 1, all the $\alpha$-cluster emitters selected in calculations are well-deformed with $\beta_{2}>0.22$ and far away from the major shell closures $Z=82$ and $N=126$. Moreover, only the $\alpha$-transitions between ground-states ($0^{+}\longrightarrow 0^{+}$) are considered in order to check the validity of 3D-TPA. There are three adjustable parameters in 3D-TPA that should be determined, namely $a_{r}$, $\kappa$, and $D$. $a_{r}$ is associated with the angular frequency $\omega$ of harmonic oscillator. $\kappa$ is the strength coefficient of orbit correction term. $D$ is the depth of harmonic oscillator potential, which can be obtained by matching the experimental decay energy $Q_{\alpha}$ with the energy $E_{0}$ once the quantum numbers $G,n_{z},m$ are determined. It is noted that different quantum numbers $G,n_{z},m$ can be chosen. Here the same parameters and quantum numbers are used for all emitters, namely, $a_{r}=0.7883,\kappa=0.0054,G=22,n_{z}=10$, and $m=0$.

We show in Table 1 the comparison between 1D and 3D $\alpha$-decay half-lives for typical $\alpha$-emitter ${}^{236}\textrm{U}$, ${}^{244}\textrm{Pu}$, ${}^{246}\textrm{Cm}$, ${}^{250}\textrm{Cf}$, and ${}^{254}\textrm{Fm}$. Note that the quadrupole and hexadecapole deformations are taken from Ref.bb24 . $Q_{\alpha}$ and experimental results of $\alpha$-decay half-lives are taken from Refs.AME1 ; AME2 . One can see from Table 1 that the experimental $\alpha$-decay half-lives vary in a quite large range from $10^{4}$ to $10^{15}$ seconds. This is also helpful for testing the validity of present 3D model for not only short-lived but also long-lived $\alpha$-emitters. The formation problem of $\alpha$-cluster on the surface of parent nucleus is not considered here and its value is assumed to be unity (P_α=1). This assumption is reasonable because the $\alpha$-cluster formation probability is known to change dramatically in the vicinity of shell closures but smoothly in the open-shell region here. One can see from Table 1 that deviations of the calculated $\alpha$-decay half-lives exist between the 1D case and 3D case. For all the $\alpha$-emitters considered in this work, the largest deviation occurs for the decay of ²³⁶U, which is possibly due to its small decay energy. In general, the 1D results are all reduced by taking nuclear deformation into account (shown also in Fig.5). Similar conclusions can also be found in the empirical approaches cxu06 . More importantly, the 3D-TPA reproduces nicely the systematics of experimental data with only one set of parameters. This is quite different from the 1D-TPA without the consideration of deformation, in which one may need to adjust the parameters such as the depth of potential for each emitter.

Finally, we discuss the possible theoretical uncertainties of present 3D-TPA calculations, which mainly come from the following aspects: a) the error associated with the approximations such as the replacement of $G_{\tilde{W}}(E)$ by $G_{\tilde{W}}(E_{0})$ in 3D-TPA. The correction due to this approximation is expected to be negligible because of very small energy shift from $E_{0}$ to $E$; b) the inner 3D wave function is approximated by the exact solutions of an anisotropic harmonic oscillator potential. We solved numerically the bound-state wave function corresponding to the inner potential in the spherical case, and found its value only deviates from the exact result by several percentages on the separation surface; c) the scattering wave function $\chi(r,\theta)$ is considered to be isotropic on the separation surface. We have checked this approximation by taking the decay of ${}^{236}U$ as an example, and found that the final results are almost not affected.

Large deformation is relevant to the $\alpha$-cluster decay of heavy nuclei and nuclides far away from the $\beta$-stability line. In this work, we apply a three-dimensional approach with the Nilsson-form Hamiltonian to calculate the decay widths of typical $\alpha$-emitters ${}^{236}\textrm{U}$, ${}^{244}\textrm{Pu}$, ${}^{246}\textrm{Cm}$, ${}^{250}\textrm{Cf}$, and ${}^{254}\textrm{Fm}$ by dividing the 3D effective potential into a bound-state inner region and a scattering outer region. The inner wave function can be well approximated by the exact solution of the anisotropic harmonic oscillator potential, and the scattering wave function can be safely considered as isotropic on the separation surface. Substantial difference is found between 1D and 3D decay width for favored transitions of these $\alpha$-emitters. The systematics of experimental $\alpha$-cluster decay half-lives is nicely reproduced. In the future, state-of-art numerical approaches can be applied to evaluate accurately the inner 3D wave function at large distances for arbitrary 3D potentials. Moreover, the combination of the present 3D-TPA and approaches of clustering such as quartetting wave function approach should be performed in order to predict reliably the $\alpha$-decay half-lives of unknown nuclei far away from the major shell closures.

$\bf{Acknowledgments}$

The work is supported by the National Natural Science Foundation of China (Grant No. 12275129) and the Fundamental Research Funds for the Central Universities (Grant No. 020414380209).