Microscopic Calculations of the First 2⁺ State: Influence of Skyrme Functional Details and nucleon Pairing Interaction

The understanding of how protons and neutrons arrange themselves into diverse collective and single-particle configurations, including low-lying excited states and rotational bands, is one of the central objectives of nuclear physics. The study of collective excitations in superheavy nuclei (SHN) has become a central theoretical pursuit, as these excitations directly probe fundamental properties like stability, deformation, shell structure, and pairing dynamics. Among these, estimating the first $ 2^+ $ state energy, $ E(2^+) $, is paramount, as it quantifies collective rigidity and deformation, key determinants of stability. Accurate theoretical prediction of $ E(2^+) $ prior to synthesis is essential, as it establishes the framework for interpreting subsequent measurements. A match between predicted and observed values provides decisive insight into whether a newly formed nucleus lies within a potential island of stability or on the brink of rapid fission, thereby serving as a critical benchmark for validating nuclear models. In well-deformed nuclei, $ E(2^+) $ is inversely proportional to the moment of inertia, characterizing their collective rotational behavior. This trend aligns with Grodzins’ empirical correlation between decreasing $ E(2^+) $ values and increasing quadrupole deformation [1]. Near shell closures, $ E(2^+) $ increases significantly due to enhanced quadrupole stiffness [2]. Microscopic extensions of Grodzins-type scaling predict $ E(2^+) $ systematics for $ Z \geq 100 $, identifying characteristic minima in mid-shell regions and sharp maxima near $ N\approx184 $ and $ Z\approx118 $ [3]. These theoretical trends are further supported by recent self-consistent HFB calculations of the $ Z = 120 $ isotopes, which reveal a minimum in $ E(2^{+}) $ near $ N = 170 $ followed by a systematic increase as the nuclei approach the spherical $ N = 184 $ shell closure [4]. These signatures are increasingly informed by and correlated with α-decay fine-structure analysis. In the absence of prompt γ-spectroscopy, transitions to the lowest $ 2^{+} $ states serve as critical structural probes [4, 5].

Predictions of $ E(2^+) $ from microscopic and macroscopic-microscopic (MM) calculations are marked by significant theoretical uncertainty, reflected in the wide spread of results across different models and, within a single model, across various Skyrme parameterizations and other forms of the effective nucleon-nucleon (NN) interaction. This wide variation stems from differing approximations of the nuclear interaction and the inclusion of many-body correlations, reflecting persistent uncertainties in the foundational modeling of collective behavior. Quantifying these uncertainties systematically across the nuclear chart is therefore critical. Only by doing so can predictions be meaningfully interpreted, distinguishing true signatures of stability in SHN and guiding targeted experimental searches toward regions where theoretical input is most needed.

Early systematic approaches to binding energy and collective behavior were anchored in the liquid-drop model (LDM), which reproduced broad trends in masses, fission, and macroscopic dynamics through its surface, Coulomb, and asymmetry terms [6]. However, while successful for gross properties and fission, the LDM remained incapable of explaining magic numbers and the binding energy irregularities that reveal shell effects. The shell model, developed in the late 1940s [7, 8], revealed that the nuclear mean field and spin-orbit coupling generate energy gaps that stabilize specific configurations. This inherent duality, between collective liquid-drop behavior and quantized single-particle structure, has since become a defining framework for modern nuclear models. The phenomenological unification of this interplay is achieved through the MM approaches [8], which augment the macroscopic liquid-drop energy with shell and pairing corrections derived from single-particle phenomenological potentials, such as Nilsson or Woods–Saxon [9]. While the macroscopic component has been enhanced for fission studies through specialized models, other approaches now derive the macroscopic energy directly from the Skyrme energy density functional (EDF), creating a stronger link to microscopic nuclear interactions [10]. Refined hybrid MM models, featuring shell corrections from self-consistent spectra, offer a unique blend of computational speed, predictive accuracy, and microscopic alignment, rendering them essential tools for systematic SHN studies [11−17].

The limitations of the MM approach motivated the transition toward fully self-consistent mean-field frameworks. This transition has been enabled by the substantial increase in computational power and by the development of efficient numerical solvers, such as HFBTHO [18], which have made large-scale, fully self-consistent global surveys feasible and allow the inclusion of several nuclear properties beyond the scope of MM models. The Hartree-Fock (HF) method laid the foundation for these approaches by providing a self-consistent single-particle potential derived from effective forms of the NN interaction. The Hartree-Fock-Bogoliubov (HFB) formalism, a generalization of HF theory, introduces pairing correlations through the Bogoliubov transformation of creation and annihilation operators. This provides a unified mean-field treatment of both particle-hole and particle-particle pairing channels [19]. This framework eliminates the need for separate pairing corrections and provides a consistent description of spherical and deformed closed-shell nuclei and open-shell systems with pairing correlations.

