Nuclear level density from relativistic density functional theory and combinatorial method

Nuclear level density is calculated with the combinatorial method based on the relativistic density functional theory including pairing correlations. The Strutinsky method is adopted to smooth the total state density in order to refine the prediction at low excitation energy. The impacts of pairing correlations and moments of inertia on the nuclear level density are discussed in detail. 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.


I. INTRODUCTION
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 [1][2][3].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 [4].At present, neutroncapture rates are usually calculated with the Hauser-Feshbach model [5], 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 [6,7].In a pioneering work [8], 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 [9,10] 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 [11][12][13] leads to a better description of the level density at low excitation energy.By using an empirical constant temperature formula, the Gilbert-Cameron model [14] 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 [15,16], 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. [17,18].
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 [19][20][21][22] solves this problem and has been applied to nuclei as heavy as lanthanides [23][24][25], but it is limited to schematic nuclear Hamiltonian due to the sign problem.More realistic interactions can be employed in the moment method [26][27][28], the stochastic estimation method [29,30], and the extrapolated Lanczos matrix approach [31] but are limited to light and medium-mass nuclei only [32].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 [33], the Skyrme Hartree-Fock-BCS plus statistical method [34,35], the Skyrme Hartree-Fock-Bogoliubov plus combinatorial method [36], and the temperature-dependent Gogny Hartree-Fock-Bogoliubov plus combinatorial method [37].
More recently, to better include correlations beyond the mean field, a new method based on the boson expansion of QRPA excitations has been proposed [38].
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 [39].
For nuclear level density, the relativistic DFT has been used to calculate the level density with both the combinatorial method [40] and the statistical method [41].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 [6,7], 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 [42] which significantly influence the cross section predictions [43].Non-relativistic calculations with the combinatorial method have been proposed [36,37,44,45] and the predictions have been used as input data in the TALYS code package [46].The relativistic DFT calculations with the combinatorial method have been reported in Ref. [40], 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. [44] 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 112 Cd is taken as an example to show the results.

II. THEORETICAL FRAMEWORK A. The relativistic Hartree-Bogoliubov theory
The relativistic Hartree-Bogoliubov (RHB) theory is employed to provide a unified and self-consistent treatment of mean field and pairing correlations.The detailed formulism of the RHB theory can be seen in Refs.[47][48][49][50].In the RHB theory, one needs to solve self-consistently the RHB equation for nucleons, where ĥD is the single-nucleon Dirac Hamiltonian, λ τ is the Fermi energy (τ = n/p for neutrons and protons), ∆ is the pairing potential, U k and V k are the quasi-particle wavefunctions, and E k is the corresponding quasi-particle energy.The single-nucleon Dirac Hamiltonian ĥD reads where, The pairing potential ∆ is determined by where V pp is the pairing force and κ is the pairing tensor.The RHB equation (1) needs to be solved self-consistently because, the scalar, vector, and pairing potentials are coupled self-consistently with the densities and, in turn, the quasi-particle wavefunctions.

