Benchmarking mean-field approximations to level densities

Benchmarking mean-field approximations to level densities

Y. Alhassid,, G.F. Bertsch, C.N. Gilbreth, and H. Nakada Center for Theoretical Physics, Sloane Physics Laboratory, Yale University, New Haven, Connecticut 06520, USA
Department of Physics and Institute for Nuclear Theory, Box 351560
University of Washington, Seattle, Washington 98915, USA
Department of Physics, Graduate School of Science,Chiba University, Inage, Chiba 263-8522, Japan
Abstract

We assess the accuracy of finite-temperature mean-field theory using as a standard the Hamiltonian and model space of the shell model Monte Carlo calculations. Two examples are considered: the nucleus Dy, representing a heavy deformed nucleus, and Sm, representing a nearby heavy spherical nucleus with strong pairing correlations. The errors inherent in the finite-temperature Hartree-Fock and Hartree-Fock-Bogoliubov approximations are analyzed by comparing the entropies of the grand canonical and canonical ensembles, as well as the level density at the neutron resonance threshold, with shell model Monte Carlo (SMMC) calculations, which are accurate up to well-controlled statistical errors. The main weak points in the mean-field treatments are seen to be: (i) the extraction of number-projected densities from the grand canonical ensembles, and (ii) the symmetry breaking by deformation or by the pairing condensate. In the absence of a pairing condensate, we confirm that the usual saddle-point approximation to extract the number-projected densities is not a significant source of error compared to other errors inherent to the mean-field theory. We also present an alternative formulation of the saddle-point approximation that makes direct use of an approximate particle-number projection and avoids computing the usual three-dimensional Jacobian of the saddle-point integration. We find that the pairing condensate is less amenable to approximate particle-number projection methods due to the explicit violation of particle-number conservation in the pairing condensate. Nevertheless, the Hartree-Fock-Bogoliubov theory is accurate to less than one unit of entropy for Sm at the neutron threshold energy, which is above the pairing phase transition. This result provides support for the commonly used “back-shift” approximation, treating pairing as only affecting the excitation energy scale. When the ground state is strongly deformed, the Hartree-Fock entropy is significantly lower than the SMMC entropy at low temperatures due to the missing contribution of rotational degrees of freedom. However, treating the rotational bands in a simple model, we find that the entropy at moderate excitation energies is reproduced to within two units, corresponding to an error in the level density of less than an order of magnitude. We conclude with a discussion of methods that have been advocated as beyond the mean field approximation, and their prospects to ameliorate the issues we have identified.

I Introduction

Nuclear level densities are important input in the theory of low-energy nuclear reactions. In situations where the reactions cannot be studied in the laboratory, the need for a reliable theory is evident. The calculation of level densities in the presence of correlations is a challenging problem. Part of the problem is the complexity and inadequate knowledge of the nuclear Hamiltonian. But even with a known Hamiltonian, it is often necessary to make approximations based on mean-field theory whose accuracy is not well understood. Our goal here is to assess these mean-field approximations by taking the Hamiltonian as known and testing them against a theory that is accurate up to well-controlled statistical errors.

Most treatments of the level density of heavy nuclei start from a mean-field theory in the form of the finite-temperature Hartree-Fock (HF) or the Hartree-Fock-Bogoliubov (HFB) approximation. The finite-temperature theory is derived from a variational principle based on the grand potential as a function of an uncorrelated trial density go81 (); ta81 (). These approximations have been widely applied and taken as a starting point for more sophisticated theories ke81 (); va83 (). It would thus be useful to validate them by comparing with a more accurate method. The auxiliary-field Monte Carlo method, known as the shell model Monte Carlo (SMMC) in nuclear physics la93 (); al94 (), fulfills this role. Given the Hamiltonian and a model space, the only inaccuracy is a controllable statistical error associated with the Monte Carlo sampling. As a finite-temperature method, SMMC is particularly suitable for the calculation of level densities na97 (); al99 ().

The SMMC starts from a Hamiltonian that is somewhat restricted but which reproduces all the long-range correlations in a realistic way. These include deformations, pairing gaps and low-energy collective excitations. It is thus well-suited to provide a benchmark for testing the validity of the mean-field treatments of nuclear thermal and statistical properties. We note that the SMMC applies to all nuclei irrespective of whether they are deformed or spherical and independently of the existence of a mean-field pairing condensate. In contrast, the formulas for level densities in HF and HFB depend on the character of the mean-field solution.

In Sec. II we summarize the statistical and thermodynamic tools we will use for the analysis and comparisions. In Sec. III we present in some detail the results of the SSMC calculations for two nuclei, namely Sm and Dy. The first is a spherical nucleus having a strong pairing condensate in HFB. The second is a well-deformed nucleus with a weak pairing condensate; here the HF is an appropriate mean-field approximation. In Sec. IV we present the HF approximation and its results for Dy, while in Sec. V we discuss the HFB approximation and its results for Sm. In Sec. VI we summarize our findings regarding the accuracy of the HF and HFB approximations with some remarks on prospects for extending the mean-field approach to level densities. Finally, the data files from the SMMC, HF and HFB calculations together with the codes to apply them are provided in the Supplementary Material depository accompanying this article.

Ii Tools of statistical theory

ii.1 The canonical entropy

A good meeting point for comparing different statistical theories is the canonical entropy function , i.e., the entropy of the canonical ensemble of states having fixed numbers of protons and neutrons and at inverse temperature . In the SMMC, the quantity that is most directly computed is the thermal energy of the canonical ensemble. The entropy can then computed by integrating the thermodynamic relation

 dSc=βdEc (1)

as

 Sc(β)=Sc(0)−∫Ec(0)Ec(β)β′dEc. (2)

This allows one to calculate the partition function from . Alternatively, can be calculated directly from by integrating the relation

 dlnZc=−Ecdβ, (3)