At the core of both the MM and self-consistent HFB frameworks lies the Skyrme EDF [18], whose parameterization has been subject to frequent refinement over decades to reproduce an ever-broader set of observables. The development of Skyrme parameter sets progressed through distinct phases. Early global Skyrme parameter sets, such as SI–SVI, demonstrated the feasibility of a self-consistent HFB description of nuclear binding and structure [20]. A subsequent generation, including SkM, SkM* and SkM1, refined surface properties and fission barriers [21], while SkP pioneered self-consistent treatment of pairing correlations [22]. Later developments, such as those in the SkT family, including SkT5 [23], built upon this foundation to provide refined descriptions of nuclear matter (NM) properties and ground-state masses using extended parameter sets [23]. Following these initial efforts, developments produced specialized fits such as SkI5 [24] for isovector properties, RATP for fission barriers [25], and globally optimized functionals like SV-min [26]. The Skyrme-Lyon (SLy) series was developed to provide a robust description of isovector properties and neutron-rich matter [27]. Functionals from this series are consequently employed in theoretical studies of SHN and nuclear decays and reactions [28−32], where the behavior of neutron skins and asymmetric matter is important. In parallel, the treatment of pairing correlations in mean-field theories advanced through regularization techniques for zero-range forces, a methodological improvement relevant for functionals like KDE0v [33] and KDEX [34].

The Skyrme EDF of the effective NN interaction describes the nuclear force by augmenting the kinetic energy and the Coulomb force with six specific interaction contributions. The central part of the force is defined by four specific terms, namely a zero-range contact interaction, a density-dependent term, an effective-mass term, and a finite-range interaction. The interaction is further refined by including non-central contributions from spin-orbit and tensor couplings. The parameters for the Skyrme interaction are determined empirically by fitting to experimental nuclear structure data and theoretical NM properties. This fitting yields unique parameter sets, each of which predicts specific, definitive values for key NM properties such as the symmetry energy, its density slope, incompressibility, and the effective mass. The Skyrme EDFs have successfully provided a semi-microscopic, self-consistent framework for describing a wide range of nuclear phenomena. This includes applications to nuclear structure [35−37], reaction processes [38−40], and decays [41−43], as well as investigations into the properties of NM [44, 45] and astrophysical systems [46]. Different parameterizations of the Skyrme EDF may yield widely divergent results for certain collective properties, even when they reproduce ground-state and NM data equally well. This sensitivity arises when the parameterization is poorly constrained for specific phenomena, highlighting that a universal EDF form for all nuclear observables remains elusive. The hundreds of proposed Skyrme EDF parameter sets enable a systematic investigation of how variations in their terms impact the results of mean-field calculations. In this work, we investigate how the parameters of the Skyrme effective interaction and its underlying NM properties influence $ E(2^{+}) $ systematics in even-even nuclei, employing ten representative EDF parameterizations in HFB calculations. Our analysis also explores how the results depend on the specific treatment of neutron and proton pairing.

This paper is structured as follows. Section II outlines the theoretical framework, giving a brief summary of the HFB scheme as implemented in the HFBTHO code, and specifies the use of the Skyrme EDF as the effective NN interaction and a zero-range density-dependent pairing functional. In Section III, we present, analyze, and discuss the results, focusing on how the force parameters, their associated NM properties, and the parameterization of the pairing functional influence the calculated $ E(2^+) $. The summary and principal conclusions are provided in Sec. IV.

The HFB method provides a microscopic foundation for the description of nuclei. It approximates the nuclear many-body ground state as a quasiparticle vacuum generated by a Bogoliubov transformation of single-particle creation and annihilation operators. By coupling particle and hole operators, this transformation defines a unified ground state characterized by a generalized density matrix, which simultaneously describes particle occupancy via the normal density ρ and pairing correlations via the pairing tensor κ. The HFB equations are derived from a constrained variational minimization of the energy density functional with respect to the generalized density matrix, under particle-number constraints and, optionally, shape constraints, ensuring self-consistency between ρ and κ. In coordinate space and for a local EDF, the HFB equations can be written in terms of the quasiparticle wave functions $ U_k({\bf{r}},\sigma) $ and $ V_k({\bf{r}},\sigma) $, quasiparticle energies $ E_k $, the single-particle HF Hamiltonian h, the pairing field Δ, and the chemical potential λ, which enforces the average particle number, as

The HFBTHO solver [18, 47, 49] enforces axial symmetry, $ [{\hat{J}}_z,\ {\hat{H}}_{HFB}]=0 $, and time-reversal invariance, $ \hat{T}{\hat{H}}_{HFB}{\hat{T}}^{-1}={\hat{H}}_{HFB} $. This implies that the quasiparticle wavefunctions are simultaneous eigenstates of $ \hat{J}_z $ and the parity operator $ {\hat{\Pi}}_z $, characterized by quantum numbers $ \Omega^\pi $. The HFB matrix then decouples into independent blocks for each $ \Omega^\pi $ value, which are diagonalized separately in a truncated axially deformed harmonic-oscillator basis $ \{\phi_{n_z,n_\bot,\Lambda,\Sigma}\} $. Optionally, a local scaling transformation is applied to construct a transformed harmonic-oscillator (THO) basis, thereby improving the asymptotic behavior of single-particle wavefunctions at large radii.

The zero-range pairing interaction used in the HFBTHO scheme has the standard density-dependent form:

Here, the neutron (proton) pairing strength is denoted by $ V_{0n(0p)} $, ρ is the local isoscalar density, $ \rho_0 $ is the nuclear saturation density, and the density dependence is scaled by the exponent β. The coupling parameter $ V_1 $ controls the spatial character of the pairing interaction, ranging from volume ($ V_1 = 0 $) to surface ($ V_1 = 1 $) contributions, with intermediate values corresponding to a mixed coupling. To avoid pairing-energy divergences associated with zero-range pairing forces, a pairing regularization is implemented: it replaces the constant input strength with a position-dependent effective strength computed locally throughout the nucleus from the density, effective mass, and the energy cutoff ($ E_{cut} $), thereby subtracting the divergent contribution [47]. The present calculations employ a linear density dependence ($ \beta=1 $) with a mixed volume–surface coupling ($ V_1 = 0.5 $) and a quasiparticle energy cutoff of $ E_{cut}=60 $ MeV. The Lipkin–Nogami (LN) prescription is employed in HFBTHO calculations to reduce the spurious particle-number fluctuations inherent in the HFB approach, thereby improving the description of pairing correlations in nuclei with closed or near-closed shells [18].

