-
Nuclear stability is a cornerstone research theme in nuclear physics, governed by nucleon-nucleon interactions and manifested through critical observables including binding energy systematics, magic number configurations, decay modes, and half-life patterns. This fundamental property not only dictates the synthesis pathways of new isotopes and superheavy elements (SHEs) but also serves as a stringent testbed for nuclear models. The nuclear shell model initially proposed in the 1960s predicted post-
208 Pb magic numbers atZ=114 andN=184 [1−3], suggesting an "island of stability" centered at (Z=114 ,N=184 ) where superheavy nuclei might exhibit enhanced stability. However, modern self-consistent mean-field calculations reveal competing predictions, proton shell closures atZ=120 or 126 and neutron shell closures spanningN=172−184 [4, 5]. These theoretical divergences highlight the critical role of SHE stability studies in validating nuclear many-body approaches. Through fusion-evaporation reactions, superheavy nuclei with104≤Z≤118 have already been successfully synthesized in experiments [6−12]. The heaviest nucleus294 Og that has been synthesized so far has a half-life on the order of milliseconds and is neutron-deficient due to the limitations in available projectile-target combinations and experimental facilities. The determination of nuclear stability exhibits mass-dependent characteristics. For light and medium nuclei (A<150 ), β-decay predominance makes binding energy minimization along isobaric chains the key stability indicator. For heavy and super-heavy nuclei (A>200 ), multi-mode decay competition disrupts simple energy-based predictions. For example, although the ground state energy of214 Po is smaller than that of214 Pb by 2.7 MeV, the half-life of214 Pb is much larger than that of214 Po, due to the different decay modes:214 Pb via β-decay in contrast to214 Po via α-decay. Therefore, new stability metrics beyond traditional binding energy considerations should be developed, particularly for SHEs. The enduring mystery surrounding the existence and location of the "island of stability" continues to drive interdisciplinary efforts combining nuclear experiments, astrophysical observations, and exascale computational modeling.The binding energy of atomic nuclei predominantly arises from the dynamic equilibrium between the nuclear strong force and Coulomb repulsion [13], which governs their β-decay stability. Quantitatively, a net energy difference is manifested between the absolute values of the attractive nuclear potential energy U and the repulsive kinetic energy T. Notably, the
T/U ratio constitutes a fundamental parameter in nuclear astrophysics [14, 15], as it correlates with two critical phenomena: nuclear matter compressibility and stability mechanisms. The ratioT/U provides insights into the stiffness of nuclear matter, a key determinant of neutron star interior structure and equation of state. In neutron stars, gravitational binding stability requires the dominance of gravitational potential energyUgrav over rotational kinetic energyTrot . A significant imbalance (|Ugrav|≫Trot ) may trigger catastrophic collapse into black holes. For atomic nuclei, theT/U ratio reflects the interplay between nucleon-nucleon interaction strength (governed by nuclear force saturation) and nucleon motion intensity (driven by Fermi energy). HigherT/U values typically indicate reduced binding coherence, thereby influencing fission barriers and decay modes. This intrinsic connection between theT/U ratio and multi-scale stability (from femtometer-scale nuclei to kilometer-scale neutron stars) underscores its significance in unified nuclear many-body theories. Systematic investigations of this ratio across nuclear chart regions can unveil universal constraints on hadronic matter under extreme conditions.The Skyrme Hartree-Fock-Bogoliubov (HFB) model is a versatile microscopic theoretical framework widely used in nuclear physics research, allowing the nuclear structure at ground state such as binding energies, charge radii, deformation and shape coexistence to be investigated self-consistently [16−20]. In this study, we performed systematic calculations for 2318 even-even nuclei using the Skyrme HFB model to analyze nuclear stability and the energy components of nuclei.
-
In HFB calculations, the Hamiltonian of a system with interacting fermions is written in terms of the annihilation and creation operators
(c,c†) as [16, 17]:H=∑μμ′Tμμ′c†μc′μ+14∑μνμ′ν′˜Vμνμ′ν′c†μc†νcν′cμ′,
(1) in which
Tμμ′=−ℏ22m⟨μ|∇2|μ′⟩ and˜Vμνμ′ν′=⟨μν|V(r12)|μ′ν′−ν′μ′⟩ are the single-particle kinetic energy matrix elements and anti-symmetrized two-body interaction matrix elements, respectively. In the HFB model, the particle operators are transformed into quasiparticle operators through the Bogoliubov transformation to handle the pairing correlationαμ=∑ν(U∗νμcν+V∗νμc†ν),α†μ=∑ν(Vνμcν+Uνμc†ν),
(2) where U and V are the transformation matrices. The quasiparticle vacuum state is defined as the HFB ground state
|Φ⟩ , satisfying the conditionαμ|Φ⟩=0 . By applying Wick's theorem, the expectation value of the Hamiltonian H can be expressed as a function of the Hermitian density matrix ρ and pairing tensor κ [16]E[ρ,κ]=⟨Φ|H|Φ⟩,
(3) where
ρμμ′=⟨Φ|c†μ′cμ|Φ⟩ andκμμ′=⟨Φ|cμ′cμ|Φ⟩ .In coordinate space, the operators
crσq andc†rσq refer to annihilation and creation of nucleons, respectively, at pointr with spinσ=±1/2 and isospinq=±1/2 , which using the quasiparticle operators are expressed ascrσq=∑ν[Uν(r,σ,q)αν+V∗ν(r,σ,q)α†ν],c†rσq=∑ν[Vν(r,σ,q)αν+U∗ν(r,σ,q)α†ν].
(4) The density matrix and the pairing tensor are expressed as
ρ(rσq,r′σ′q′)=⟨Φ|c†r′σ′q′crσq|Φ⟩=∑νV∗ν(r,σ,q)Vν(r′,σ′,q′),κ(rσq,r′σ′q′)=⟨Φ|cr′σ′q′crσq|Φ⟩=∑νV∗ν(r,σ,q)Uν(r′,σ′,q′).
(5) where κ can be replaced by the pairing density matrix via
˜ρ(rσq,r′σ′q′)=−2σ′κ(rσq,r′−σ′q′) for convenience to describe time-even quasiparticle states when both ρ and˜ρ are Hermitian and time-even [16].With the Skyrme force characterized by zero-range and non-local interaction, the HFB energy in Eq. (3) can be expressed as the volume integral of the energy density
E[ρ,˜ρ]=∫d3rH(r),
(6) H(r)=HS(r)+˜HP(r)+HC(r) includes the energy densities of Skyrme interactionHS(r) , pairing˜HP(r) , and Coulomb interactionHC(r) for protons. The mean field with Skyrme forceHS(r) is a functional of the particle densityρ=∑qρq , kinetic-energy densityτ=∑qτq , and spin-current tensorJij=∑qJq,ij , with indexq=p,n corresponding to proton and neutron, respectively:HS(r)=ℏ22mτ+12t0[(1+12x0)ρ2−(12+x0)∑qρ2q]+12t1[(1+12x1)ρ(τ−34Δρ)−(12+x1)∑qρq(τq−34Δρq)]+12t2[(1+12x2)ρ(τ+14Δρ)−(12+x2)∑qρq(τq+14Δρq)]+112t3ρα[(1+12x3)ρ2−(x3+12)∑qρ2q]−18(t1x1+t2x2)∑ijJ2ij+18(t1−t2)∑q,ijJ2q,ij−12W0∑ijkεijk[ρ∇kJij+∑qρq∇kJq,ij],
(7) where
εijk is the Levi-Civita symbol withi,j,k= (1, 2, 3), and(x0,x1,x2,x3,t0,t1,t2,t3,W0,α) are the parameters of the Skyrme force. The pairing energy density is expressed as˜HP(r)=12V0[1−12ρρ0]∑q˜ρ2q
(8) with
V0 the pairing strength andρ0 the saturation density determined by the Skyrme parameters. The Coulomb interaction should be included into the energy densities for the case of protonHC(r)=e22ρp∫d3r′ρp(r′)|r−r′|−34e2(3π)1/3ρ4/3p.
(9) The Skyrme HFB equations can be obtained by varying the HFB energy in Eq. (6) with respect to ρ and
˜ρ . Detailed expressions are available in the literature [16−18]. After obtaining self-consistent solutions to the ground-state density matrix ρ and pairing density matrix˜ρ using the Skyme HFB equations, the kinetic energy can be expressed as:T=ℏ22m∫τd3r
(10) and total energy
Etot in the ground state can be obtained. The contribution of the non-local term of the Skyrme interaction to the kinetic energy is incorporated via solving the density matrix under single-particle potential. The potential energy U is calculated by subtracting the kinetic energy from the total energyU=Etot−T,
(11) which includes contributions from the mean field, pairing energy, and Coulomb energy. Based on the predicted kinetic and potential energies for a nucleus, the ratio
T/U for the nucleus in its ground state can be obtained.The axially deformed solution of the Skyrme HFB equations can be obtained by using the transformed harmonic oscillator basis [16−18, 21−23], and the corresponding code HFBTHO (v2.00d) was used in the calculations. A series of Universal Nuclear Energy Density Functionals (UNEDFs) [24] are proposed to provide a more accurate description of the properties of the ground state nuclei. The HFBTHO (v2.00d) also implements HFB calculations adopting the UNEDF functional. In this work, the UNEDF0 functional and standard Skyrme functionals SLy4 [25], SkM* [26] and SIII [27] were adopted. The number of oscillator shells was set as
Nshellsmax=20 and the total number of states in the basis wasNstates=1771 . The default valueb0=−1.0 fm was used for the oscillator length, which means that the code automatically setsb0=√ℏ/mω withℏω=1.2×41/A1/3 [17]. The pairing force was assumed to beVn,ppair=Vn,p0(1−12ρ(→r)ρ0)δ(→r−→r′) withρ0=0.16fm−3 , and a pre-defined pairing strengthVn,p0 was used for each Skyrme force for simplicity. The default value adopted for the quasi-particle cutoff energy wasEcut=60.0 MeV. -
Based on the Skyrme HFB code HFBTHO (v2.00d) [17], we systematically investigated the ground state properties of even-even nuclei with
8≤Z≤130 . First, we checked the calculated binding energies by adopting different Skyrme parameter sets. In Fig. 1, we show the calculated ground state energiesEtot for the isobaric chains withA=40 , 100, 150, 208, 228, and 256 using UNEDF0, SLy4, SkM* and SIII. The experimental data from AME2020 [28] are also presented for comparison. The data are reasonably well reproduced for medium mass nuclei with all four Skyrme forces. For heavy nuclei withA>200 , the uncertainties of the predictions from these Skyrme forces significantly increase, especially for neutron rich nuclei. Note that the calculated results with UNEDF0 are the best for nuclei shown in Fig. 1. UNEDF0 demonstrates better predictive capabilities in the region of heavy nuclei, possibly due to its refined treatment of isospin dependence. We therefore adopted UNEDF0 for the following analysis of the total energies and kinetic to potential energyT/U ratios.Figure 1. (color online) Nuclear ground state energies
Etot for nuclei withA=40 , 100, 150, 208, 228, 256. The curves denote the calculated results with UNEDF0, SLy4, SIII, and SkM*. The red circles denote the experimental data [28].The nucleus with the minimal energy for certain light and medium isobars can be considered as the most stable nuclide on this isobaric chain as mentioned previously. According to the available data [28], in Fig. 1, the corresponding known stable nuclei (e.g.
40 Ar,40 Ca,100 Ru,100 Mo,150 Nd,150 Sm and208 Pb) are located around the valleys of the ground state energies for the isobaric chains withZ≤82 . In addition to the nuclear ground state energyEtot , we concurrently investigated the corresponding ratioT/U . Fig. 2 shows the ground state energyEtot andT/U ratio for isobaric chains withA=40 , 48, 208, 298, 90, and 120. The nuclei with minimal energies are generally those with maximal values ofT/U , which indicates that the ratioT/U also has a close relationship with nuclear stability in addition to nuclear binding energy. The values ofT/U essentially show the relationship between the average kinetic and potential energies per nucleon, reflecting the relevance of nuclear force strength to the intensity of nucleonic motion. A larger average kinetic energy per nucleonT/A may indicate a more vigorous nucleon motion. Similarly, a larger absolute value of the average potential energy per nucleon|U|/A indicates a stronger attractive potential acting on the nucleon. Nuclear stability is therefore significantly influenced by the balance of T and U.Figure 2. (color online) Comparison of the ground state energy
Etot and the ratioT/U for isobaric chains withA=40 , 48, 208, 298, 90, and 120. The dashed lines serve as a guide.In Fig. 2, notably the isobars with
A=40 and 48 and maximal values ofT/U are located at the double magic nuclei40 Ca and48 Ca, respectively, whereas the corresponding values of minimalEtot are found at nuclides40 Ar and48 Ti. To investigate the influence of shell closure onEtot andT/U , we further analyzed theT/U andEtot for nuclei around208 Pb. The predicted ground state energyEtot and correspondingT/U ratio for isobaric chains withA=202 , 204, 206, and 210 are shown in Fig. 3. The peaks ofT/U for these isobars are all located at proton magic numberZ=82 . For210 Pb with proton shell closure and210 Po with neutron shell closure, the predicted values ofT/U are very close to each other. Apparently, the known magic numbers can be more clearly observed from the ratioT/U , which can be helpful for exploring the magic numbers in the super-heavy region.Figure 3. (color online) The same as Fig. 2 but for isobaric chains with
A=202 , 204, 206, and 210.To further analyze the difference between
Etot andT/U , we systematically investigated theT/U ratios for even-even nuclei with8≤Z≤130 . Fig. 4 shows the nuclei with maximalT/U on a given isobaric chain predicted by the HFBTHO calculations with UNEDF0. The solid curve denotes the predicted β-stability line according to Green's formula [29],Z=A/2[1−0.4A/(A+200)] . In the regions with known stable nuclei, the predicted nuclei with maximalT/U are generally located around the β-stability line. In addition, the calculated nuclei with minimal total energyEtot are also presented for comparison. In the super-heavy mass region, the number of nuclei obtained fromT/U is much more than that fromEtot atN=184 . Note that the nuclei withZ=108 andZ=118 obtained fromT/U are far more numerous than the neighboring nuclei.Figure 4. (color online) Nuclear β-stability line. The gray squares denote the known even-even stable nuclei. The circles and red squares denote the calculated even-even nuclei with maximal
T/U and those with minimalEtot , respectively. The solid curve denotes the results of Green's formula.Notably, the central position of the "island of stability" is closely related to the microscopic shell correction energies of nuclei which are usually estimated by subtracting the smooth macroscopic part from
Etot . In this work, we studied the corresponding microscopic energies of nuclei by combining the modified Bethe-Weizsäcker mass formula [30, 31]ELD(A,Z)=avA+asA2/3+acZ(Z−1)A1/3(1−Z−2/3)+asymI2A.
(12) Here, the parameters in Eq. (12) are determined by fitting the calculated ground-state energies of 2138 even-even nuclei with UNEDF0. In Fig. 5, we show the contour plots of the calculated microscopic energies
Etot−ELD for nuclei withZ≥90 . In the regionN>126 , one can see two islands: one located aroundN=152−162 ,Z=102−108 and the other located aroundN=184 , which is consistent with the predictedT/U ratios. The predicted microscopic energy for the super-heavy nucleus (Z=120,N=184 ) is approximately 7 MeV, which is comparable with the results of macroscopic-microscopic mass models [32, 33].Figure 5. (color online) Microscopic energy (in MeV)
Etot−ELD based on the ground state energy from UNEDF0 and smoothed liquid-drop energyELD from the modified Bethe-Weizsäcker mass formula. The dashed lines serve as a guide.In Fig. 6(a), the peak values of
T/U are displayed for each isobaric chain as a function of mass number A. For heavy nuclei, the maximum ofT/U in general decreases linearly with A. The oscillations in the HFBTHO+UNEDF0 results are mainly due to the microscopic shell effects. For light stable nuclei withA<50 , the contribution of nuclear surface energy is evident from Eq. (12), and the average kinetic energy per nucleon is relatively low, which results in a largerT/U ratio compared with those of the heavy nuclei. In Fig. 6 (b), the average kinetic energy per nucleon significantly increases with mass number A in the regionA<150 . In the heavy and superheavy mass region (A>150 ), the average kinetic energy per nucleon generally approaches a constant due to nuclear density saturation. Coulomb energy is known to significantly increase with charge number Z in the heavy mass region, which results in a reduction of the total potential energy in absolute value. For stable nuclei in the medium-heavy mass region, the nearly linear dependence ofT/U on mass number A is clearly observed in Fig. 6(a), which is mainly due to the competition between the significant increase in Coulomb energy and saturation behavior of the nuclear force.Figure 6. (color online) (a) Maximum of the ratio
T/U for each isobaric chain as a function of mass number A. (b) The corresponding average kinetic energy per nucleonT/A (dashed curve) and average potential energy per nucleonU/A (solid curve).The virial theorem states that for an inverse square force field (such as Coulomb potential), the ratio of kinetic to potential energy is
T/|U|=0.5 . According to the HFBTHO calculations, the absolute values ofT/U for bound nuclei are approximately0.65∼0.75 , which implies that the short-range characteristic and complexity of the nuclear force may lead to a different virial theorem. The linear dependence ofT/U on mass number A for heavy nuclei can be helpful in further exploring the virial theorem of nuclear forces and equation of state for neutron stars. -
In this study, we investigated the nuclear ground state energies and the ratios of kinetic energy T to potential energy U using HFBTHO. Based on the Skyrme energy density functional UNEDF0, the ratios
T/U for 2318 even-even nuclei were systematically calculated. We note that the nuclei with maximal value ofT/U for a certain isobaric chain are generally stable or long-lived. Also note that the known magic numbers can be more clearly observed from the ratioT/U than from nuclear binding energy, particularly for isobaric chains with semi-magic nuclei (with magic numbers for either proton or neutron). Combining nuclear binding energies with the values ofT/U from the Skyrme HFB calculations, the magic numbers in super-heavy mass region were studied. In the super-heavy mass region, the neutron magic numberN=184 was clearly observed from the obtainedT/U ratios and extracted microscopic energy of nuclei using UNEDF0. We also found that the nuclei withZ=108 andZ=118 obtained fromT/U are much more numerous than the neighboring nuclei, which is helpful for exploring the shell structure of superheavy nuclei. For super-heavy nuclei withZ=118 andN=178−196 , the calculated deformations of the nuclei were spherical or oblate in shape, whereas the predicted292−300 Fl were all spherical in shape. Therefore, further investigations into the influence of nuclear oblate deformations on the stability of super-heavy nuclei are required in the future. TheT/U peak value can be used as a supplementary indicator in selecting potential long-lived superheavy nuclei. In addition, we note that the obtainedT/U ratios decreased almost linearly with mass number A for heavy nuclei around the β-stability line due to the competition between strong Coulomb repulsion and the saturation property of the nuclear force, which implies that the short-range characteristic and complexity of the nuclear force may lead to a different virial theorem.