and the canonical entropy is then calculated from . Finally, the density of states at given energy and particle numbers is obtained from by an inverse Laplace transform, carried out in the saddle-point approximation

 ρ(E,Np,Nn) = 12πi∫i∞−i∞dβ′eβ′EZc (4) ≈ (2π∣∣∣∂E∂β∣∣∣)−1/2eSc(β),

where in the above expression is determined as a function of from the saddle-point condition

 E=−∂lnZc∂β=Ec(β). (5)

The entropy of a system whose Hamiltonian is defined in a finite-dimensional model space satisfies a sum rule obtained from (2) in the limit

 ∫E(β=0)E(β=∞)βdE=S(0)−S(∞). (6)

For a finite-dimensional model space, the entropy at both end points and is finite and can be determined analytically. In particular, for even-even nuclei, the entropy at zero temperature must be zero.

In the configuration-interaction (CI) shell model, we have a certain number of nucleons in single-particle model spaces of dimensions giving a canonical entropy at of

 Sc(0)=Sp(0)+Sn(0)=ln(DpNp)+ln(DnNp). (7)

We have found the sum rule in Eq. (6) useful for testing the computer codes we have employed in this study, and also for setting end points on the entropy plots we show later.

ii.2 The grand canonical entropy

The finite-temperature HF (FTHF) and finite-temperature HFB (FTHFB) approximations are defined in the framework of the grand canonical ensemble which depends on the additional independent variables (), related to the chemical potentials by .

The grand canonical entropy , when expressed as a function of the energy and the average number of particles of type in the grand canonical ensemble, satisfies

 dSgc=βdEgc−∑i=p,nαidNi,gc, (8)

and can be calculated as in Eq. (2) for fixed , i.e., when the grand canonical energy is expressed as a function of and the given particle numbers .

The grand canonical entropy at is also a simple function of the dimension of the single-particle shell-model space in the CI shell model, since all correlations disappear in that limit. Choosing values for (where are chemical potentials) to produce average particle numbers , we have

 Sgc(0,Np,Nn)=−∑i=p,nDi[filnfi+(1−fi)ln(1−fi)], (9)

where are occupation factors for .

The grand canonical partition satisfies

 dlnZgc=−Egcdβ+∑iNi,gcdαi. (10)

The Legendre transform of with respect to is a function of and defined by and satisfies

 dln~Zgc=−Egcdβ−∑iαidNi,gc. (11)

Thus we can alternatively calculate the grand canonical entropy by integrating with respect to at fixed to determine and then find the entropy from .

The state density is related to the partition function by the three-dimensional inverse Laplace transform

 ρ(E,Np,Nn)=1(2πi)3∫i∞−i∞dβ∫i∞−i∞dαp∫i∞−i∞dαn ×eβE−αpNp−αnNnZgc(β,αp,αn). (12)

Normally one carries out the integration over all three variables in the three-dimensional (3-D) saddle-point approximation, resulting in the formula for the state density BM ()

 ρ(E,Np,Nn)=1(2π)3/2∣∣∣∂(E,Np,Nn)∂(β,αp,αn)∣∣∣−1/2eSgc, (13)

where the values of are determined from by the saddle-point conditions

 E =−∂lnZgc∂β=Egc(β,αp,αn), Ni =∂lnZgc∂αi=Ni,gc(β,αp,αn). (14)

Using Eqs. (II.2), the Jacobian in (13) can be written as the determinant of the matrix of second derivatives of the logarithm of the grand canonical partition function with respect to .

ii.3 Ensemble reduction

To compare the mean-field entropies to the canonical SMMC entropies we need to reduce the grand canonical ensemble of the mean-field formalism to the canonical ensemble. For approximate theories such as the HF and HFB, the only consistent way (in the sense of Appendix II) to carry out the reduction is the variation-after-projection method (VAP), but this is difficult andm as far as we know, has never been put into practice for calculating level densities in heavy nuclei. We will therefore only consider simpler reduction methods, recognizing that they cannot be free from ambiguity.

A straightforward way to determine a canonical entropy is to separate the 3-D saddle-point integration Eq. (II.2) into two steps, integrating first over the chemical variables . This yields the following expression for the integrand of the integration:

 ζ−1Zgc(β,αp,αn)eβE−∑iαiNi, (15)

where

 ζ=2π∣∣∣∂(Np,Nn)∂(αp,αn)∣∣∣1/2. (16)

and are determined by the 2-D saddle-point conditions (). Comparing with the integrand in Eq. (4), we can identify the approximate canonical partition function as

 lnZc≈lnZgc−∑iαiNi−lnζ. (17)

If we carry out in a second step the integration of Eq. (15), we obtain an expression of the same form as in Eq. (4) where the approximate canonical entropy is given by

 Sc≈Sgc−lnζ. (18)

The expression we find for the state density is equivalent to Eq. (13) of the 3-D saddle-point approximation det ().

However, the above result does not take into account the variation of the prefactor with respect to . If we consider this dependence explicitly when we perform the integration, the saddle-point condition that determines in terms of the energy becomes , where is an approximate canonical energy given by

 Eζ=Egc−δE,withδE=−dlnζdβ. (19)

The level density is then given by the canonical form (4), where the canonical entropy is

 Sc(β,Np,Nn)=Sgc−lnζ+βδE. (20)

A simple model is presented in Appendix I, showing that the saddle-point shift in Eq. (19) improves the accuracy of the calculated and .

ii.4 Discrete Gaussian model and the particle-number fluctuation

Another source of error may arise from the Gaussian approximation in the 2-D saddle-point integration used to derive Eq. (16). For example, suppose that the grand canonical partition function were dominant by a single nuclide in some range of . Then we could approximate . Treating this in the saddle-point Gaussian integration gives (since are independent of ), rather than the correct answer of . The problem can be repaired by recognizing that and are discrete integers, not continuous variables. We calculate the matrix as before but now we calculate as a discrete sum over particle numbers ,

 ζ=∑N′i,N′jexp(−12∑i,j∂N∂α∣∣∣−1ij(N′i−Ni)(N′j−Nj)). (21)