To describe collective quadrupole dynamics beyond the static HFB minimum, HFBTHO implements the generator coordinate method (GCM) within the Gaussian overlap approximation (GOA). This involves superposing constrained HFB states labeled by a collective coordinate q, typically the axial quadrupole moment $ Q_{20} $. The description of large-amplitude collective motion within the GCM–GOA framework relies on calculating key quantities, including the collective mass tensor—derivable from either the GCM or Adiabatic Time-Dependent HFB (ATDHFB) approximations—and the zero-point energy correction [47]. A critical step in this process is the self-consistent readjustment of the corresponding Lagrange parameters, accomplished efficiently using the cranking approximation to the quasiparticle random-phase approximation (QRPA) matrix. This yields closed-form expressions for the GCM metric.

and the collective mass tensor (inverse inertia, B),

as well as the zero-point energy (ZPE)

These terms are determined by the underlying HFB solution and the associated local QRPA matrix, which is represented by the moment matrices.

where the summation extends over all two-quasiparticle states $ |\mu\nu\rangle $. For a given deformation characterized by the collective variables q, the collective potential incorporates quantum fluctuations and is given by

This computationally tractable, self-consistent framework is implemented in the HFBTHO code, where it enables the many-body problem to be recast into a collective Schrödinger-like equation driven by these derived parameters.

In the collective rotational approximation for well-deformed nuclei exhibiting axial symmetry, nuclear rotational motion is described by assuming a fixed, axially deformed intrinsic density distribution that rotates rigidly about an axis perpendicular to its symmetry axis. Consequently, the excitation energy of the first $ 2^{+} $ state can be approximated in terms of the perpendicular moment of inertia, $ \mathfrak{I}_\bot $, obtained as a specific component of the collective inertia tensor, via the $ K=0 $ rotor formula $ E(2^{+}) = 3\hbar^2/\mathfrak{I}_\bot $. Within the GCM-GOA framework, this moment of inertia can be derived microscopically from its dependence on the collective mass parameter ($ M_{22} $) associated with the axial quadrupole deformation coordinate and on the intrinsic quadrupole moment $ Q_{20} $. $ \mathfrak{I}_\bot/\hbar^2=M_{22}\times Q_{20}^2 $. Within the cranking approximation to the QRPA matrix, the GCM metric and collective inertia reduce to the Inglis-Belyaev (IB) expression [48]. Hence, the cranking-based adjustment of Lagrange parameters in the GCM-GOA framework parallels the rotational formula using the IB moment of inertia, to which the rotational energy is inversely proportional, computed from two-quasiparticle sums over the HFB vacuum. Both approaches rely on the ATDHFB cranking framework. The QRPA also mirrors the IB structure by supplying mass-tensor elements via Thouless-Valatin restoration. However, this approximation holds only for well-deformed nuclei with established rotational motion. In transitional regions, where vibrational modes dominate the low-energy spectrum, IB-like expressions overestimate excitation energies because the rigid cranking assumption breaks down. For such transitional nuclei, a full QRPA calculation provides more reliable results for collective excited states.

The Skyrme EDF defines the nuclear energy density ($ {\cal{H}}_{Sk} $) as a functional of local neutron and proton densities ($ \rho_{n,p} $), kinetic-energy densities ($ \tau_{n,p} $), and spin densities ($ {\bf{J}}_{n,p} $). It is constructed as a sum of specific terms representing different physical contributions to the total energy of the system.

Here, the total kinetic energy density ($ {\cal{H}}_K $) is given by the sum of the proton ($ \tau_p $) and neutron ($ \tau_n $) contributions.

The central part of the interaction, which governs bulk nuclear properties, is modeled by four distinct contributions: the zero-range $ {\cal{H}}_0 $, density-dependent $ {\cal{H}}_3 $, effective-mass $ {\cal{H}}_{eff} $, and finite-range $ {\cal{H}}_{fin} $ terms. $ {\cal{H}}_{so}\left(W_0,{\bf{J}}_i,{\bf{\nabla}} \rho_i\right) $ and $ {\cal{H}}_{sg} $ respectively represent the spin-orbit term and the tensor coupling, including the spin and gradient contributions of the EDF. The different parts of the Skyrme effective interaction are expressed in terms of the parameters $ t_{i=0,1,2,3} $ and $ x_i $, which correspond to the strengths of the standard and exchange terms, respectively, and the exponent σ that governs the density dependence. The strength of the zero-range spin-orbit force is given by $ W_0 $. The total spin-orbit density, J, is the sum of the proton, $ J_p $, and neutron, $ J_n $, contributions. These Skyrme parameters are typically fitted to reproduce properties of finite nuclei and infinite NM, as well as various nuclear phenomena. The Coulomb energy is the sum of its direct (dir) and exchange (exch) contributions, both calculated using the proton density distribution.