B. The combinatorial method
In the combinatorial method, the nuclear level density is calculated based on nuclear single-particle levels, masses, radii, deformations, and moments of inertia, and the detailed formulism can be seen in Refs.[36,37,44,45].The single-particle levels are used to calculate the number of incoherent particle-hole (ph) states as functions of the single-particle excitation energy, the angular momentum projection onto the symmetry axis, and the parity.
For a given excitation energy of the nucleus, the total number of the states is deduced by folding the incoherent ph states and the collective vibration states whose number is counted by introducing a generalized boson partition function [36].In principle, one can derive the total state density at any excitation energy by counting the total number of states.However, one should keep in mind that the excitation energies are continuous so there are infinitely many of them.Therefore, the numbers of the folded states are counted only at a series of equally spaced excitation energies which are expressed as integer multiples of an energy discretization unit ε 0 .The resultant total state densities ρ at each excitation energy are then given by Here, U is the excitation energy, K is angular momentum projection onto the symmetry axis, P is parity, and C(U, K, P ) is the number of the folded states in an energy interval ε 0 centered around U. It turns out that the obtained total state densities ρ in Eq. ( 6) strongly depend on the energy discretization unit ε 0 , and one cannot obtain a smooth total state density as a function of the excitation energy U even at very small ε 0 .Consequently, a smoothing method is required to obtain smooth total state densities against the excitation energy.
In the conventional smoothing method [44], the logarithm of the cumulated number of states is first obtained with ln N(U, K, P ) = ln and then the total state densities are obtained with ρ(U, K, P ) ≡ dN(U, K, P ) dU = N(U, K, P ) where N(U, K, P ) and d ln N(U, K, P )/dU are deduced from a linear interpolation, over an energy interval δU centered around U, of the ln N(U, K, P ) values.This method gives a smooth total state density against the excitation energy, but is found to present poor results at low excitation energy.
To improve the predictions at low excitation energy, the Strutinsky method [51][52][53][54][55], is adopted to directly smooth the total state density ρ given by Eq. ( 6) with ρ(U, where the function f (x) is constructed as f (x) = 1 √ π e −x 2 P (x), where P (x) is a generalized Here, γ 0 is the smoothing range, and M 0 is the order of the generalized Laguerre polynomial.The Strutinsky method can be justified with two important facts.Firstly, it would keep the smoothed function unchanged if smoothed again with the same procedure.Secondly, it strictly fulfills the conservation of the number of states.The second fact is essential for the accuracy of the calculation, especially at low excitation energy.
Based on the obtained smooth total state densities, the nuclear level densities can be deduced.For a deformed nucleus, collective rotation effects must be included, by building up rotational bands consistently on each of the folded states.The nuclear level density is, thus, given by [44]: where ρ is the smooth total state density, U is the excitation energy, J is angular momentum, K is angular momentum projection onto the symmetry axis, and P is parity.The δ (x) equals to 1 if x holds true and 0 otherwise, which restricts the rotational bands established on K = 0 states to the level sequences J = 0, 2, 4, ... for P = +, and J = 1, 3, 5, ... for P = −.The rotation energy E J,K rot is obtained with where J ⊥ is the moment of inertia of a nucleus rotating around an axis perpendicular to the symmetry axis.
Three formulas of moments of inertia (MOI) are examined in the work, i.e., the MOI of a rigid rotor the MOI of irrotational fluid and the microscopic Inglis-Belyaev formula [57, 58] Here, A is the mass number, and ζ denotes the axis of rotation.In Eq. ( 12) and ( 13), the nuclear mass M, quadrupole deformation parameter β 2 , and nuclear radius R are determined by solving the RHB equation.In Eq. ( 14), the occupation probabilities v i , and single-particle wavefunctions ψ i are obtained from quasi-particle wavefunctions by canonical transforming [54].The energies in Eq. ( 14) are E i = (ǫ i − λ) 2 + ∆ 2 i , in which λ is Fermi energy, ǫ i is single-particle energy, and ∆ i is energy gap.The summation in Eq. ( 14) runs over the proton and neutron single-particle states.

III. NUMERICAL DETAILS
In the present work, the relativistic density functional PC-PK1 [59] is employed in the RHB calculation and, in the pairing channel, a finite-range separable pairing force with the pairing strength G = 728 MeV fm 3 [60] 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 [61], which contains 14 major shells.The obtained ground-state deformation of 112 Cd is β 2 = 0.145.
In the calculation of nuclear level density, the cut-off of the excitation energy and angular momentum are taken as 10 MeV and 49 respectively, and the energy discretization unit is taken as ε 0 = 0.005, 0.01, 0.05 MeV respectively.The energy interval in the conventional smoothing method is δU = 0.2 MeV [44].The parameters in the Strutinsky method are chosen to be γ 0 = 0.2 MeV and M 0 = 1.The dependence of the total state densities ρ given by Eq. ( 6) on the energy discretization unit ε 0 in the combinatorial method is illustrated in Fig. 1, where the total state densities ρ of 112 Cd are shown as functions of excitation energies.A dramatic dependency on ε 0 is found, i.e., the smaller ε 0 , the more significant the fluctuations.It could be foreseen that if ε 0 were to be reduced to zero, the ρ 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 ε 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.

IV. RESULTS AND DISCUSSION
The smooth total state densities ρ, obtained with the conventional smoothing method and the Strutinsky method, are shown in Fig. 2. The total state densities ρ 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 ρ barely change when ε 0 is taken from 0.005 MeV to 0.01 MeV, and even to 0.05 MeV.Therefore, in the present work, we take ε 0 = 0.01 MeV. of inertia.The data from shape method [42] are displayed with black dots and the data from low-lying levels in the RIPL3 database [62] are displayed with solid black line.
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 [42] and low-lying levels in RIPL3 database [62].The moments of inertia rotating around an axis perpendicular to the symmetry axis J ⊥ calculated by these three formulas fulfill the relation J irr ⊥ < J ing ⊥ < J 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% [63].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.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 112 Cd.

V. SUMMARY
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 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.

FIG. 2 .
FIG. 2. (color online)The smooth total state densities ρ of 112 Cd obtained with the conventional smoothing method (a) and the Strutinsky method (b).

FIG. 3 .
FIG. 3. (color online) The nuclear level densities of 112 Cd based on different formulas of moments

FIG. 4 .FIG. 5 .
FIG. 4. (color online) Comparison of the nuclear level densities of 112 Cd calculated with (solid red line) and without pairing correlations (dashed red line).