Expression (21) for reduces to the saddle-point result (16) in the limit when the r.h.s. of (16) is large, but has the advantage that it is always larger than and it approaches when the r.h.s. of (16) goes to .

To see more physically how the approximation (21) works, consider the case when there is only one type of particles, say neutrons, and the Hamiltonian used in the Gibbs density operator is independent of and . The required derivative is then

 ∂2lnZgc∂α2n=∂Nn∂αn=⟨(ΔNn)2⟩, (22)

where is the neutron-number fluctuation in the grand canonical ensemble. Carrying out the Gaussian saddle-point integration, in the 2-D saddle-point approximation (15) becomes

 ζ−1n=(2π⟨(ΔNn)2⟩)−1/2. (23)

is just the ratio of states with particle number to the total number of states in an ensemble in which the particle number is distributed as a discrete Gaussian

 PN′n=ζ−1ne−(N′n−Nn)2/2⟨(ΔNn)2⟩ (24)

in the limit that .

The finite-temperature mean-field approximation also provides a many-particle density matrix, so that the particle-number fluctuation can be calculated directly from this density matrix (see Eqs. (37) and (47) below). In fact, this direct calculation is much easier to carry out than calculating numerically the matrix of second derivatives of the logarithm of the partition function. However, because the canonical reduction is not carried out in a variational way, the two methods are not guaranteed to give the same answer. In the sections below, we will examine and compare both methods of carrying out the canonical reduction.

In the finite-temperature mean-field approximations we use here, the off-diagonal particle-number correlations vanish for and factorizes into two separate factors for protons and neutrons

 ζ=ζpζn,whereζi=∑N′ie−(N′i−Ni)2/2⟨(ΔNi)2⟩. (25)

in Eq. (25) can be considered as the partition function which describes the fluctuations in the number of particles of type . The reduction from the grand canonical to the canonical partition function is then given by

 Zgc(β,αp,αn)e−∑iαiNi≈Zc(β,Np,Nn)ζpζn. (26)

Relation (26) describes the factorization of the grand canonical partition function into a canonical partition function and particle-number fluctuation partition functions. It is exact in the simple model presented in Appendix I.

In summary, the saddle-point approximation breaks down when . However, in the discrete Gaussian model (Eqs. (25) and (26) or Eq. (21)) always satisfies and can be used even when the particle-number fluctuation is small. We will demonstrate the improvement to the saddle-point formula in this limit in the case of Dy  (Sec. IV.2) where pairing correlations are weak.

ii.5 Spin-parity projected level density

The ultimate goal is to calculate the spin-parity projected densities , defined as the number of levels of given angular momentum and parity per unit energy, not counting the magnetic degeneracy of the levels. The spin-dependent level densities can be calculated through an angular momentum projection. The present paper is mainly focussed on the total state density and we will not examine in details spin-projection methods. However, to make at least a tentative comparison of the level density at the neutron resonance threshold we will calculate them taking a simplified model for the spin-parity projection. We follow common practice and assume the spin distribution is Gaussian in the three components of the angular momentum vector  er60 (). Then the fraction of levels having angular momentum is given by

 PJ≈√12πJ+12σ3exp⎛⎝−(J+12)22σ2⎞⎠. (27)

Here a pre-exponential factor of , arising from the three-dimensional volume element of , is reduced to the first power of by dividing out the magnetic degeneracy factor. The parameter , known as the spin cutoff parameter, is estimated from the second moment of

 σ2=⟨J2z⟩. (28)

The normalization condition of in (27) is . Assuming equal positive- and negative-parity level densities (usually justified at the neutron binding energy), we have

 ρJπ(E)≈12PJρ(E). (29)

Iii SMMC results

The SMMC method is formulated in the framework of a CI spherical shell-model Hamiltonian. The CI Hamiltonians shown to be amenable to Monte-Carlo sampling contain one- plus two-body operators, with the two-body part restricted to interactions that have a “good sign” in the grand canonical formulation sign (). In finite nuclei, the method is implemented with particle-number projection or94 () for both protons and neutrons, and the calculated observables are the expectation values in the canonical density matrix. In particular, we consider here the nuclei Sm and Dy. The nucleus Sm is an example of a heavy spherical nucleus whose ground state has significant correlation energy associated with pairing, while the nucleus Dy  has a deformed ground state with weak pairing.

The parameterization of the Hamiltonian and other aspects of the SMMC calculation have been published elsewhere al08 (); oz13 (). For reference, Table I and its caption describes the model space employed in the calculations.

The canonical energy and the mean square angular momentum at fixed numbers of protons and neutrons are calculated directly in SMMC as a function of . Table 2 shows the SMMC energies at (i.e., the infinite temperature limit) and at high extrapolated to infinity. The energy at is largely determined by the one-body part of the Hamiltonian in the grand canonical ensemble. There is no contribution from the direct component of the interaction but there is a small contribution from the exchange terms.

The variation of the excitation energy with is shown in Fig. 1 using a logarithmic scale for the energy. The excitation energy of Dy is higher than that of Sm  from to and is then lower up to . The higher excitation energy in Sm near is likely due to the collapse of strong pairing in that nucleus. Similarly, the higher Dy excitation energy at may be ascribed to the loss of deformation energy in that temperature region.

We also need the ensemble averages to calculate the spin-dependent level densities. These are shown in Fig. 2 for Sm  and Dy. The higher values of of Dy  at high are largely due to its deformation and its low-lying first excited state at MeV. In contrast, for Sm decreases dramatically at high , as expected for a nucleus with a ground state and a gap of MeV to the first excited state. At low , the remaining enhancement for Dy is due to its larger number of active valence nucleons in the model space. The errors shown in Fig. 2 are statistical errors from the Monte Carlo sampling.