The exchange part is evaluated within the local Slater approximation [50]. The energy per nucleon in infinite asymmetric nuclear matter (ANM) is, based on the Skyrme EDF, given by $ E_A={\cal{H}}_{Sk}(\rho)/\rho $ [51, 52]. Expanding ϵ in terms of ρ and the isospin asymmetry (I) allows one to derive the symmetry energy, which quantifies the isospin dependence of the NN interaction. Key quantities that characterize the equation of state (EOS) of ANM and govern its behavior are the symmetry energy coefficient (J), the density slope (L) of the symmetry energy, and the incompressibility coefficient $ K_0=9\rho\left.\partial^2E_A\left(\rho\right)/\partial\rho^2\right|_{\rho_0} $, all defined at the saturation density ($ \rho_0 $). Additionally, the effective nucleon mass ($ m^* $) characterizes the modified inertia of a nucleon in NM. It is not a free parameter but a consequence of interactions, encoded in the dispersion relation. It is defined relative to the free nucleon mass (m) as $ m^*/m = (m/p) d\varepsilon_N/dp $, which ensures the correct nucleon velocity $ v=d\varepsilon_N/dp=p/m* $. Within the Skyrme EDF, the isoscalar (IS) and isovector (IV) effective masses, which reflect the symmetric and asymmetric properties of NM, respectively, are given as

and

Each Skyrme EDF is defined by a specific set of parameters that determine its predicted NM properties. These characteristics, in turn, shape the predicted isospin asymmetry and density dependence of the structure of finite nuclei and their reaction dynamics when the EDF is applied [30, 45, 52−55].

As the first stage of our analysis, we employ several variants of the Skyrme effective interaction to examine the isotopic and isotonic trends in binding energy and the excitation energy of the first $ 2^{+} $ rotational state in the even-even nuclei along the Th and $ N=134 $ chains. We perform microscopic, self-consistent mean-field calculations using the HFB method, implemented via the HFBTHO code [18, 49, 56], in which pairing correlations are treated with a density-dependent zero-range pairing interaction that is regularized to eliminate ultraviolet divergences and to yield a well-defined pairing field [18, 57]. To compute rotational excitation energies, the Generator Coordinate Method (GCM) combined with the Local Gaussian Overlap Approximation (GCM+GOA) is applied to the collective inertia tensor [47], as outlined in Sec. II. To this end, we examine a diverse set of Skyrme EDFs to cover wide ranges of the Skyrme parameters and the associated saturation properties of asymmetric NM. Our analysis specifically employs ten Skyrme parameterizations: SkI5 ($ m^*/m=0.58 $) [24], SGOI (0.61) [58], RATP (0.67) [25], SLy4 (0.69) [27], KDE0v (0.72) [33], SkM1 (0.79) [59], KDEX (0.81) [34], Es (0.84) [60], SV-min (0.95) [26], and SkT5 (1.00) [23].

In Fig. 1, panels (a) and (b) show, respectively, the binding energies and the energies of the $ 2^{+} $ rotational states for the Th isotopes, while panels (c) and (d) display the corresponding quantities for the $ N=134 $ isotones. The HFBTHO framework utilizes a volume-surface, zero-range pairing force with linear density dependence and distinct fitted strengths for neutrons and protons. In Skyrme EDFs, pairing is not an intrinsic part of the functional but is often added as a simpler, standard supplement—such as volume ($ V_1=0 $), surface ($ V_1=1 $), or mixed terms—without being consistently refitted. When considering a given EDF parameterization originally fitted without an explicit pairing component, or that uses a pairing form different from the one implemented in HFBTHO or similar solvers, the strengths of the added pairing functional must be adjusted separately to reproduce relevant empirical observables such as odd–even mass staggering and rotational properties, to ensure realistic pairing correlations in the final calculation. In the calculations presented in Fig. 1, the neutron and proton coupling strengths in the HFBTHO code are set to the average values $ V_{0n}=-300 $ MeV fm³ and $ V_{0p}=-310 $ MeV fm³. This aims to investigate how the strength of this pairing force compares to the effective parameters of the EDF that determine a nucleus's rotational behaviour.

Figure 1.  (color online) Panels (a) and (c) show the calculated ground-state binding energies, and panels (b) and (d) show the energies of the first $ 2^{+} $ rotational state, for Th isotopes (a, b) and $ N=134 $ isotones (c, d), computed with ten Skyrme EDFs and compared with the corresponding experimental data [63]. Calculations were performed within the HFBTHO framework, following the methodology detailed in Secs. II and III, using average pairing strengths of $ V_{0n}=-300 $ MeV fm³ and $ V_{0p}=-310 $ MeV fm³.

Fig. 1 shows that most of the considered parameterizations reproduce the binding energies within acceptable agreement with the data, with SGOI exhibiting the largest deviation. This can be traced to its notably large symmetry-energy ($ J = 45.2 $ MeV) and incompressibility ($ K_0 = 361.59 $ MeV) coefficients relative to the other forces. In addition, it has the second-weakest zero-range spin-orbit strength in the set, $ W_0=102.6 $ MeV fm⁵, second only to KDEX, and it simultaneously features the least negative ratio, $ t_2/t_1=-0.15 $, i.e., the smallest magnitude among the forces with $ t_2/t_1<0 $, as well as the largest surface-energy coefficient, 18.5 MeV. The $ t_2/t_1 $ ratio is a key indicator of the nonlocality and momentum dependence of the EDF, with important consequences for single-particle and collective nuclear properties.

