Nuclear level density from relativistic density functional theory and combinatorial method

The origin of elements is one of the most fundamental scientific problems. Currently, the rapid neutron capture process (r-process) is believed to be responsible for the nucleosynthesis of about half of the elements heavier than Iron Burbidge et al. (1957); Kajino et al. (2019); Cowan et al. (2021). In r-process simulations, neutron capture rates of thousands of neutron-rich nuclei are required, and most of these nuclei lie far beyond the experimentally known region. Therefore, the r-process models rely on the theoretical predictions for the required neutron capture rates Larsen et al. (2019). At present, neutron-capture rates are usually calculated with the Hauser-Feshbach model Hauser and Feshbach (1952), in which the nuclear level density plays a crucial role.

Accurate prediction of the nuclear level density is challenging because of the exponential growth of the nuclear level density with the excitation energy, and the complex nuclear structure and dynamics such as shell structures, pairing correlations, and collective motions Ericson (1960); Bohr and Mottelson (1969). In a pioneering work Bethe (1936), Bethe proposed a level density formula based on the Fermi gas model with equidistant spacing of the single-particle levels near the Fermi level, in which pairing correlations and collective excitations were not considered. To include the pairing correlations, the conventional shifted Fermi gas model Newton (1956); Cameron (1958) is proposed by shifting the excitation energy according to the semiempirical pairing energies. By treating the energy shift as a free parameter, the back-shifted Fermi gas model Vonach et al. (1968); Huizenga et al. (1969); Vonach and Hille (1969) leads to a better description of the level density at low excitation energy. By using an empirical constant temperature formula, the Gilbert-Cameron model Gilbert and Cameron (1965) also improves the nuclear level density at low excitation energy. For better consideration of the pairing correlations, the Bardeen-Cooper-Schrieffer (BCS) theory is adopted in the generalized superfluid model Ignatyuk et al. (1979, 1993), in which the collective excitations are taken into account with the enhancement factors. By solving a pairing Hamiltonian, along with the single-particle levels in a Woods-Saxon well plus a spin-orbit interaction, the partition function and nuclear level density are calculated in Ref. Hung et al. (2017); Dang et al. (2017).

The conventional nuclear shell model exactly solves a many-body Hamiltonian within a given model space, and could in principle provide an exact prediction of the nuclear level density. However, the computational cost is highly demanding. The shell-model Monte Carlo approach Nakada and Alhassid (1997); Alhassid et al. (1999, 2007, 2015) solves this problem and has been applied to nuclei as heavy as lanthanides Alhassid et al. (2008); Özen et al. (2013, 2015), but it is limited to schematic nuclear Hamiltonian due to the sign problem. More realistic interactions can be employed in the moment method Sen’kov and Horoi (2010); Sen’kov et al. (2013); Sen’kov and Zelevinsky (2016), the stochastic estimation method  Shimizu et al. (2016); Chen et al. (2023), and the extrapolated Lanczos matrix approach Ormand and Brown (2020) but are limited to light and medium-mass nuclei only Fanto and Alhassid (2021). Self-consistent mean-field methods can be used to describe nuclear level density across the entire nuclear chart. Based on nuclear properties provided by mean-field calculations, the nuclear level density is usually calculated using either the statistical method or the combinatorial method. Many calculations of this kind have been reported, including the extended Thomas-Fermi approximation with Skyrme forces plus statistical method Kolomietz et al. (2018), the Skyrme Hartree-Fock-BCS plus statistical method Demetriou and Goriely (2001); Minato (2011), the Skyrme Hartree-Fock-Bogoliubov plus combinatorial method Goriely et al. (2008), and the temperature-dependent Gogny Hartree-Fock-Bogoliubov plus combinatorial method Hilaire et al. (2012). More recently, to better include correlations beyond the mean field, a new method based on the boson expansion of QRPA excitations has been proposed Hilaire et al. (2023).