We next apply Eq. (2) to compute the canonical SMMC entropy. We start from with the initial value of the canonical entropy given by Eq. (7) and use the relation

 ∫Ec(0)Ec(β)β′dEc=−βEc(β)+∫β0Ec(β′)dβ′, (30)

where is the canonical thermal energy calculated in SMMC as the thermal expectation value of the Hamiltonian. The results are shown in Fig. 3, with the main figure showing the low to intermediate values of , and the inset showing the large values of .

The Sm  and Dy  entropies are nearly equal for in the range MeV. We do not know any obvious reason why that should be the case. At higher values of , shown in the inset, one observes a difference between the two nuclei. For Sm at MeV the entropy is essentially zero, as expected at low temperatures for even-even nucleus with a pairing gap. On the other hand, the entropy of Dy remains a couple of units higher than zero up to at least  MeV. This is due to the low excitations associated with the ground-state rotational band, together with the weak pairing in this nucleus.

The state densities calculated from Eq. (4) are shown in Fig. 4. As expected, the state density is higher for Dy  than Sm  at low excitation energy. Interestingly, they become much closer at excitation energies in the range 5-10 MeV.

As a check on the quality of the CI Hamiltonian, we compare the calculated level densities with the experimental -wave resonance spacings , measured at an energy that corresponds to the neutron separation energy. These are calculated from the spin-parity dependent level density (29) and (27), taking the spin cutoff parameter from Eq. (28). The results are shown in Table III. The good agreement gives us confidence that the Hamiltonian is realistic enough to provide useful tests of the HF and HFB approximations.

Iv The finite-temperature HF approximation

We first consider the FTHF approximation. It is derived by minimizing the grand potential in terms of a trial uncorrelated many-particle density matrix. Such a density is uniquely characterized by its one-body density matrix , where label single-particle orbitals in the model space.

For simplicity, we consider only one type of particles and the results are easily generalized to both protons and neutrons. At the FTHF minimum, the one-body density satisfies the self-consistent equation

 ϱ=1eβhϱ−α+1, (31)

where is the one-body HF Hamiltonian, expressed respectively in terms of the one- and two-body matrices and of the configuration-interaction shell model Hamiltonian. The value of at the HF minimum is given by

 βΩHF=−lnZHF=βEHF−SHF−αNHF, (32)

where is the HF approximation to the grand canonical partition function. The thermal HF energy is calculated as

 EHF=tr(tϱ)+12tr(ϱvϱ), (33)

the entropy is given by

 SHF = −tr(ϱlnϱ)−tr[(1−ϱ)ln(1−ϱ)] (34) = −∑λfklnfk−∑k(1−fk)ln(1−fk),

and the average number of particles is computed as

 NHF=trϱ. (35)

The occupation probabilities are the usual Fermi-Dirac occupations where are the single-particle HF energies at temperature .

It is not obvious from Eqs. (32), (33) and (35) that satisfies the thermodynamic derivative relations (II.2) due to the dependence of the HF single-particle energies and density on and . Nevertheless, the contributions to the derivatives that originate in the implicit dependences on and vanish. In particular, the quantities defined in the HF theory satisfy

 −∂lnZHF∂β=EHF,∂lnZHF∂αi=Ni,HF. (36)

The proof is provided in Appendix II. The validity of these relations in FTHF guarantee that the grand canonical HF entropy satisfies the relation at fixed average particle numbers , and thus we can compute the entropy in the same way as we did in the SMMC.

iv.1 Approximate canonical projectors in FTHF

The canonical partition function can be approximated either in the saddle-point approximation of Sect. II.3 or in the discrete Gaussian model of Sect. II.4. However, the connection to particle-number fluctuations is more tenuous because the second derivative expressions for in terms of particle-number fluctuations such as Eq. (22) no longer holds for . Nevertheless, it is interesting to compare computed with the full matrix to that computed from the particle-number fluctuations. We make such a comparison for Dy in the next section.

iv.2 Application to 162Dy

Here we discuss the strongly deformed nucleus Dy. The pairing in this nucleus is weak, so that the FTHF is the appropriate mean-field theory.

iv.2.1 Number partition function ζ

We first examine the number partition function obtained from the approximations we presented in Sec. II. The diagonal particle-number fluctuations in FTHF are given by

 ⟨(Δ^Ni)2⟩=tr[ϱi(1−ϱi)] (37)

and the off-diagonal ones are zero. The corresponding obtained from Eq. (25) is shown as the dashed line in Fig. 5. This is compared to calculated from Eq. (21) (using the matrix ), shown as the solid line. They differ by less than 10%, except for a tiny region near the spherical-to-deformed phase transition. There the calculated from the Jacobian of diverges (when approached from the deformed side). Thus, it appears to be a very good approximation to calculate in term of the individual particle-number fluctuations as in Eq. (25). In the next section, we shall see that this is not the case for the FTHFB theory in the presence of strong pairing correlations. In any case, we will use Eq. (21) in the results shown below.

iv.2.2 Thermal excitation energy

The range of FTHF energies as a function of is shown in the fourth column of Table II. The high-temperature limit is very close to the SMMC value, since all correlation energies disappear in that limit. At the other limit of large , at or near the ground state, the SMMC energy is about 3.5 MeV lower than its HF value. The HF ground-state energy has the correlations associated with the static deformation, but is missing the rotational energy and other correlation effects. It can also be seen from the fifth column of Table II that the HFB approximation hardly lowers the energy.

The excitation energy is shown in Fig. 6, comparing its HF value (solid line) with the SMMC thermal energy (solid circles) from Fig. 1. The HF density is spherical for and becomes deformed above that value. The energy of the spherical HF solution at higher values of is shown in Fig. 6 by the dashed-double dotted line. We observe that the HF energy has a cusp at the onset of deformation. Extrapolating the energy of the spherical solution to large , we find that the deformed ground-state solution is lower in energy than the spherical solution by 11.4 MeV. The SMMC energy, shown by the solid circles, is remarkably close to the HF energy in the region MeV except in the vicinity of the shape phase transition. The cusp in the HF energy function disappears in the SMMC energy, leaving no trace of a shape phase transition. The fact that mean-field theory over-emphasizes phase transitions in finite systems is well-known in nuclear theory; see, e.g., in Refs. es93, ; ro98, ; ka06, ; ga13, . The inset in Fig. 6 shows at higher values of using a finer energy scale. We observe that at large the grand canonical HF energy (dashed line) overestimates the excitation energy but that the approximate canonical energy of Eq. (19) (solid line) is in overall good agreement with the SMMC excitation energy.