For the energy of the first $ 2^{+} $ state, the SLy4, RATP, KDE0v, and SGOI parameterizations, using the provided proton and neutron pairing coupling strengths, provide the best agreement with the data. Among the considered forces, these four parameterizations have an exchange parameter of the zero-range term ($ x_0 $) between 0.42 and 0.63 and yield effective masses in the range $ 0.61\,m $ to $ 0.72\,m $. The remaining force parameters do not systematically distinguish these four from the others. For SkI5, a reduced effective mass ($ 0.58\,m $) and a large symmetry-energy slope ($ L=129.33 $ MeV), together with a large, unbalanced $ t_0/t_3 $ ratio of $ -0.22 $ fm^−3/4, lead to an underestimated $ E(2^{+}) $. The SLy4, RATP, and KDE0v forces share symmetry-energy slopes in the range $ L \approx 32 $–45 MeV and have $ t_0/t_3 $ ratios around −0.18 fm^−3σ. Increasing the effective mass in the other parameterizations, or shifting $ x_0 $ outside its characteristic range for the four highlighted interactions, produces an overestimated $ E(2^{+}) $ by disrupting the quadrupole deformation. The poorest results are associated with the lowest spin-orbit strengths, specifically $ W_0 $ values of 99–102 MeV fm⁵ for the KDEX and SV-min parameterizations. For the Th isotopes shown in Fig. 1, the four interactions SLy4, RATP, KDE0v, and SGOI produce pairing gaps under the Lipkin-Nogami scheme in the ranges 0.55–1.26 MeV, 0.55–1.24 MeV, 0.56–1.21 MeV, and 0.47–1.26 MeV, respectively. The corresponding ranges for the N=134 isotones are 0.55–0.78 MeV, 0.57–0.82 MeV, 0.59–0.83 MeV, and 0.49–0.70 MeV, respectively. The results based on SLy4, RATP, and KDE0v are most consistent with the three-point pairing gaps derived from experimental nuclear masses [63], which yield a range of 0.73–1.27 MeV (0.85–1.04 MeV). The three-point pairing gaps for these three interactions yield smaller pairing values than those obtained directly from the HFB+LN calculations. For instance, the SLy4 (RATP) interaction yields pairing gaps of 1.26 MeV (1.24 MeV) and 0.81 MeV (0.78 MeV) for the ²¹⁶Th isotope based on the direct HFB+LN calculations and the three-point difference calculations, respectively, while its empirical three-point pairing gap is 1.27 MeV [63].

To improve the systematic description of $ E(2^{+}) $ energies, we selected the SLy4, KDE0v and SGOI interactions based on their performance in Fig. 1, with $ V_{0n}=-300 $ MeV fm³ and $ V_{0p}=-310 $ MeV fm³. This framework was then applied to a broader set of nuclei. Figure 2(a) shows results for even-even ^{128−140,144−160}Nd, ^236−250Cm, and ^242−254Cf isotopes alongside ^220−236Th isotopes, while Fig. 2(c) displays $ N=60 $ ($ A=98-114 $), 100 (160, 164−182), 144 (232−242), 148 (240−248), and 150 (242 − 254) isotones alongside $ N=134 $ ($ A=218-226 $). Additional calculations for ^248−256Fm, ^252,254No, and ^254,256Rf isotopes and for $ N=152 $ ($ A=246-256 $), 154 (250 − 254), and 156 (254, 256) isotones were performed but omitted from the figure as they are entirely coincident with the displayed data, which would obscure the figure. For these sets of isotopes and isotones, the standard deviation (SD) between calculated and observed $ E(2^{+}) $ values decreases from $ \sigma=0.197 $ MeV (SLy4) and 0.181 MeV (KDE0v) to 0.170 MeV (SGOI). Standard deviations $ \sigma=\sqrt{\sum_{i=1}^{n}\left(x_i-\bar{x}\right)^2/(n-1)} $ are calculated from the residuals between calculated and observed values where $ x_i = x_{\rm{cal}} - x_{\rm{obs}} $, $ \bar{x}=\sum x_i/n $ is the mean residual, and n is the number of data points. However, the SGOI EDF performs less favourably in terms of both mean and median absolute errors and in predicting low $ E(2^{+}) $ values. It yields a mix of over- and under-predictions, a pattern clearly visible in Fig. 2(a) where Cm isotopes are systematically underestimated and Cf isotopes overestimated. The smallest SD for SGOI indicates that its errors are consistent. Yet, it shows a pronounced positive bias at higher $ E(2^{+}) $ values, pointing to systematic overprediction in that region, while tending to underestimate lower $ E(2^{+}) $ values. Consequently, SGOI can be described as precise but not accurate. In contrast, KDE0v achieves greater overall accuracy on average, even though it is slightly less precise than SGOI. The SLy4 EDF yields intermediate mean and median errors relative to the SGOI and KDE0v functionals. For binding energies, the agreement is markedly better for SLy4 ($ \sigma = 3 $ MeV) than for KDE0v (31 MeV) or SGOI (28 MeV). SLy4 thus delivers the most accurate results for both $ E(2^{+}) $ excitation energies and binding energies among these EDFs. The quadrupole deformations ($ \beta_2 $) obtained with the SLy4 and KDE0v interactions, which maintain simultaneous best average accuracy for the calculated binding energy and $ E(2^{+}) $ for the isotopes in Fig. 2(a) and isotones in Fig. 2(c), are presented in Figs. 2(b) and 2(d), respectively. The deformations from the two interactions agree with each other. They reveal the same regions of prolate and oblate deformations and their trend behavior with both Z and N, typical maximum and minimum deformations, as well as consistent candidate regions for prolate/oblate shape coexistence, particularly within Z=44–46 and Z=78–80 in Fig. 2(d).