During the past decades, the relativistic density functional theory (DFT) has been successfully applied to describe both ground-state and excited properties of atomic nuclei Meng (2016). For nuclear level density, the relativistic DFT has been used to calculate the level density with both the combinatorial method Geng et al. (2023) and the statistical method Zhao et al. (2020). It should be emphasized that the differences between the statistical method and the combinatorial method are fundamental. In the statistical method, also known as the Darwin-Fowler method and partition function method, the intrinsic level density is obtained with an inverse Laplace transform of the partition function with the saddle-point approximation Ericson (1960); Bohr and Mottelson (1969), and the collective excitations are taken into account by enhancement factors. In the combinatorial method, however, the states of single-particle excitations are exactly counted by expanding a generating function constructed from the single-particle levels. The combinatorial method provides the energy-, spin-, and parity-dependent nuclear level densities and describes their nonstatistical behaviors Goriely et al. (2022) which significantly influence the cross section predictions Goko et al. (2006). Non-relativistic calculations with the combinatorial method have been proposed Hilaire et al. (2001); Hilaire and Goriely (2006); Goriely et al. (2008); Hilaire et al. (2012) and the predictions have been used as input data in the TALYS code package Koning and Rochman (2012). The relativistic DFT calculations with the combinatorial method have been reported in Ref. Geng et al. (2023), while the pairing correlations are neglected. Furthermore, to remove the large fluctuations of state density at low excitation energy, a smoothing procedure is required. However, the conventional smoothing method described in Ref. Hilaire et al. (2001) faces difficulties in correctly describing total state densities at low excitation energy, leading to unsatisfactory reproductions of the corresponding experimental data. Therefore, it is necessary to improve the combinatorial method with a proper smoothing method.

In this work, we calculate the nuclear level density with the combinatorial method based on the relativistic DFT. The relativistic Hartree-Bogoliubov (RHB) theory is employed to calculate nuclear properties. The pairing correlations ignored in the previous investigation are taken into account at the mean-field level and in single-particle excitations as well. A more mathematically appropriate method, the Strutinsky method, is adopted to remove the large fluctuations of state density at low excitation energy. The nucleus $\mathrm{{}^{112}Cd}$ is taken as an example to show the results.

In the present work, the relativistic density functional PC-PK1 Zhao et al. (2010) is employed in the RHB calculation and, in the pairing channel, a finite-range separable pairing force with the pairing strength $G=728\ \mathrm{MeV\ fm^{3}}$ Tian et al. (2009) is adopted. The RHB equation is solved by expanding the quasi-particle wavefunctions in terms of a three-dimensional harmonic oscillator basis in Cartesian coordinates Nikšić et al. (2014), which contains 14 major shells. The obtained ground-state deformation of $\mathrm{{}^{112}Cd}$ is $\beta_{2}=0.145$.

In the calculation of nuclear level density, the cut-off of the excitation energy and angular momentum are taken as $10\ \mathrm{MeV}$ and $49\hbar$ respectively, and the energy discretization unit is taken as $\varepsilon_{0}=0.005,0.01,0.05\ \mathrm{MeV}$ respectively. The energy interval in the conventional smoothing method is $\delta U=0.2\ \mathrm{MeV}$ (Hilaire et al., 2001). The parameters in the Strutinsky method are chosen to be $\gamma_{0}=0.2\ \mathrm{MeV}$ and $M_{0}=1$.

The dependence of the total state densities $\rho$ given by Eq. (6) on the energy discretization unit $\varepsilon_{0}$ in the combinatorial method is illustrated in Fig. 1, where the total state densities $\rho$ of $\mathrm{{}^{112}Cd}$ are shown as functions of excitation energies. A dramatic dependency on $\varepsilon_{0}$ is found, i.e., the smaller $\varepsilon_{0}$, the more significant the fluctuations. It could be foreseen that if $\varepsilon_{0}$ were to be reduced to zero, the $\rho$ would become a combination of a series of delta functions lying on the exact excitation energies of the folded states. It is therefore understandable that, in Eq. (6), a small $\varepsilon_{0}$ does not lead to smooth total state densities. Consequently, a smoothing method is required to remove the large fluctuations and obtain smooth total state densities.

The smooth total state densities $\tilde{\rho}$, obtained with the conventional smoothing method and the Strutinsky method, are shown in Fig. 2. The total state densities $\rho$ are significantly smoothed by both methods. However, the Strutinsky method better presents the details at low excitation energy, in particularly the ground state. This would influence the nuclear level density. Moreover, smooth total state densities $\tilde{\rho}$ barely change when $\varepsilon_{0}$ is taken from $0.005\ \mathrm{MeV}$ to $0.01\ \mathrm{MeV}$, and even to $0.05\ \mathrm{MeV}$. Therefore, in the present work, we take $\varepsilon_{0}=0.01\ \mathrm{MeV}$.