iv.2.3 Entropies

The grand canonical HF entropy is shown in Fig. 7 as a function of (dashed line) and compared with the approximate canonical entropy (20) with from the discrete Gaussian formula Eq. (21) (solid line). At , the grand canonical HF entropy is larger than the canonical entropy due to particle-number fluctuations. The entropies at large values of are shown in the inset. The grand canonical HF entropy vanishes in the limit , as expected. However, the saddle-point canonical entropy calculated from Eqs. (15) and (16) increases at large values of (dotted line in the inset), indicating the breakdown of the saddle-point approximation to the particle-number projection. In contrast, the discrete Gaussian treatment (21) gives an entropy that approaches zero at high , thus satisfying the sum rule Eq. (6). For moderate and small values of , the entropy (20) of the discrete Gaussian model (21) essentially coincides with the saddle-point canonical entropy. As noted earlier, the SMMC entropy remains nonzero in the range . We examine this further in the next paragraph.

To compare the projected HF and SMMC entropies in more detail, we have replotted them in Fig. 8 with some additional information. Both curves start at the same value at because the model spaces are identical. In the limit of large , the SMMC entropy does not approach zero as fast as the canonical HF entropy. This is because the SMMC entropy includes a contribution from the ground-state rotational band, most visible at  MeV. We can estimate this contribution as follows. The moment of inertia of the ground-state band of Dy was determined to be  MeV by fitting the low-temperature SMMC values of to  al08 (). For  MeV, , and we can treat the rotational motion classically. The classical partition function of the rotational band is given by . We can then calculate its entropy from , where is the free energy of the ground-state rotational band. We find

 Srot=1+ln(Igsℏ2T). (38)

This contribution is described by the dotted line in the inset of Fig. 8, and is in good agreement with the SMMC entropy at large with no adjustable parameters.

We also show in Fig. 8 the canonical entropy of the spherical HF solution (dashed-double dotted line). This entropy approaches a finite non-zero value in the limit (see inset) because of the large degeneracy of the spherical HF solutions at . There are 2 valence protons in the orbitals and 6 valence neutrons in the orbitals, leading to a highly degenerate ground state with a canonical entropy of . The grand canonical HF entropy in this limit is even larger. It can be calculated assuming the uniform filling of the valence degenerate orbitals in the limit of the HF approximation. The corresponding formula has the form of Eq. (9), where , and , and gives a grand canonical entropy of .

iv.2.4 Angular momentum fluctuations

In HF, the variance of the angular momentum components () can be calculated using Wick’s theorem as

 ⟨(ΔJq)2⟩=⟨^J2q⟩−⟨^Jq⟩2=tr[jq(1−ϱ)jqϱ], (39)

where is the matrix representing in the single-particle space. When the HF equilibrium ensemble is axially symmetric around the -axis, is invariant under rotations around the axis and . Assuming time-reversal invariance, we also have . It follows that the variances of the angular momentum components are the same as the mean square moments. In Fig. 9, we compare the HF mean square moments of and (solid lines) with the SMMC moments (solid circles) in Dy . The HF mean square moments of and coincide above the shape transition temperature, where the HF solution is spherical. However, at large values , the HF mean square moment of is much larger than the respective moment of . Since the deformed intrinsic ground state has good , approaches zero in the limit , while remains finite and large in this limit. We also show by dashed line for the spherical HF solution.

iv.2.5 State density and level spacing

In Fig. 10 we show the HF density vs.  in the saddle-point approximation (4) (solid line) (where the approximate canonical entropy and energy include the correction in Eqs. (20) and (19) respectively), and compare it with the SMMC state density (solid circles). The kink in the HF density at MeV signifies the shape transition from a deformed to spherical shape. At lower excitation energies, the HF state density underestimates the SMMC values; the SMMC density includes a contribution from rotational bands that are built on top of intrinsic states, and are not captured in the HF approximation. Above the shape transition energy, the equilibrium shape is spherical and no longer supports rotational bands. The HF density is then very close to the SMMC density.

We can try to repair the HF approximation by recognizing that the each of the deformed HF configurations represents a band bj74 (). The angular momentum corresponds to the -quantum number of the band. Assuming that is Gaussian distributed, the -dependent HF state density can be expressed

 ρK(E)=PKρHF(E). (40)

where

 PK=1(2πσ2)1/2exp(−K22σ2). (41)

and

 σ2=⟨J2z⟩. (42)

For each positive there will be a rotational band with . For the sequence may skip odd or even values. For a complete treatment of the band model one next introduces a moment of inertia of the band to calculate the -dependent level density. However, we do not wish to introduce new parameters that take us beyond the HF theory so we assume that all the levels in a band are degenerate. The level density is then

 ρJπ(E)≈12J∑K=0ρK(E), (43)

treating the bands the same as the others. The factor of is for parity projection. The resulting average resonance spacing at the neutron threshold of MeV for Dy  is reported in Table III. It underestimates the SMMC value by a factor of . This is a substantial disagreement; in our view uncovers a serious problem with the HF theory of level densities in deformed nuclei.

iv.2.6 Independent-particle model

A common approximation in the calculation of state densities is to take the HF ground-state mean field, and assume the excited states can be calculated in the independent-particle model (IPM) with single-particle energies given by that potential field. For an axially deformed nucleus such as Dy, the single-particle HF levels come in doubly degenerate time-reversed pairs and for an even number of particles the ground state is non-degenerate, so that the entropy is zero. It is rather easy to carry out the exact paritcle-number projection in the IPM be14 (), so we will use it in the comparison.