Figure 2.  (color online) HFB calculations of $ E(2^{+}) $ energies for even-even (a) Nd, Th, Cm, and Cf isotopes and (c) $ N=60,100,134,144,148, $ and 150 isotones, using the KDE0v, SGOI, and SLy4 parameterizations, with pairing strengths $ V_{0n} = -300 $ MeV fm³ and $ V_{0p} = -310 $ MeV fm³. (b) and (d) show the corresponding quadrupole deformations obtained with the SLy4 and KDE0v interactions for the isotopes in panel (a) and the isotones in panel (c), respectively.

Since the SLy4 EDFs, with their initially considered neutron and proton coupling strengths, exhibit a pronounced positive bias and systematically overpredict $ E(2^{+}) $ values, we performed an optimization of their coupling parameters to better align with experimental data. Optimization led to revised values of $ V_{0n}=-285 $ MeV fm³ and $ V_{0p}=-300 $ MeV fm³, which significantly improved the accuracy of the model for $ E(2^{+}) $ predictions. The calculations performed based on the two coupling sets for the SLy4 EDF are expanded to encompass additional even-even nuclei, including ^120,138,140Te and ^{172,174,180−202}Pt isotopes, as well as the $ N=110 $ ($ A=182,186,188,194 $) and $ N=132 $ (216, 218) isotones, supplementing the previously considered sets of isotopes and isotones. Figure 3(a) shows results for the considered Te, Nd, Pt, Po, Th, Cm and Cf isotopes, and Fig. 3(b) presents analogous results for isotonic sequences at $ N = 60,100,132,134,144,148, $ and 150. To maintain clarity in the two figures, the remaining overlapping data points that coincide entirely with the presented results have been omitted to prevent visual overcrowding. Using the optimized coupling strengths, the overall SD across all nuclei improved to $ \sigma = 0.199 $ MeV compared to $ \sigma = 0.227 $ MeV for the earlier coupling strengths. The SD for the earlier subset of nuclei is also reduced by employing the optimized coupling strengths, achieving $ \sigma = 0.161 $ MeV. A modest reduction in binding energy uncertainty was also observed across the full set of nuclei, with σ decreasing from 3.187 MeV to 3.177 MeV. The optimized pairing strengths exhibit improved agreement with experimental $ E(2^{+}) $ values across the entire dataset. This is reflected in smaller, more frequently near-zero residuals compared to the earlier parameterization. Specifically, the refined strengths produce lower errors (< 50%) in 43% of cases for moderate deviations and 13% for the larger deviations, while also reducing the mean absolute error from 59.3% to 44.0% and the root-mean-square error from 86.8% to 70.7%, indicating a significant reduction in average prediction error. Both parameterizations exhibit a moderate positive correlation with the experimental values, with the refined set showing a marginally weaker correlation. However, the original strengths systematically overpredict the measurements, often with larger deviations. Below the shell closures at $ Z=82 $ or $ N=126 $, the two parameterizations alternate in performance. Above these closures, the reduced pairing strengths provide significantly improved agreement with experiment, reducing the substantial overpredictions seen with the preliminary strengths. This improvement is evident for the Th, Cm, and Cf isotopes in Fig. 3(a) and for the isotones of $ N=132,134,144 $ and 148 in Fig. 3(b).

Figure 3.  (color online) HFB calculations of the $ E(2^{+}) $ energies for even-even (a) Te, Nd, Pt, Po, Th, Cm, and Cf isotopes and (b) $ N = 60,100,132,134,144,148, $ and 150 isotones, computed using the SLy4 EDF with two sets of pairing strengths $ (V_{0n}, V_{0p}) = (-300 $ MeV fm³, −310 MeV fm³) and (−285 MeV fm³, −300 MeV fm³).

Increasing pairing strength triggers two distinct model responses. The decrease in calculated $ E(2^{+}) $ with reduced pairing coupling strengths, evident in Figs. 3(a) and 3(b), is consistent with the pairing-driven softening response. The QRPA matrix elements governing the $ 2^{+} $ excitation depend on the unperturbed 2qp energies and the residual interactions derived from pairing and particle-hole channels. As the pairing field weakens, the resultant reduction in 2qp energies shifts the balance in the QRPA eigenvalue problem toward more collective, lower-energy modes. Consequently, the QRPA eigenvalue for the lowest $ 2^{+} $ mode decreases, indicating a softer collective response despite a modest enhancement in the underlying single-particle level density and in the HFB shell gaps. An alternative deformation-driven influence produces $ E(2^{+}) $ increase with reduced pairing, where weaker pairing decrease the pairing gap, increasing the energy spacing between the single particle levels near the Fermi surface. As a result, shell gaps are enhanced, which suppresses the intruder orbitals responsible for driving deformation, resulting in a smaller intrinsic quadrupole moment ($ Q_2 $) and decreasing in the moment of inertia of the collective rotor. This yields higher $ E(2^{+}) $. In the vicinity of shell closures, the QRPA pairing-softening effect generally exceeds the influence of static deformation adjustments. This prevalence steadily diminishes as the number of neutrons (protons) increases toward mid-shell occupation. For example, lowering the pairing strengths to their determined optimum values for the Nd, Pt, Cm, and Cf isotopes results in a greater estimated reduction in $ E(2^{+}) $ energies as N approaches 82,126, and 152, and it yields a smaller reduction as N approaches 88,106, and 146. The corresponding reduction for the $ N=60,100,134,144, $ and 150 isotones increases as Z approaches the vicinity of 50, 82, and 104, while a smaller reduction is obtained as Z approaches 44, 74, and 92.