The performances of different formulas of moments of inertia, i.e., rigid rotor, irrotational fluid, and the Inglis-Belyaev formula, in the nuclear level density prediction are compared in Fig. 3, together with the data from shape method Goriely et al. (2022) and low-lying levels in RIPL3 database Capote et al. (2009). The moments of inertia rotating around an axis perpendicular to the symmetry axis $\mathcal{J}_{\perp}$ calculated by these three formulas fulfill the relation $\mathcal{J}_{\perp}^{\text{irr}}<\mathcal{J}_{\perp}^{\text{ing}}<\mathcal{J}_{\perp}^{\text{rig}}$. Naturally, a larger MOI would lead to larger nuclear level densities. This is because a larger MOI leads to smaller rotational excitation energies as in Eq. (11), and thus leads to denser levels in rotation bands. In this sense, the rigid rotor assumption may lead to too large of MOI and the irrotational fluid assumption may lead to too small of MOI in the calculations of nuclear level densities. The Inglis-Belyaev formula provides a proper MOI and it reproduces the experimental data quite well. This conclusion remains true even after taking into account the Thouless-Valatin dynamical rearrangement contributions by enhancing the Inglis-Belyaev MOI by 30% Libert et al. (1999). Therefore, the Inglis-Belyaev formula is suggested to be used in the microscopic calculations of nuclear level densities and would be adopted in the following discussions.

Pairing correlations are certainly important in the predictions of nuclear properties, which would certainly influence the predictions of nuclear level densities. The nuclear level densities of $\mathrm{{}^{112}Cd}$ calculated with and without pairing correlations are presented in Fig. 4. It is found that the calculation with pairing correlations reproduces the experimental data for nuclear level densities quite well. In contrast, the calculation without pairing correlations provides nuclear level densities of about one order of magnitude higher than the ones with pairing correlations. The pairing correlations reduce the MOI and the density of incoherent ph states and, thus, lower the nuclear level density.

Finally, in Fig. 5, the nuclear level densities of $\mathrm{{}^{112}Cd}$ calculated by the combinatorial method improved with Strutinsky method and based on RHB with PC-PK1 and DD-PC1 are compared with the ones calculated by the statistical method based on Skyrme Hartree-Fock-BCS with MSk7 Demetriou and Goriely (2001), statistical method based on finite-temperature RHB with DD-PC1 Zhao et al. (2020), combinatorial method based on Skyrme Hartree-Fock-Bogoliubov with BSk14 Goriely et al. (2008), and combinatorial method based on temperature-dependent Gogny Hartree-Fock-Bogoliubov with D1M Hilaire et al. (2012). The result based on PC-PK1 in this work not only reproduces the nuclear level density at high excitation energy, but also well describes the fluctuations at low excitation energy ($\lesssim 2\ \mathrm{MeV}$). With DD-PC1, both statistical method and combinatorial method predict slightly higher results at high excitation energy and similar fluctuations at low excitation energy. The excitation energy for the first peak predicted by statistical method with DD-PC1 is different from the ones by combinatorial method with DD-PC1 and PC-PK1 because the collective effects are taken into account in different ways. In contrast, the D1M result well reproduces the experimental data at high excitation energy but fails at low excitation energy for the fluctuations. The MSk7 and BSk14 results both significantly overestimate the nuclear level density at high excitation energy and also fail to describe the fluctuations. The results in this work can be regarded as a benchmark of nuclear level density based on relativistic DFT, which well reproduces experimental data for nuclear level density of $\mathrm{{}^{112}Cd}$.

In summary, the nuclear level density is investigated with the combinatorial method based on relativistic density functional theory including pairing correlations for the first time. In the combinatorial method, the Strutinsky method effectively removes the large fluctuations at low excitation energy and smoothes the total state density. The microscopic Inglis-Belyaev formula provides a proper MOI and it reproduces the experimental data quite well. The calculation with pairing correlations well reproduces the experimental data, while the calculation without pairing correlations provides nuclear level densities of about one order of magnitude higher.

Taking $\mathrm{{}^{112}Cd}$ as an example, it is demonstrated that the nuclear level density based on the relativistic density functional PC-PK1 can reproduce the experimental data at the same level as or even better than the previous approaches. The successful descriptions of the nuclear level density with relativistic DFT pave a road for future developments of systematic calculations applied to the nuclear chart and further applications to neutron capture rates and r-process.