In Fig. 11, we compare the IPM state density with the FTHF density determined by using Eq. (20) for . They agree very well at low excitation energy. At the neutron separation energy, MeV, the IPM state density is lower than the FTHF density by less than a factor of 2. However, the discrepancy increases with excitation energy, exceeding one unit entropy beyond  15 MeV excitation energy. At the shape phase transition, the IPM underestimates the HF density by more than two orders of magnitude. We conclude that IPM is a good approximation at excitation energies that are small compared to the shape transition energy, but not near and above this transition energy.

V The finite-temperature HFB approximation

The HFB is the preferred mean-field approximation for nuclei exhibiting strong pairing correlations. Like the FTHF, the FTHFB is based on the grand canonical ensemble. However, unlike the FTHF, the simple approximate particle-number projection onto a canonical ensemble is not expected to be a good approximation at low temperatures. This will be evident as we go through the steps to calculate the level density of Sm, starting from the HFB energy function . To simplify the notation, we consider only one type of particles as we did in the FTHF. The HFB thermal energy is expressed in terms of the normal and anomalous densities as

 EHFB=tr(tϱ)+12tr(ϱvϱ)+14tr(κ†vκ). (44)

The HFB entropy is given by

 SHFB=−∑kfklnfk−∑k(1−fk)ln(1−fk), (45)

where

 fk=11+eβEk (46)

are the quasi-particle occupations expressed in terms of the quasi-particle energies .

The HFB partition function satisfies relations similar to Eqs. (36), and thus the entropy can be computed from the energy function by an integral similar to Eq. (2). All the expressions regarding the level density in Sec. IV carry over to the HFB approximation, except that the particle-number variance in Eq. (37) is now calculated using all three Wick contraction terms in the expectation value . This leads to an additional contribution from the anomalous density

 ⟨(Δ^N)2⟩=∑kϱkk−∑klϱ2kl+∑kl|κkl|2. (47)

v.1 Application to 148Sm

The mean-field ground state of Sm is spherical and has a substantial pairing condensate. Thus FTHFB is the proper mean-field theory for this nucleus. The pairing correlation energy, defined as the difference between the HF and HFB ground-state energies is 1.82 MeV (see Table II). Also shown in Table II is the correlation energy of the SMMC ground state with respect to the HFB solution (i.e., the difference between the HFB and SMMC ground-state energies). This correlation energy is labeled “missing” in the table and is about MeV for Sm. The pairing transitions for protons and neutrons occur at MeV ( MeV) and MeV ( MeV), respectively.

We first examine , the factor used to convert the grand canonical quantities to canonical in the HFB. Naively, we may try using the HFB particle-number fluctuation in Eq. (47) in Eq. (25). The resulting is shown as the dashed line in Fig. 12. This approximation is bound to fail at low temperatures because the HFB ground state has a nonphysical particle-number fluctuation. We have also calculated from Eq. (21) using the full matrix suggested by the 3-D saddle-point approximation. This gives a better result for MeV, as may be seen by the solid line in Fig. 12. We will therefore use this method for the canonical quantities calculated below.

v.1.1 Excitation energies

The Sm thermal excitation energy function vs.  in FTHFB is compared with the SMMC results in Fig. 13. The HFB energy is shown by the dashed line. The excitation energies at are nearly equal, differing only by the missing correlation energy. The two kinks in the HFB curve, visible in the inset, are associated with the neutron and proton pairing phase transitions. The HFB excitation energy is higher than that of the SMMC, as to be expected from the higher limiting entropy of the grand canonical ensemble. We next calculate the approximate canonical projection of the HFB energy using Eq. (21). The result is closer to the SMMC for high , We note that for MeV(i.e., for temperatures above the pairing transitions) , the absolute HF and HFB energies coincide, so that the HF excitation energy is lower than the HFB excitation energy by exactly the amount of pairing correlation energy in the ground state.

v.1.2 Entropies

The grand canonical HFB entropy (dashed line) and SMMC entropy (open squares) functions in Sm are shown in Fig. 14 vs. . Their absolute values are set by integrating from the point, where their respective values are given in Table I. Both entropy functions approaches zero at large , confirming the sum rule (6). The neutron and proton pairing transitions are also visible as kinks in the HFB entropy curve Cv ().

We also show in Fig. 14 the approximate canonical HFB entropy in Eq. (20) where is given by Eq. (21) in the discrete Gaussian model. Since the particle-number fluctuations in FTHFB remain relatively large even at low temperatures, similar results for are found in the saddle-point approximation Eq. (16).

This approximate canonical entropy coincides with the SMMC at low values of and overestimates the SMMC entropy around MeV , i.e., in the vicinity of the proton pairing transition. At larger values of , for which a nonzero pairing condensate exists, the approximate canonical entropy is in overall agreement with the SMMC entropy up to MeV but at lower temperatures it becomes negative with a value of about -2 at MeV when the system reaches its HFB ground state. We note that if were to be calculated using the HFB particle-number fluctuations, the large- entropy would have been even more negative at units. A negative entropy at zero temperature is unphysical since there is only one state in the ensemble at zero temperature and the entropy should be zero. The HFB ground state violates particle-number conservation; the probability of having the proper proton and neutron numbers for Sm at MeV is only , hence the unphysical negative entropy at low temperatures.

To summarize the results of this section, we have not found a simple acceptable procedure to project the HFB onto a canonical ensemble if we require the correct entropy at high temperature and an error of less that one unit at low temperatures. We will comment further on this situation in the conclusion.

v.1.3 Angular momentum fluctuations

Eq. (39) which describes the angular momentum fluctuations in the FTHF approximation, has an additional contribution in FTHFB from a contraction in Wick’s theorem that involves the anomalous density

 ⟨(ΔJq)2⟩ = ⟨J2q⟩−⟨Jq⟩2 (48) = tr[jq(1−ϱ)jqϱ]−tr[jqκj∗qκ∗].