The large overestimated value of $ E(2^{+}) $ based on the optimized pairing strengths are obtained for the nuclei ^246,248Cf_148,15, ¹⁶⁶Dy₁₀₀, ^254,256Rf_150,152, ¹⁸²Hf₁₁₀, ¹⁸⁶Os₁₁₀ and ^190−196Pt_112−118, with deviations from experiment ranging from 0.05 MeV to 1.1 MeV. These nuclides occupy strongly deformed regions beyond $ N=100 $, where the theoretical framework overestimates quadrupole moments due to several inherent deficiencies. Key among these is the inadequate treatment of tensor and spin-orbit couplings, Skyrme functionals lacking full tensor terms lead to the underbinding of neutron intruder orbitals, such as $ \nu i_{{13/2}^+} $ and $ \nu j_{{15/2}^-} $ and the deformed states arising from them, which in turn artificially enhances the predicted deformation. This is further amplified by the strong correlation between the predicted deformation and occupations of the deformed single-particle shell orbits. A contributing factor to this discrepancy could be the "pairing collapse" endemic to the superheavy region, for instance in ^246,248Cf and ^254,256Rf. This is similar to what occurs in standard BCS calculations, which spuriously enhance the occupation of deformed orbits across the $ N=152 $ neutron shell. Additional shape degrees of freedom are also relevant. The present calculations indicate a non-negligible negative hexadecapole deformation ($ \beta_4 $) as a defining characteristic of ¹⁸²Hf and ¹⁸⁶Os in their ground-state, whereas the ^192,196Pt isotopes exhibit a less pronounced negative $ \beta_4 $ value. Such negative hexadecapole deformation alters the underlying single-particle level structure and the wavefunctions of the ground and excited states, thereby influencing the calculated $ E(2^{+}) $ value. Furthermore, nuclear shape complexities such as rapid shape coexistence and γ-softness, or octupole-driven correlations, are also probable in isotopes of Dy, Hf, Os, and Pt. The possibility of shape coexistence aligns with the indicated prolate-to-oblate ground-state shape transition for nuclei above $ N=108 $ [61, 62]. For instance, the present calculations indicate shape coexistence in ¹⁹⁰Pt between oblate ($ \beta_2=-0.159 $) and prolate ($ \beta_2=0.185 $) configurations, with an energy difference of less than 0.09 MeV. The estimated energy separation between the two configurations increases to approximately 0.46 MeV for the isotopes ^192,196Pt.

In addition to the pairing strength, other critical parameters of the zero-range pairing force given by Eq. (2) that influence the calculated $ E(2^{+}) $, and its overestimation in the mentioned 10 nuclei, are the volume-surface coupling parameter ($ V_1 $), and the quasiparticle energy cutoff ($ E_{cut} $) that drives the pairing energy regulation. For instance, further reduction in the pairing strength from ($ V_{0n}=-285 $ MeV fm³, $ V_{0p}=-300 $ MeV fm³) to (−250 MeV fm³, −260 MeV fm³), causes an additional lowering of the estimated $ E(2^{+}) $ for ²⁵⁶Rf from 0.179 MeV to 0.1 MeV, with the predicted deformation $ \beta_2 $ remaining nearly constant at 0.281. Increasing the energy cutoff $ E_{\rm{cut}} $ from 60 MeV to 80 MeV, with pairing strengths (−285 MeV fm³, −300 MeV fm³), results in a modest reduction of the estimated $ E(2^{+}) $ value to 0.176 MeV. Here, increasing the quasiparticle energy cutoff expands the single-particle configuration space included in the pairing calculation. This enables the pairing field to develop more fully, leading to a marginal increase in the average pairing gap and a consequent slight reduction in $ E(2^{+}) $ due to enhanced core softness.

On the other hand, increasing the volume-surface mixing parameter ($ V_1 $) from 0.5 to 0.75 reduces the estimated value of $ E(2^{+}) $ for ²⁵⁶Rf from 0.179 MeV to 0.073 MeV. Larger surface pairing favors soft surface modes over volume stiffness. Increasing the volume-surface mixing parameter enhances surface-peaked pairing correlations in the low-density surface region, which strengthen pairing in deformed configurations by softening the low-energy collective modes through greater density-gradient sensitivity and reduced central pairing, leading to significantly reduced $ E(2^+) $ values. While multi-parameter optimization remains an ideal long-term goal, the targeted refinement of $ V_{0n(p)} $ remains a physically motivated and universally justifiable strategy for improving the description of low-energy collective excitations. This single-parameter adjustment provides a sufficient reduction in the excitation energy that closely matches experimental data across the entire set of investigated nuclei preserving the global predictive power of the EDF for properties beyond $ E(2^+) $. It directly scales the pairing gap and enhances rotational softening, without the structural side effects introduced by altering the volume-surface parameter. This avoids the overcorrection produced by increasing the volume-surface mixing parameter, which not only risks mismatches in shell gaps and deformation but also contradicts the physical expectation near shell closures where reduced level density favors a decreased mixing parameter toward volume-dominant pairing over surface enhancement. For physical consistency in Skyrme-HFB models, multi-parameter optimization, particularly tuning the volume-surface mixing coefficient alongside pairing strength, must incorporate a broader set of nuclear structure observables beyond $ E(2^+) $, such as charge radii, surface diffuseness, odd-even mass staggering, pairing gaps, nucleon separation energies, and full low-energy collective rotational spectra, in addition to neutron skin effects across isotopic chains.

