Benchmarking meanfield approximations to level densities
Abstract
We assess the accuracy of finitetemperature meanfield 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 finitetemperature HartreeFock and HartreeFockBogoliubov 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 wellcontrolled statistical errors. The main weak points in the meanfield treatments are seen to be: (i) the extraction of numberprojected 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 saddlepoint approximation to extract the numberprojected densities is not a significant source of error compared to other errors inherent to the meanfield theory. We also present an alternative formulation of the saddlepoint approximation that makes direct use of an approximate particlenumber projection and avoids computing the usual threedimensional Jacobian of the saddlepoint integration. We find that the pairing condensate is less amenable to approximate particlenumber projection methods due to the explicit violation of particlenumber conservation in the pairing condensate. Nevertheless, the HartreeFockBogoliubov 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 “backshift” approximation, treating pairing as only affecting the excitation energy scale. When the ground state is strongly deformed, the HartreeFock 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 lowenergy 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 meanfield theory whose accuracy is not well understood. Our goal here is to assess these meanfield approximations by taking the Hamiltonian as known and testing them against a theory that is accurate up to wellcontrolled statistical errors.
Most treatments of the level density of heavy nuclei start from a meanfield theory in the form of the finitetemperature HartreeFock (HF) or the HartreeFockBogoliubov (HFB) approximation. The finitetemperature 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 auxiliaryfield 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 finitetemperature 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 longrange correlations in a realistic way. These include deformations, pairing gaps and lowenergy collective excitations. It is thus wellsuited to provide a benchmark for testing the validity of the meanfield 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 meanfield pairing condensate. In contrast, the formulas for level densities in HF and HFB depend on the character of the meanfield 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 welldeformed nucleus with a weak pairing condensate; here the HF is an appropriate meanfield 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 meanfield 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
(1) 
as
(2) 
This allows one to calculate the partition function from . Alternatively, can be calculated directly from by integrating the relation
(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 saddlepoint approximation
(4)  
where in the above expression is determined as a function of from the saddlepoint condition
(5) 
The entropy of a system whose Hamiltonian is defined in a finitedimensional model space satisfies a sum rule obtained from (2) in the limit
(6) 
For a finitedimensional model space, the entropy at both end points and is finite and can be determined analytically. In particular, for eveneven nuclei, the entropy at zero temperature must be zero.
In the configurationinteraction (CI) shell model, we have a certain number of nucleons in singleparticle model spaces of dimensions giving a canonical entropy at of
(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 finitetemperature HF (FTHF) and finitetemperature 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
(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 singleparticle shellmodel 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
(9) 
where are occupation factors for .
The grand canonical partition satisfies
(10) 
The Legendre transform of with respect to is a function of and defined by and satisfies
(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 threedimensional inverse Laplace transform
(12) 
Normally one carries out the integration over all three variables in the threedimensional (3D) saddlepoint approximation, resulting in the formula for the state density BM ()
(13) 
where the values of are determined from by the saddlepoint conditions
(14) 
ii.3 Ensemble reduction
To compare the meanfield entropies to the canonical SMMC entropies we need to reduce the grand canonical ensemble of the meanfield 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 variationafterprojection 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 3D saddlepoint integration Eq. (II.2) into two steps, integrating first over the chemical variables . This yields the following expression for the integrand of the integration:
(15) 
where
(16) 
and are determined by the 2D saddlepoint conditions (). Comparing with the integrand in Eq. (4), we can identify the approximate canonical partition function as
(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
(18) 
The expression we find for the state density is equivalent to Eq. (13) of the 3D saddlepoint 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 saddlepoint condition that determines in terms of the energy becomes , where is an approximate canonical energy given by
(19) 
The level density is then given by the canonical form (4), where the canonical entropy is
(20) 
A simple model is presented in Appendix I, showing that the saddlepoint shift in Eq. (19) improves the accuracy of the calculated and .
ii.4 Discrete Gaussian model and the particlenumber fluctuation
Another source of error may arise from the Gaussian approximation in the 2D saddlepoint 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 saddlepoint 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 ,
(21) 
Expression (21) for reduces to the saddlepoint 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
(22) 
where is the neutronnumber fluctuation in the grand canonical ensemble. Carrying out the Gaussian saddlepoint integration, in the 2D saddlepoint approximation (15) becomes
(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
(24) 
in the limit that .
The finitetemperature meanfield approximation also provides a manyparticle density matrix, so that the particlenumber 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 finitetemperature meanfield approximations we use here, the offdiagonal particlenumber correlations vanish for and factorizes into two separate factors for protons and neutrons
(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
(26) 
Relation (26) describes the factorization of the grand canonical partition function into a canonical partition function and particlenumber fluctuation partition functions. It is exact in the simple model presented in Appendix I.
In summary, the saddlepoint 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 particlenumber fluctuation is small. We will demonstrate the improvement to the saddlepoint formula in this limit in the case of Dy (Sec. IV.2) where pairing correlations are weak.
ii.5 Spinparity projected level density
The ultimate goal is to calculate the spinparity 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 spindependent 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 spinprojection 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 spinparity 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
(27) 
Here a preexponential factor of , arising from the threedimensional 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
(28) 
The normalization condition of in (27) is . Assuming equal positive and negativeparity level densities (usually justified at the neutron binding energy), we have
(29) 
Iii SMMC results
The SMMC method is formulated in the framework of a CI spherical shellmodel Hamiltonian. The CI Hamiltonians shown to be amenable to MonteCarlo sampling contain one plus twobody operators, with the twobody part restricted to interactions that have a “good sign” in the grand canonical formulation sign (). In finite nuclei, the method is implemented with particlenumber 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.
Sm n  16  34.38  36.55  
Sm p  12  22.44  24.43  
Dy n  26  41.95  44.25  
Dy p  16  24.86  26.92 
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 onebody 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.
nucleus  SMMC  HF  HFB  correlation energy  

HF/HFB  missing  
Sm  0  119.15  119.0  119.0  
230.69  232.51  1.82  3.14  
Dy  0  238.35  238.12  238.12  
371.78  371.91  3.48 
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 spindependent 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 lowlying 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
(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 eveneven 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 groundstate 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 510 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 spinparity 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.
Nucleus  (MeV)  (eV)  


SMMC  HF  HFB  Exp.  
Sm  8.1  4.1  5.7  
Dy  8.2  0.5  2.4 
Iv The finitetemperature HF approximation
We first consider the FTHF approximation. It is derived by minimizing the grand potential in terms of a trial uncorrelated manyparticle density matrix. Such a density is uniquely characterized by its onebody density matrix , where label singleparticle 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 onebody density satisfies the selfconsistent equation
(31) 
where is the onebody HF Hamiltonian, expressed respectively in terms of the one and twobody matrices and of the configurationinteraction shell model Hamiltonian. The value of at the HF minimum is given by
(32) 
where is the HF approximation to the grand canonical partition function. The thermal HF energy is calculated as
(33) 
the entropy is given by
(34)  
and the average number of particles is computed as
(35) 
The occupation probabilities are the usual FermiDirac occupations where are the singleparticle 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 singleparticle 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
(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 saddlepoint approximation of Sect. II.3 or in the discrete Gaussian model of Sect. II.4. However, the connection to particlenumber fluctuations is more tenuous because the second derivative expressions for in terms of particlenumber 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 particlenumber fluctuations. We make such a comparison for Dy in the next section.
iv.2 Application to Dy
Here we discuss the strongly deformed nucleus Dy. The pairing in this nucleus is weak, so that the FTHF is the appropriate meanfield 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 particlenumber fluctuations in FTHF are given by
(37) 
and the offdiagonal 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 sphericaltodeformed 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 particlenumber 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 hightemperature 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 groundstate 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 dasheddouble 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 groundstate 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 meanfield theory overemphasizes phase transitions in finite systems is wellknown 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 particlenumber 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 saddlepoint canonical entropy calculated from Eqs. (15) and (16) increases at large values of (dotted line in the inset), indicating the breakdown of the saddlepoint approximation to the particlenumber 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 saddlepoint 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 groundstate rotational band, most visible at MeV. We can estimate this contribution as follows. The moment of inertia of the groundstate band of Dy was determined to be MeV by fitting the lowtemperature 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 groundstate rotational band. We find
(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 (dasheddouble dotted line). This entropy approaches a finite nonzero 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
(39) 
where is the matrix representing in the singleparticle space. When the HF equilibrium ensemble is axially symmetric around the axis, is invariant under rotations around the axis and . Assuming timereversal 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 saddlepoint 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
(40) 
where
(41) 
and
(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
(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 Independentparticle model
A common approximation in the calculation of state densities is to take the HF groundstate mean field, and assume the excited states can be calculated in the independentparticle model (IPM) with singleparticle energies given by that potential field. For an axially deformed nucleus such as Dy, the singleparticle HF levels come in doubly degenerate timereversed pairs and for an even number of particles the ground state is nondegenerate, so that the entropy is zero. It is rather easy to carry out the exact paritclenumber 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 finitetemperature HFB approximation
The HFB is the preferred meanfield 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 particlenumber 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
(44) 
The HFB entropy is given by
(45) 
where
(46) 
are the quasiparticle occupations expressed in terms of the quasiparticle 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 particlenumber 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
(47) 
v.1 Application to Sm
The meanfield ground state of Sm is spherical and has a substantial pairing condensate. Thus FTHFB is the proper meanfield theory for this nucleus. The pairing correlation energy, defined as the difference between the HF and HFB groundstate 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 groundstate 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 particlenumber 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 particlenumber fluctuation. We have also calculated from Eq. (21) using the full matrix suggested by the 3D saddlepoint 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 particlenumber fluctuations in FTHFB remain relatively large even at low temperatures, similar results for are found in the saddlepoint 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 particlenumber 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 particlenumber 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
(48)  
The additional contribution is negative, leading to a reduction in the meansquare 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 meanfield 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 saddlepoint 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 manybody correlations increase the level density. This is seen as an increase in the effective mass of the quasiparticles 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 backshifted 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 particlenumber 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 finitetemperature HF and HFB approximations for level densities in heavy nuclei provides a quantitative assessment of the limitations of these meanfield 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 meanfield 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 meanfield 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 particlenumber fluctuation is small and the standard saddlepoint 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 finitesize nucleus, in contrast to the sharp phase changes that are present in the meanfield 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 lowenergy 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 wellestablished approximation for groundstate properties. Unfortunately, we found no simple way to project from the grand canonical to the canonical ensemble, due to the particlenumber mixing that is inherent in the HFB wave function. Furthermore, the oddeven 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 meanfield 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 . wellknown as the “backshift” 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. (810) 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 meanfield theory will lower the base groundstate 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 toohigh 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 meanfield treatments discussed here, short of doing the full SMMC sampling.
The meanfield theory is the first approximation in the systematic manybody perturbation theory. The next approximation adds the secondorder 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 noncollective correlations within the shell model Hamiltonian spaces in use.
The next systematic correction to meanfield theory is the randomphase 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 groundstate 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 numberparity 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 HFBCS 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. DEFG0291ER40608 and DEFG0200ER411132, and by GrantinAid 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. DEAC0205CH11231. 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 saddlepoint approximation
Here we use a simple model to assess the accuracy of the saddlepoint approximations for the state density in Sec. II.1 and the approximate particlenumber 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 singleparticle levels is large) is the spacing of the singleparticle states . The Hamiltonian of this system is
(49) 
where is the singleparticle level spacing and is the total number of singleparticle 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 singleparticle meanlevel 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 groundstate energy of the particles (i.e., in terms of excitation energy) is independent of . Changing shifts the Fermi energy but since the singleparticle spectrum is invariant under such a shift, the particlehole excitations remain unchanged. A typical particlehole 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 groundstate energy is given by where is the partition calculated using the excitation energies. The particle groundstate energy of the Hamiltonian (49) is given by , and thus the particle canonical partition is
(50) 
We expand the grand canonical partition function for a value that gives an average number of particles
(51) 
and use Eq. (50) to find
(52) 
The quasiequality “” is a reminder that the formula is valid in for . Using (51) we have
(53) 
With the help of Eq. (50) for , we can rewrite the last relation in the form