The additional contribution is negative, leading to a reduction in the mean-square moment of the angular momentum. This is just what one would expect as an effect of the pairing correlations.

In Fig. 15 we show the angular momentum fluctuations in Sm as a function of . The HFB solutions in Sm is spherical at all temperatures so all directions are equal and we only need examine one of them. We compare for HFB (solid line) with its SMMC (open squares) values. Below the pairing transition temperature, the HFB values are strongly suppressed compared to the spherical HF solution (dashed line), a known effect of pairing correlations. The SMMC values are further suppressed compared with HFB, in particular in the vicinity of the pairing transition. We observed substantial suppression also above the pairing transition temperature, indicating the persistence of pairing correlations in the SMMC results, even when the mean-field condensate no longer exists.

While the difficiency of the HFB around the phase transition is interesting to see, the magnitude of the error is not large enough to be of concern in calculating level densities.

v.1.4 State densities and level spacing

In Fig. 16 we show the HFB density of Sm (solid line) in comparison with the SMMC state density (open squares). The HFB density was calculated with the canonical saddle-point approximation (4), taking the canonical entropy from Eq. (20) and from Eq. (21). The result is practically indistinguishable when the term in Eq. (20) is omitted. We observe good agreement between the HFB and SMMC densities for excitation energies above the pairing transitions MeV. At those energies, the HFB solution coincides with the HF solution, and the only role of the pairing is to reset the origin of the excitation energy scale by the pairing correlation energy. This good agreement may be fortuitous in view of two compensating errors in the HFB: the missing correlation energy in the ground state resets the excitation scale to lower the level density, while the many-body correlations increase the level density. This is seen as an increase in the effective mass of the quasi-particles ga58 (). Beyond that, for attractive interactions, the RPA correlation energy further raises the level density va96 (). Nevertheless, we can take the present agreement as support for a popular model of the level density, namely the back-shifted Fermi gas. That model assumes that pairing correlations only affect the origin of the excitation energy scale.

At the same time, the calculated HFB level density is too small at low excitation energies. As with the calculated canonical entropy coming out negative, the problem relates to the violation of particle-number conservation in HFB.

Lastly we compute the neutron resonance spacing at threshold. For that we need the spin dependence of the level density. We use Eqs. (29) and (27), as was done for the SMMC, except that now is taken to be the HFB density. This is justified since the HFB solution is spherical. The spin cutoff factor is taken from the HFB variance of (see Fig. 15). As discussed previously, these fluctuations are larger in the HFB than in SMMC. However, this differences will affect the level density by less than a factor of two. We find in HFB a neutron resonance spacing at the neutron threshold of 4.1 eV, in very good agreement with the SMMC value of 3.7 eV (see Table III).

Vi Conclusion and outlook

Our benchmarking of the finite-temperature HF and HFB approximations for level densities in heavy nuclei provides a quantitative assessment of the limitations of these mean-field theories, but also justifies their use under fairly broad sets of conditions. We have emphasized the relation between the grand canonical and canonical statistics because the mean-field theories are formulated in the grand canonical ensemble, but the actual level densities and the SMMC theory used for benchmarking are canonical. In the FTHF (which provides an appropriate mean-field approximation for a nucleus with weak pairing correlations), we found a simpler way to approximate the projection from the grand canonical to the canonical ensemble when the particle-number fluctuation is small and the standard saddle-point approximation breaks down. The corresponding formulas, Eqs. (20) and (25), are accurate to much less than a unit of entropy, which is entirely acceptable in view of other sources of error. However, we found no simple way to project the HFB to the canonical ensemble with accurate entropies at temperatures below the pairing transition. We will come back to this later.

A fundamental problem is how to treat broken symmetries. As has been long known, both the pairing and shape transitions are rather smooth in the finite-size nucleus, in contrast to the sharp phase changes that are present in the mean-field theory. Otherwise, the issues that arise in the context of level densities are rather different for the shape and the pairing transitions. The deformed shape of the nucleus Dy is quite robust, persisting to excitation energy of MeV, so the change of shape with temperature can more or less be ignored in the statistical regime accessible with low-energy neutrons. The HF state density of Dy is too low by about an order of magnitude at low excitation energy, but the physics is easy to understand: the HF is describing band heads and not the individual rotational states in each band. This may be seen in the calculated canonical entropy (see, e.g., inset of Fig. 8). The SMMC entropy is units at MeV due to the rotational band contribution. In contrast, the HF entropy is close to zero in that region. We attempted to take rotational band physics into account for level densities at small by treating each HF state as a rotational band head, and neglecting the rotational energies as in Ref. bj74, . This gave a level density at neutron threshold that is times larger than the SMMC benchmark. We conclude that a better treatment of rotational band structure is required to achieve a good accuracy for applications.

The HFB approximation for nuclei with strong pairing correlations is a well-established approximation for ground-state properties. Unfortunately, we found no simple way to project from the grand canonical to the canonical ensemble, due to the particle-number mixing that is inherent in the HFB wave function. Furthermore, the odd-even energy staggering of paired systems requires that the variational principle must take into account the number parity of the ensemble. According to Ref. ba99 (), the effects may persist in a nuclear context up to temperatures approaching the mean-field transition point. Fortunately for many applications, the pairing condensate disappears at rather low excitation energy, and the theory above this energy reduces to the HF approximation. The only effect of pairing correlations at these higher energies is to reset the excitation energy scale by adding the pairing correlation energy, as in the . well-known as the “back-shift” models of level densities ra97 (); eg05 (); ga13 (); po14 ().