To examine isotopic and isotonic trends in the first $ 2^{+} $ state energy for even-even nuclei from Sr to Rf ($ N=60 $–156), we performed HFB calculations using the HFBTHO code. Pairing was treated with a regularized density-dependent zero-range interaction to prevent ultraviolet divergences, while the collective inertia tensor was obtained via the GCM+GOA method. Our study utilized ten different Skyrme EDFs, including SkI5, SGOI, RATP, SLy4, KDE0v, SkM1, KDEX, Es, SV-min, and SkT5, encompassing a wide spectrum of parameterizations and NM properties to test the robustness of the observed trends, considering $ V_{0n}=-300 $ MeV fm³ and $ V_{0p}=-310 $ MeV fm³. Most parameterizations reproduce binding energies reasonably well, though the SGOI EDF shows the largest deviation due to its high symmetry energy, high incompressibility, weak spin-orbit strength, large surface energy, and a $ t_2/t_1 $ ratio closest to zero among negative values. The magnitude of the reduction in the $ t_2/t_1 $ ratio to a negative value governs the non-locality and momentum dependence of the effective interaction, thereby influencing single-particle energies and collective excitations. Agreement with the experimental $ E(2^{+}) $ is best achieved, with the provided pairing coupling strengths, by the interactions sharing a constrained exchange parameter ($ x_0 \sim 0.42 $–0.63) and yield effective masses of 0.61–0.72 m. Underestimation is associated with a lower effective mass, a high symmetry-energy slope, and a larger unbalanced magnitude of negative $ t_0/t_3 $ ratio, while overestimation results from raising the effective mass or shifting $ x_0 $ beyond its indicated characteristic range. The discrepancy from the data increases as the spin-orbit interaction weakens, particularly for $ W_0 $ values below 120 MeV fm⁵.

Discrepancies in $ E(2^{+}) $ serve as indicators of pairing strength, overestimation implies stronger pairing is needed, while underestimation implies weaker pairing. Optimized pairing strengths of $ V_{0n}=-285 $ MeV fm³ and $ V_{0p}=-300 $ MeV fm³ significantly improve the agreement between calculated $ E(2^{+}) $ energies using the SLy4 parameterization and experimental values. The observed reduction in the calculated $ E(2^{+}) $ with diminished pairing strength is consistent with a pairing-driven softening of the collective response. This occurs because a weaker pairing field lowers 2qp energies, making the QRPA solution more collective and shifting the lowest $ 2^{+} $ mode to lower energy, which yields a softer response. A deformation-driven influence competes to raise $ E(2^{+}) $ as pairing weakens, where the reduction in pairing strength decreases the pairing gap and widens the energy spacing of single-particle levels near the Fermi surface. This amplifies shell gaps, which in turn suppresses the intruder orbitals that drive deformation. This yields a smaller intrinsic quadrupole moment and a lower moment of inertia for the collective rotor, producing a higher $ E(2^{+}) $ energy. The pairing-softening effect is strongest near shell closures, surpassing the role of static deformation. This dominance decreases progressively as the nucleon numbers approaches mid-shell. Also, the calculated $ E(2^{+}) $ decreases slightly when the quasiparticle energy cutoff is increased. This occurs because expanding the configuration space for pairing regularization allows the pairing field to develop more fully, yielding a marginally larger average pairing gap and thereby enhancing core softness. Furthermore, the predicted $ E(2^{+}) $ decreases substantially with an increased volume-surface mixing parameter of the pairing interaction. This is a direct consequence of enhanced surface-peaked pairing, which softens collective modes by strengthening pairing in deformed configurations. The increased density-gradient sensitivity and reduced central pairing associated with the surface component are the key mechanisms behind this softening and the resultant drop in $ E(2^{+}) $.

The overestimation of $ E(2^{+}) $ for certain nuclei in strongly deformed regions beyond $ N=100 $ can be attributed to several theoretical deficiencies. These include the inadequate treatment of tensor and spin-orbit couplings in Skyrme functionals, which lack complete tensor terms, and the "pairing collapse" phenomenon in the SHN region, for instance in ^246,248Cf and ^254,256Rf. Further related issues are the presence of non-negligible negative hexadecapole deformation or more complex shape effects such as rapid shape coexistence, γ-softness, and octupole correlations, particularly in Dy, Hf, Os, and Pt isotopes. While multi-parameter optimization is a long-term goal, targeted pairing-strength refinement remains a justified strategy. It sufficiently lowers $ E(2^{+}) $ toward experimental values while preserving the global predictive power of the EDF. Crucially, it avoids the structural overcorrection and physical contradictions associated with altering the volume-surface mixing, which can distort shell gaps, deformation, and the physical expectation near shell closures. By identifying and quantifying several critical sources of uncertainty in predictions of collective rotation and $ E(2^{+}) $, this work will contribute to the accuracy needed to interpret previously uncharacterized nuclei, thereby guiding future experiments.