In view of the comparative benchmarking of a deformed and a spherical nucleus, it is interesting to revisit the arguments presented by Bjornholm, Bohr and Mottelson for effects of breaking the rotational symmetry on the level density bj74 (). First of all, there is no visible effect when comparing the experimental level densities or the benchmark level densities of the two nuclei we examined, Dy and Sm. In fact, we followed the prescription described in Eqs. (8-10) of Ref. bj74, for extracting the level densities from band head densities, and obtained a level density that is too high. Also, the benchmark entropies of both nuclei are very similar for in the range MeV. We conclude that, except for the very lowest excitation energies, the deformation is much less important than commonly assumed.

In a theory of the statistical properties of nuclei, an important source of error is calculating the baseline for the excitation energy of the thermal ensemble. Any correlations in the ground state beyond mean-field theory will lower the base ground-state energy and thus tend to increase the excitation energy of a thermal ensemble. On the other hand, if similar correlations are also present at the excitation energies of interest, their effect will partially cancel the shift of the baseline. The more problematic correlations are those that are present in the ground state but are suppressed in the thermal ensemble at the excitation energies of interest. A possible explanation of the too-high level density calculated for Dy might be that the rotational energy is in this category. On the other hand, we saw that the missing correlation energy of Sm is similar to that of Dy (see Table II). It cannot change very much in Sm and still keep good agreement with the SMMC.

We conclude with some remarks on the possibilities of improving and extending the mean-field treatments discussed here, short of doing the full SMMC sampling.

The mean-field theory is the first approximation in the systematic many-body perturbation theory. The next approximation adds the second-order corrections to the energy or the grand potential . In the infinite Fermi gas, these corrections increase the level density irrespective of the sign of the interaction; it would be interesting to determine if this is the case for the non-collective correlations within the shell model Hamiltonian spaces in use.

The next systematic correction to mean-field theory is the random-phase approximation (RPA). It provides a powerful method to treat correlation energies, even in the presence of degeneracies that are associated with the broken symmetries. With some exceptions va83 (); qu10 (), the RPA has mostly been applied in the framework of the static path approximation (SPA) ke81 (). In early studies it was found that the ground-state energy obtained in the SPA was not accurate enough to be useful for setting the excitation energy scale for the level densities. However, later model studies that included the RPA corrections to pairing have been quite successful pu91 (); ro98 () and the method has been applied to physical systems ka06 (). We note also that the SPA+RPA with the inclusion of number-parity projection describes well thermodynamic properties of superconducting nanoscale metallic grains ne13 (). In the nuclear context, up to now there have not been any systematic study of the accuracy of the SPA+RPA in the framework of the shell model Hamiltonian such as the one used in the SMMC; we plan to examine it in the future.

There is also an interesting recent study within the HF-BCS theory using the combinatorial method rather than the grand canonical ensemble  uh13 (). This method combines the IPM of the the HF mean field together with the pairing energy from BCS with specific blocked orbitals. This requires a large number of BCS calculations, but it was still possible the carry out a systematic survey of the level densities of deformed nuclei that extends up to the actinides. We found that the IPM based on a deformed HF ground state was accurate to within one unit of entropy up to excitation energy of MeV in Dy.

A final issue is the lack of a systematic theory simpler than the SMMC that applies equally well to deformed and paired spherical nuclei. Ref. uh13, limits itself to deformed nuclei only; the global studies reported in Ref. hi06, and successor papers uses different formulas for spherical and deformed nuclei. In principle, the SPA could be a basis for more systematic theory. However, that would require keeping explicit integrations over the two pairing fields (for protons and neutrons) and the two intrinsic quadrupole fields. Furthermore, the problem of setting the excitation energy baseline remains an obstacle to using the SPA+RPA as a global theory.

We hope that the present availability of realistic Hamiltonians and accurate SMMC calculations of their thermal and statistical properties will provide guidance to a renewed search for better methodologies for approximate theories that are less computationally intensive.

Vii Acknowledgments

This work was supported in part by the U.S. DOE grant Nos. DE-FG02-91ER40608 and DE-FG02-00ER411132, and by Grant-in-Aid for Scientific Research (C) No. 25400245 by the JSPS, Japan. The research presented here used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. It also used resources provided by the facilities of the Yale University Faculty of Arts and Sciences High Performance Computing Center.

Appendix I: accuracy of the saddle-point approximation

Here we use a simple model to assess the accuracy of the saddle-point approximations for the state density in Sec. II.1 and the approximate particle-number projections in Secs. II.3 and II.4 for the entropy. The model describes independent fermions populating equidistant energy levels. The only parameter in the model (in the limit when the number of single-particle levels is large) is the spacing of the single-particle states . The Hamiltonian of this system is

 H=δΩ−1∑i=0(i+1/2)a†iai (49)

where is the single-particle level spacing and is the total number of single-particle levels.

We first show that the factorization (26) of the grand canonical partition function is essentially exact in this model when , the temperature in units of the single-particle mean-level spacing, is much smaller than both the number of particles and . The key observation is that under these conditions the canonical partition function , calculated with respect to the ground-state energy of the particles (i.e., in terms of excitation energy) is independent of . Changing shifts the Fermi energy but since the single-particle spectrum is invariant under such a shift, the particle-hole excitations remain unchanged. A typical particle-hole excitation energy is of order so the number of excited particles is typically smaller than (under the above conditions).

The canonical partition function of a system with ground-state energy is given by where is the partition calculated using the excitation energies. The particle ground-state energy of the Hamiltonian (49) is given by , and thus the particle canonical partition is

 Zc(β,N)=e−12βN2δZc(β). (50)

We expand the grand canonical partition function for a value that gives an average number of particles

 α0=βN0δ, (51)

and use Eq. (50) to find

 Zgc(β,α0)≈∑Neα0N−12βN2δZc(β). (52)

The quasi-equality “” is a reminder that the formula is valid in for . Using (51) we have

 Zgc(β,α0)≈∑Ne−βδ(N−N0)22e12βN20δZc(β). (53)

With the help of Eq. (50) for , we can rewrite the last relation in the form

 Zgc(β,α0)e−α0N0=Zc(β,N0)(∑Ne−βδ(