Benchmarking mean-field approximations to level densities
aa r X i v : . [ nu c l - t h ] D ec 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 351560University of Washington, Seattle, Washington 98915, USA Department of Physics, Graduate School of Science,Chiba University, Inage, Chiba 263-8522, Japan
We assess the accuracy of finite-temperature mean-field theory using as a standard the Hamilto-nian 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 heavyspherical nucleus with strong pairing correlations. The errors inherent in the finite-temperatureHartree-Fock and Hartree-Fock-Bogoliubov approximations are analyzed by comparing the entropiesof the grand canonical and canonical ensembles, as well as the level density at the neutron reso-nance threshold, with shell model Monte Carlo (SMMC) calculations, which are accurate up towell-controlled statistical errors. The main weak points in the mean-field treatments are seen tobe: (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 pairingcondensate, we confirm that the usual saddle-point approximation to extract the number-projecteddensities is not a significant source of error compared to other errors inherent to the mean-field the-ory. We also present an alternative formulation of the saddle-point approximation that makes directuse of an approximate particle-number projection and avoids computing the usual three-dimensionalJacobian of the saddle-point integration. We find that the pairing condensate is less amenable toapproximate particle-number projection methods due to the explicit violation of particle-numberconservation in the pairing condensate. Nevertheless, the Hartree-Fock-Bogoliubov theory is accu-rate to less than one unit of entropy for
Sm at the neutron threshold energy, which is above thepairing phase transition. This result provides support for the commonly used “back-shift” approx-imation, treating pairing as only affecting the excitation energy scale. When the ground state isstrongly deformed, the Hartree-Fock entropy is significantly lower than the SMMC entropy at lowtemperatures due to the missing contribution of rotational degrees of freedom. However, treatingthe rotational bands in a simple model, we find that the entropy at moderate excitation energies isreproduced to within two units, corresponding to an error in the level density of less than an orderof magnitude. We conclude with a discussion of methods that have been advocated as beyond themean field approximation, and their prospects to ameliorate the issues we have identified.
I. INTRODUCTION
Nuclear level densities are important input in the the-ory of low-energy nuclear reactions. In situations wherethe reactions cannot be studied in the laboratory, theneed for a reliable theory is evident. The calculation oflevel densities in the presence of correlations is a chal-lenging problem. Part of the problem is the complexityand inadequate knowledge of the nuclear Hamiltonian.But even with a known Hamiltonian, it is often neces-sary to make approximations based on mean-field theorywhose accuracy is not well understood. Our goal here isto assess these mean-field approximations by taking theHamiltonian as known and testing them against a theorythat is accurate up to well-controlled statistical errors.Most treatments of the level density of heavy nu-clei start from a mean-field theory in the form of thefinite-temperature Hartree-Fock (HF) or the Hartree-Fock-Bogoliubov (HFB) approximation. The finite-temperature theory is derived from a variational prin-ciple based on the grand potential as a function of anuncorrelated trial density [1, 2]. These approximationshave been widely applied and taken as a starting pointfor more sophisticated theories [3, 4]. It would thus be useful to validate them by comparing with a more accu-rate method. The auxiliary-field Monte Carlo method,known as the shell model Monte Carlo (SMMC) in nu-clear physics [5, 6], fulfills this role. Given the Hamilto-nian and a model space, the only inaccuracy is a control-lable statistical error associated with the Monte Carlosampling. As a finite-temperature method, SMMC isparticularly suitable for the calculation of level densi-ties [7, 8].The SMMC starts from a Hamiltonian that is some-what restricted but which reproduces all the long-rangecorrelations in a realistic way. These include deforma-tions, pairing gaps and low-energy collective excitations.It is thus well-suited to provide a benchmark for testingthe validity of the mean-field treatments of nuclear ther-mal and statistical properties. We note that the SMMCapplies to all nuclei irrespective of whether they are de-formed or spherical and independently of the existenceof a mean-field pairing condensate. In contrast, the for-mulas for level densities in HF and HFB depend on thecharacter of the mean-field solution.In Sec. II we summarize the statistical and thermo-dynamic tools we will use for the analysis and compar-isions. In Sec. III we present in some detail the re-sults of the SSMC calculations for two nuclei, namely
Sm and
Dy. The first is a spherical nucleus havinga strong pairing condensate in HFB. The second is a well-deformed nucleus with a weak pairing condensate; herethe HF is an appropriate mean-field approximation. InSec. IV we present the HF approximation and its resultsfor
Dy, while in Sec. V we discuss the HFB approxi-mation and its results for
Sm. In Sec. VI we summa-rize our findings regarding the accuracy of the HF andHFB approximations with some remarks on prospects forextending the mean-field approach to level densities. Fi-nally, the data files from the SMMC, HF and HFB cal-culations together with the codes to apply them are pro-vided in the Supplementary Material depository accom-panying this article.
II. TOOLS OF STATISTICAL THEORYA. The canonical entropy
A good meeting point for comparing different sta-tistical theories is the canonical entropy function S c ( β, N p , N n ), i.e., the entropy of the canonical ensembleof states having fixed numbers of protons and neutrons N p , N n and at inverse temperature β . In the SMMC, thequantity that is most directly computed is the thermalenergy E c ( β ) of the canonical ensemble. The entropy canthen computed by integrating the thermodynamic rela-tion dS c = β dE c (1)as S c ( β ) = S c (0) − Z E c (0) E c ( β ) β ′ dE c . (2)This allows one to calculate the partition function Z c from Z c = exp( − βE c + S c ). Alternatively, Z c can becalculated directly from E c by integrating the relation d ln Z c = − E c dβ , (3)and the canonical entropy is then calculated from S c =ln Z c + βE c . Finally, the density of states ρ ( E, N p , N n ) atgiven energy E and particle numbers N p , N n is obtainedfrom Z c by an inverse Laplace transform, carried out inthe saddle-point approximation ρ ( E, N p , N n ) = 12 πi Z i ∞− i ∞ dβ ′ e β ′ E Z c ≈ (cid:18) π (cid:12)(cid:12)(cid:12)(cid:12) ∂E∂β (cid:12)(cid:12)(cid:12)(cid:12)(cid:19) − / e S c ( β ) , (4)where β in the above expression is determined as a func-tion of E from the saddle-point condition E = − ∂ ln Z c ∂β = E c ( β ) . (5) The entropy of a system whose Hamiltonian is definedin a finite-dimensional model space satisfies a sum ruleobtained from (2) in the limit β → ∞ Z E ( β =0) E ( β = ∞ ) β dE = S (0) − S ( ∞ ) . (6)For a finite-dimensional model space, the entropy at bothend points β = 0 and β = ∞ is finite and can be deter-mined analytically. In particular, for even-even nuclei,the entropy at zero temperature S c ( ∞ ) must be zero.In the configuration-interaction (CI) shell model, wehave a certain number N p , N n of nucleons in single-particle model spaces of dimensions D p , D n giving acanonical entropy at β = 0 of S c (0) = S p (0) + S n (0) = ln (cid:18) D p N p (cid:19) + ln (cid:18) D n N p (cid:19) . (7)We have found the sum rule in Eq. (6) useful for testingthe computer codes we have employed in this study, andalso for setting end points on the entropy plots we showlater. B. The grand canonical entropy
The finite-temperature HF (FTHF) and finite-temperature HFB (FTHFB) approximations are definedin the framework of the grand canonical ensemble whichdepends on the additional independent variables α i ( i = p, n ), related to the chemical potentials µ i by µ i = α i /β .The grand canonical entropy S gc , when expressed asa function of the energy E gc and the average number ofparticles N i,gc of type i in the grand canonical ensemble,satisfies dS gc = βdE gc − X i = p,n α i dN i,gc , (8)and can be calculated as in Eq. (2) for fixed N p,gc , N n,gc ,i.e., when the grand canonical energy E gc is expressed asa function of β and the given particle numbers N i,gc .The grand canonical entropy at β = 0 is also asimple function of the dimension of the single-particleshell-model space in the CI shell model, since all cor-relations disappear in that limit. Choosing values for α p = βµ p , α n = βµ n (where µ p , µ n are chemical poten-tials) to produce average particle numbers N p , N n , wehave S gc (0 , N p , N n ) = − X i = p,n D i [ f i ln f i + (1 − f i ) ln(1 − f i )] , (9)where f i = N i /D i are occupation factors for i = p, n .The grand canonical partition Z gc = Z gc ( β, α p , α n )satisfies d ln Z gc = − E gc dβ + X i N i,gc dα i . (10)The Legendre transform of ln Z gc with respect to α i isa function of β and N i,gc defined by ln ˜ Z gc = ln Z gc − P i α i N i,gc and satisfies d ln ˜ Z gc = − E gc dβ − X i α i dN i,gc . (11)Thus we can alternatively calculate the grand canonicalentropy by integrating − E gc ( β, N i,gc ) with respect to β at fixed N i,gc to determine ln ˜ Z gc and then find the en-tropy from S gc ( β, N i,gc ) = ln ˜ Z gc + βE gc .The state density ρ ( E, N p , N n ) is related to the par-tition function Z gc ( β, α p , α n ) by the three-dimensionalinverse Laplace transform ρ ( E, N p , N n ) = 1(2 πi ) Z i ∞− i ∞ dβ Z i ∞− i ∞ dα p Z i ∞− i ∞ dα n × e βE − α p N p − α n N n Z gc ( β, α p , α n ) . (12)Normally one carries out the integration over all threevariables in the three-dimensional (3-D) saddle-point ap-proximation, resulting in the formula for the state den-sity [9] ρ ( E, N p , N n ) = 1(2 π ) / (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( E, N p , N n ) ∂ ( β, α p , α n ) (cid:12)(cid:12)(cid:12)(cid:12) − / e S gc , (13)where the values of β, α p , α n are determined from E, N p , N n by the saddle-point conditions E = − ∂ ln Z gc ∂β = E gc ( β, α p , α n ) ,N i = ∂ ln Z gc ∂α i = N i,gc ( β, α p , α n ) . (14)Using Eqs. (14), the Jacobian in (13) can be writtenas the determinant of the matrix of second derivatives ofthe logarithm of the grand canonical partition functionwith respect to β, α p , α n . C. Ensemble reduction
To compare the mean-field entropies to the canonicalSMMC entropies we need to reduce the grand canoni-cal ensemble of the mean-field formalism to the canoni-cal ensemble. For approximate theories such as the HFand HFB, the only consistent way (in the sense of Ap-pendix II) to carry out the reduction is the variation-after-projection method (VAP), but this is difficult andmas far as we know, has never been put into practice forcalculating level densities in heavy nuclei. We will there-fore only consider simpler reduction methods, recognizingthat they cannot be free from ambiguity.A straightforward way to determine a canonical en-tropy is to separate the 3-D saddle-point integrationEq. (12) into two steps, integrating first over the chemicalvariables ( α p , α n ). This yields the following expressionfor the integrand of the β integration: ζ − Z gc ( β, α p , α n ) e βE − P i α i N i , (15) where ζ = 2 π (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( N p , N n ) ∂ ( α p , α n ) (cid:12)(cid:12)(cid:12)(cid:12) / . (16)and α p , α n are determined by the 2-D saddle-point condi-tions N i = ∂ ln Z gc /∂α i ( i = p, n ). Comparing with theintegrand in Eq. (4), we can identify the approximatecanonical partition function asln Z c ≈ ln Z gc − X i α i N i − ln ζ . (17)If we carry out in a second step the β integration ofEq. (15), we obtain an expression of the same form asin Eq. (4) where the approximate canonical entropy isgiven by S c ≈ S gc − ln ζ . (18)The expression we find for the state density ρ is equiv-alent to Eq. (13) of the 3-D saddle-point approxima-tion [10].However, the above result does not take into accountthe variation of the prefactor ζ with respect to β . If weconsider this dependence explicitly when we perform the β integration, the saddle-point condition that determines β in terms of the energy E becomes E = E ζ , where E ζ is an approximate canonical energy given by E ζ = E gc − δE, with δE = − d ln ζdβ . (19)The level density is then given by the canonical form (4),where the canonical entropy is S c ( β, N p , N n ) = S gc − ln ζ + βδE . (20)A simple model is presented in Appendix I, showing thatthe saddle-point shift δE in Eq. (19) improves the accu-racy of the calculated S c ( β ) and ρ ( E ). D. Discrete Gaussian model and theparticle-number fluctuation
Another source of error may arise from the Gaussianapproximation in the 2-D saddle-point integration usedto derive Eq. (16). For example, suppose that the grandcanonical partition function were dominant by a singlenuclide ( N p , N n ) in some range of α p , α n . Then we couldapproximate Z gc ( β, α p , α n ) ≈ Z c ( β, N p , N n ) exp( α p N p + α n N n ). Treating this in the saddle-point Gaussian in-tegration gives ζ = 0 (since N i are independent of α j ),rather than the correct answer of ζ = 1. The problemcan be repaired by recognizing that N p and N n are dis-crete integers, not continuous variables. We calculate thematrix ∂N i /∂α j as before but now we calculate ζ as adiscrete sum over particle numbers N i , ζ = X N ′ i ,N ′ j exp − X i,j ∂N∂α (cid:12)(cid:12)(cid:12)(cid:12) − ij ( N ′ i − N i )( N ′ j − N j ) . (21)Expression (21) for ζ reduces to the saddle-point result(16) in the limit when the r.h.s. of (16) is large, but hasthe advantage that it is always larger than 1 and it ap-proaches 1 when the r.h.s. of (16) goes to 0.To see more physically how the approximation (21)works, consider the case when there is only one type ofparticles, say neutrons, and the Hamiltonian used in theGibbs density operator is independent of β and α i . Therequired derivative is then ∂ ln Z gc ∂α n = ∂N n ∂α n = h (∆ N n ) i , (22)where h (∆ N n ) i = h ˆ N n i − N n is the neutron-numberfluctuation in the grand canonical ensemble. Carryingout the Gaussian saddle-point integration, ζ − in the 2-D saddle-point approximation (15) becomes ζ − n = (cid:0) π h (∆ N n ) i (cid:1) − / . (23) ζ − n is just the ratio of states with particle number N n tothe total number of states in an ensemble in which theparticle number N ′ n is distributed as a discrete Gaussian P N ′ n = ζ − n e − ( N ′ n − N n ) / h (∆ N n ) i (24)in the limit that h (∆ N n ) i >> h ∆ N i ∆ N j i vanish for i = j and ζ factorizes into two sep-arate factors for protons and neutrons ζ = ζ p ζ n , where ζ i = X N ′ i e − ( N ′ i − N i ) / h (∆ N i ) i . (25) ζ i in Eq. (25) can be considered as the partition functionwhich describes the fluctuations in the number of parti-cles of type i . The reduction from the grand canonicalto the canonical partition function is then given by Z gc ( β, α p , α n ) e − P i α i N i ≈ Z c ( β, N p , N n ) ζ p ζ n . (26)Relation (26) describes the factorization of the grandcanonical partition function into a canonical partitionfunction and particle-number fluctuation partition func-tions. It is exact in the simple model presented in Ap-pendix I. In summary, the saddle-point approximation breaksdown when 2 π h (∆ N i ) i ≤
1. However, ζ in the discreteGaussian model (Eqs. (25) and (26) or Eq. (21)) alwayssatisfies ζ ≥ Dy (Sec. IV B) where pairing correlationsare weak.
E. Spin-parity projected level density
The ultimate goal is to calculate the spin-parity pro-jected densities ρ J π ( E ), defined as the number of levelsof given angular momentum J and parity π per unit en-ergy, not counting the 2 J + 1 magnetic degeneracy ofthe levels. The spin-dependent level densities ρ J π ( E )can be calculated through an angular momentum projec-tion. The present paper is mainly focussed on the totalstate density and we will not examine in details spin-projection methods. However, to make at least a tenta-tive comparison of the level density at the neutron reso-nance threshold we will calculate them taking a simplifiedmodel for the spin-parity projection. We follow commonpractice and assume the spin distribution is Gaussian inthe three components of the angular momentum vector ~J [11]. Then the fraction of levels having angular mo-mentum J is given by P J ≈ r π J + σ exp (cid:18) − ( J + ) σ (cid:19) . (27)Here a pre-exponential factor of ( J + ) , arising fromthe three-dimensional volume element of ~J , is reduced tothe first power of ( J + ) by dividing out the (2 J + 1)magnetic degeneracy factor. The parameter σ , known asthe spin cutoff parameter, is estimated from the secondmoment of J z σ = h J z i . (28)The normalization condition of P J in (27) is P J (2 J +1) P J = 1. Assuming equal positive- and negative-paritylevel densities (usually justified at the neutron bindingenergy), we have ρ J π ( E ) ≈ P J ρ ( E ) . (29) III. SMMC RESULTS
The SMMC method is formulated in the framework ofa CI spherical shell-model Hamiltonian. The CI Hamil-tonians shown to be amenable to Monte-Carlo samplingcontain one- plus two-body operators, with the two-bodypart restricted to interactions that have a “good sign”in the grand canonical formulation [12]. In finite nu-clei, the method is implemented with particle-numberprojection [13] for both protons and neutrons, and thecalculated observables are the expectation values in thecanonical density matrix. In particular, we consider herethe nuclei
Sm and
Dy. The nucleus
Sm is anexample of a heavy spherical nucleus whose ground statehas significant correlation energy associated with pairing,while the nucleus
Dy has a deformed ground statewith weak pairing.The parameterization of the Hamiltonian and otheraspects of the SMMC calculation have been publishedelsewhere [14, 15]. For reference, Table I and its captiondescribes the model space employed in the calculations. N i d i S c,i (0) S gc,i (0) Sm n 16 8 . · Sm p 12 5 . · Dy n 26 1 . · Dy p 16 6 . · Dy and
Sm and corresponding canonicaland grand canonical entropies at β = 0. The single-particle basis employed in the CI spherical shell model consistsof the orbitals 0 g / , d / , d / , s / , h / , f / for pro-tons and 0 h / , h / , f / , f / , p / , p / , i / , g / for neutrons. The numbers of the single-particle states (in-cluding their magnetic degeneracy) are D p = 40 and D n = 66. N i are the number of valence particles of type i (where theindex i distinguishes neutrons and protons), and d i is thedimension of the many-particle model space for particles oftype i . The canonical and grand canonical entropies S c,i (0)and S gc,i (0) are calculated from Eqs. (7) and (9), respectively. The canonical energy E c = h ˆ H i N p ,N n and the meansquare angular momentum h ˆ J i N p ,N n at fixed numbersof protons and neutrons are calculated directly in SMMCas a function of β . Table II shows the SMMC energies at β = 0 (i.e., the infinite temperature limit) and at high β extrapolated to infinity. The energy at β = 0 is largelydetermined by the one-body part of the Hamiltonian inthe grand canonical ensemble. There is no contributionfrom the direct component of the interaction but there isa small contribution from the exchange terms.The variation of the excitation energy E with β isshown in Fig. 1 using a logarithmic scale for the energy.The excitation energy of Dy is higher than that of
Sm from β = 0 to β ≈ . β ≈ .
5. The higher excitation energy in
Sm near β = 3 is likely due to the collapse of strong pairing inthat nucleus. Similarly, the higher Dy excitation en-ergy at β ≈ h ˆ J i ensemble averages to calculatethe spin-dependent level densities. These are shown inFig. 2 for Sm and
Dy. The higher values of h ˆ J i of Dy at high β are largely due to its deformationand its low-lying first excited 2 + state at ∼ .
08 MeV. Incontrast, h ˆ J i for Sm decreases dramatically at high nucleus β SMMC HF HFB correlation energyHF/HFB missing
Sm 0 -119.15 -119.0 -119.0 ∞ − . ± .
015 -230.69 -232.51 1.82 3.14
Dy 0 -238.35 -238.12 -238.12 ∞ − . ± .
02 -371.78 -371.91 11 .
41 3.48TABLE II: Limiting values of the energies (MeV) calculatedby SMMC, HF, and HFB for
Sm and
Dy. The HF/HFBcorrelation energy is the energy difference between HF andHFB ground states for
Sm, and the difference betweenspherical and deformed ground states for
Dy. The term“missing” denotes the differences between the HFB energies(with both pairing and deformation) and the SMMC ground-state energies. The SMMC energies at β = ∞ include extrap-olation and statistical sampling errors [14]. E x ( M e V ) β (MeV -1 ) Dy Sm 0 0.3 5 10 15 20 E x ( M e V ) β (MeV -1 ) Sm FIG. 1: Canonical excitation energies E x (on a logarith-mic scale) versus β calculated by the SMMC for Sm (opensquares) and for
Dy (solid circles), with lines drawn toguide the eye. The inset shows the large β values using alinear scale for the excitation energy. The Monte Carlo sta-tistical errors are about 0 . β , as expected for a nucleus with a J = 0 ground stateand a gap of ∼ . J = 2 state.At low β , the remaining enhancement for Dy is due toits larger number of active valence nucleons in the modelspace. The errors shown in Fig. 2 are statistical errorsfrom the Monte Carlo sampling.We next apply Eq. (2) to compute the canonicalSMMC entropy. We start from β = 0 with the initialvalue of the canonical entropy given by Eq. (7) and usethe relation Z E c (0) E c ( β ) β ′ dE c = − βE c ( β ) + Z β E c ( β ′ ) dβ ′ , (30)where E c ( β ) is the canonical thermal energy calculatedin SMMC as the thermal expectation value of the Hamil- < J z > β (MeV -1 ) < J z > β (MeV -1 ) FIG. 2: SMMC values of h J z i = h ~J i / β in the canonical ensembles for Sm (open squares) and
Dy (solid circles). The insetshows h J z i = h ~J i / β values. S β (MeV -1 ) Dy Sm 0 6 12 2 4 6 8 10 S β (MeV -1 ) FIG. 3: SMMC entropies of
Sm (open squares) and
Dy (solid circles) for β < − . In this range, theMonte Carlo errors are the size of the symbols or smaller.The inset shows the entropies for the larger β values. tonian. The results are shown in Fig. 3, with the mainfigure showing the low to intermediate values of β , andthe inset showing the large values of β .The Sm and
Dy entropies are nearly equal for β in the range 1 . − . − . We do not know anyobvious reason why that should be the case. At highervalues of β , shown in the inset, one observes a differencebetween the two nuclei. For Sm at β ≥ − the entropy is essentially zero, as expected at low tem-peratures for even-even nucleus with a pairing gap. Onthe other hand, the entropy of Dy remains a couple ofunits higher than zero up to at least β ∼
10 MeV − . This is due to the low excitations associated with the ground-state rotational band, together with the weak pairing inthis nucleus.The state densities calculated from Eq. (4) are shownin Fig. 4. As expected, the state density is higher for Dy than
Sm at low excitation energy. Interest-ingly, they become much closer at excitation energies inthe range 5-10 MeV. As a check on the quality of the
0 5 10 15 20 ρ ( M e V - ) E x (MeV) FIG. 4: State densities vs. excitation energy E x calcu-lated in SMMC using the saddle-point expression Eq. (4) for Dy (solid circles) and
Sm (open squares). The range istruncated at the lower end because the saddle-point approx-imation breaks down when there are only a few levels in arange β − around E . CI Hamiltonian, we compare the calculated level den-sities with the experimental s -wave resonance spacings D = [ P J ρ J π ( E )] − , measured at an energy E that cor-responds to the neutron separation energy. These arecalculated from the spin-parity dependent level density(29) and (27), taking the spin cutoff parameter σ fromEq. (28). The results are shown in Table III. The goodagreement gives us confidence that the Hamiltonian is re-alistic enough to provide useful tests of the HF and HFBapproximations. Nucleus E x (MeV) J π D (eV)SMMC HF HFB Exp. Sm 8.1 (3 − , − ) 3 . ± . Dy 8.2 (2 + , + ) 2 . ± . s -wave resonance spacings D at the excitation en-ergy E x that corresponds to the neutron binding energy. J π are the values of spin and parity of the relevant neutron reso-nance levels. The SMMC, HF and HFB results are comparedwith the experimental values from Ref. 16. The HF spacingis calculated using the model of Sec. IV B 5 to estimate thecontribution of rotational bands. IV. THE FINITE-TEMPERATURE HFAPPROXIMATION
We first consider the FTHF approximation. It is de-rived by minimizing the grand potential Ω in terms of atrial uncorrelated many-particle density matrix. Such adensity is uniquely characterized by its one-body densitymatrix ̺ kl , where k, l label single-particle orbitals in themodel space.For simplicity, we consider only one type of particlesand the results are easily generalized to both protons andneutrons. At the FTHF minimum, the one-body density ̺ satisfies the self-consistent equation ̺ = 1 e βh ̺ − α + 1 , (31)where h ρ = t + v̺ is the one-body HF Hamiltonian, ex-pressed respectively in terms of the one- and two-bodymatrices t and v of the configuration-interaction shellmodel Hamiltonian. The value of Ω at the HF minimumis given by β Ω HF = − ln Z HF = βE HF − S HF − αN HF , (32)where Z HF is the HF approximation to the grand canon-ical partition function. The thermal HF energy E HF iscalculated as E HF = tr( t̺ ) + 12 tr( ̺v̺ ) , (33)the entropy S HF is given by S HF = − tr( ̺ ln ̺ ) − tr[(1 − ̺ ) ln(1 − ̺ )]= − X λ f k ln f k − X k (1 − f k ) ln(1 − f k ) , (34)and the average number of particles N HF is computed as N HF = tr ̺ . (35)The occupation probabilities f k = [1+ e β ( ǫ k − µ ) ] − are theusual Fermi-Dirac occupations where ǫ k are the single-particle HF energies at temperature T .It is not obvious from Eqs. (32), (33) and (35) thatln Z HF satisfies the thermodynamic derivative relations(14) due to the dependence of the HF single-particle en-ergies and density on α and β . Nevertheless, the contri-butions to the derivatives that originate in the implicitdependences on α and β vanish. In particular, the quan-tities defined in the HF theory satisfy − ∂ ln Z HF ∂β = E HF , ∂ ln Z HF ∂α i = N i,HF . (36)The proof is provided in Appendix II. The validity ofthese relations in FTHF guarantee that the grand canon-ical HF entropy S HF satisfies the relation dS HF = βdE HF at fixed average particle numbers N i,HF = N i ,and thus we can compute the entropy in the same wayas we did in the SMMC. A. Approximate canonical projectors in FTHF
The canonical partition function can be approximatedeither in the saddle-point approximation of Sect. II C orin the discrete Gaussian model of Sect. II D. However, theconnection to particle-number fluctuations is more tenu-ous because the second derivative expressions for ln Z gc in terms of particle-number fluctuations such as Eq. (22)no longer holds for ln Z HF . Nevertheless, it is interestingto compare ζ computed with the full matrix ∂N i /∂α j to that computed from the particle-number fluctuations.We make such a comparison for Dy in the next section.
B. Application to Dy Here we discuss the strongly deformed nucleus
Dy.The pairing in this nucleus is weak, so that the FTHF isthe appropriate mean-field theory.
1. Number partition function ζ We first examine the number partition function ζ ob-tained from the approximations we presented in Sec. II.The diagonal particle-number fluctuations in FTHF aregiven by h (∆ ˆ N i ) i = tr[ ̺ i (1 − ̺ i )] (37)and the off-diagonal ones are zero. The corresponding ζ obtained from Eq. (25) is shown as the dashed line inFig. 5. This is compared to ζ calculated from Eq. (21)(using the matrix ∂N i /∂α j ), shown as the solid line.They differ by less than 10%, except for a tiny regionnear the spherical-to-deformed phase transition. Therethe ζ calculated from the Jacobian of ∂N i /∂α j diverges(when approached from the deformed side). Thus, it ap-pears to be a very good approximation to calculate ζ interm of the individual particle-number fluctuations as inEq. (25). In the next section, we shall see that this is notthe case for the FTHFB theory in the presence of strongpairing correlations. In any case, we will use Eq. (21) inthe results shown below.
2. Thermal excitation energy
The range of FTHF energies as a function of β is shownin the fourth column of Table II. The high-temperaturelimit is very close to the SMMC value, since all correla-tion energies disappear in that limit. At the other limitof large β , at or near the ground state, the SMMC en-ergy is about 3.5 MeV lower than its HF value. The HFground-state energy has the correlations associated withthe static deformation, but is missing the rotational en-ergy and other correlation effects. It can also be seen l n ζ β l n ζ β FIG. 5: ln ζ vs. β in FTHF for Dy. Solid line: ln ζ calculated from Eq. (21). Dashed line: ln ζ calculated fromEq. (25). E x ( M e V ) β (MeV -1 ) E x ( M e V ) β (MeV -1 ) FIG. 6: The HF excitation energy of
Dy vs. β . The grandcanonical HF energy (dashed line) is compared with the ap-proximate canonical energy E ζ (solid line) from Eq. (19). Thelatter omits the region near the shape transition point, where ζ becomes singular (see solid line in Fig. 5). We also show theenergy for the spherical HF solution (dashed-double dottedline) with respect to the deformed ground-state energy. Thesolid circles are the SMMC excitation energies from Fig. 1.Inset: expanded energy scale for higher β values. from the fifth column of Table II that the HFB approx-imation hardly lowers the energy. The excitation energy E x ( β ) is shown in Fig. 6, comparing its HF value (solidline) with the SMMC thermal energy (solid circles) fromFig. 1. The HF density is spherical for β < . β is shown inFig. 6 by the dashed-double dotted line. We observe that the HF energy has a cusp at the onset of deforma-tion. Extrapolating the energy of the spherical solutionto large β , we find that the deformed ground-state so-lution is lower in energy than the spherical solution by11.4 MeV. The SMMC energy, shown by the solid cir-cles, is remarkably close to the HF energy in the region0 < β < − except in the vicinity of the shapephase transition. The cusp in the HF energy functiondisappears in the SMMC energy, leaving no trace of ashape phase transition. The fact that mean-field the-ory over-emphasizes phase transitions in finite systems iswell-known in nuclear theory; see, e.g., in Refs. 17–20.The inset in Fig. 6 shows E ( β ) at higher values of β us-ing a finer energy scale. We observe that at large β thegrand canonical HF energy (dashed line) overestimatesthe excitation energy but that the approximate canoni-cal energy E ζ of Eq. (19) (solid line) is in overall goodagreement with the SMMC excitation energy.
3. Entropies
The grand canonical HF entropy S HF is shown inFig. 7 as a function of β (dashed line) and comparedwith the approximate canonical entropy (20) with ζ fromthe discrete Gaussian formula Eq. (21) (solid line). At β = 0, the grand canonical HF entropy is larger than thecanonical entropy due to particle-number fluctuations.The entropies at large values of β are shown in the in-set. The grand canonical HF entropy vanishes in thelimit β → ∞ , as expected. However, the saddle-pointcanonical entropy calculated from Eqs. (15) and (16)increases at large values of β (dotted line in the inset),indicating the breakdown of the saddle-point approxima-tion to the particle-number projection. In contrast, thediscrete Gaussian treatment (21) gives an entropy thatapproaches zero at high β , thus satisfying the sum ruleEq. (6). For moderate and small values of β , the entropy(20) of the discrete Gaussian model (21) essentially coin-cides with the saddle-point canonical entropy. As notedearlier, the SMMC entropy remains nonzero in the range8 < β <
15. We examine this further in the next para-graph.To compare the projected HF and SMMC entropiesin more detail, we have replotted them in Fig. 8 withsome additional information. Both curves start at thesame value at β = 0 because the model spaces are iden-tical. In the limit of large β , the SMMC entropy doesnot approach zero as fast as the canonical HF entropy.This is because the SMMC entropy includes a contribu-tion from the ground-state rotational band, most visibleat β > − . We can estimate this contribution asfollows. The moment of inertia I gs of the ground-stateband of Dy was determined to be I gs / ~ = 35 . ± . − by fitting the low-temperature SMMC valuesof h ~J i to h ~J i = 2( I gs / ~ ) T [14]. For β <
20 MeV − , T > ~ / I gs , and we can treat the rotational motion clas-sically. The classical partition function of the rotational S β (MeV -1 ) S β (MeV -1 ) FIG. 7: Entropy of
Dy in the FTHF approximation, com-paring the grand canonical entropy (dashed line) with thecanonical entropy defined in Eqs. (20) and (19) (solid line).The dashed-dotted line is the approximate canonical entropy(18) in the 3-D saddle-point approximation, i.e., without thecorrection term in Eq. (19). The inset shows the large β value.The calculations use the discrete Gaussian model formula (21)for ζ . The dotted line in the inset uses the saddle-point ex-pression Eq. (16) for β > − . S β (MeV -1 ) S β (MeV -1 ) FIG. 8: Approximate canonical HF entropy defined byEqs. (20) and (19) for
Dy (solid line) is compared with theSMMC entropy (solid circles). The dashed-double dotted lineis the canonical entropy of the spherical HF solution in thesame approximation. The inset shows the entropies at largevalues of β . The dotted line in the inset is the ground-staterotational band contribution (38). band J = 0 , , . . . is given by Z rot = I gs T / ~ . We canthen calculate its entropy from S rot = − ∂F rot /∂T , where F rot = − T ln Z rot is the free energy of the ground-state rotational band. We find S rot = 1 + ln (cid:18) I gs ~ T (cid:19) . (38)This contribution is described by the dotted line in theinset of Fig. 8, and is in good agreement with the SMMCentropy at large β with no adjustable parameters.We also show in Fig. 8 the canonical entropy of thespherical HF solution (dashed-double dotted line). Thisentropy approaches a finite non-zero value in the limit β → ∞ (see inset) because of the large degeneracy of thespherical HF solutions at T = 0. There are 2 valenceprotons in the 0 h / orbitals and 6 valence neutrons inthe 0 h / orbitals, leading to a highly degenerate groundstate with a canonical entropy of ln (cid:2)(cid:0) (cid:1)(cid:0) (cid:1)(cid:3) = 9 . T → p = 12, f p = 2 /
12 andΩ n = 10, f n = 6 /
10 and gives a grand canonical entropyof 13 .
4. Angular momentum fluctuations
In HF, the variance of the angular momentum com-ponents J q ( q = x, y, z ) can be calculated using Wick’stheorem as h (∆ J q ) i = h ˆ J q i − h ˆ J q i = tr[ j q (1 − ̺ ) j q ̺ ] , (39)where j q is the matrix representing ˆ J q in the single-particle space. When the HF equilibrium ensemble isaxially symmetric around the z -axis, ̺ is invariant underrotations around the z axis and h ˆ J x i = h ˆ J y i = 0. As-suming time-reversal invariance, we also have h ˆ J z i = 0.It follows that the variances of the angular momentumcomponents are the same as the mean square moments.In Fig. 9, we compare the HF mean square momentsof ˆ J x and ˆ J z (solid lines) with the SMMC moments h ˆ J x i = h ˆ J y i = h ˆ J z i = h ˆ J i / Dy .The HF mean square moments of ˆ J x and ˆ J z coincideabove the shape transition temperature, where the HFsolution is spherical. However, at large values β , theHF mean square moment of ˆ J x is much larger than therespective moment of ˆ J z . Since the deformed intrinsicground state has good K = 0, h ˆ J z i approaches zero inthe limit β → ∞ , while h ˆ J x i remains finite and largein this limit. We also show by dashed line h ˆ J x i for thespherical HF solution.
5. State density and level spacing
In Fig. 10 we show the HF density vs. E x in the saddle-point approximation (4) (solid line) (where the approx-0 < J x , z > β (MeV -1 )
Dy.The solid lines are the HF results which exhibit a kink at theshape transition point. The dashed line describes the spheri-cal HF solution for temperatures where the lowest equilibriumsolution is deformed. These HF moments may be comparedwith the SMMC moments shown by solid circles. The SMMCmoments satisfy h J x,y i = h J z i = h ~J i /
0 10 20 30 40 50 ρ ( M e V - ) E x (MeV)
0 2 4 6 8 ρ ( M e V - ) E x (MeV) FIG. 10: The HF density of
Dy calculated in the saddle-point approximation (4) using Eqs. (19) and (20) for the ap-proximate energy and canonical entropy (solid line) is com-pared with the SMMC state density (solid circles) as a func-tion of excitation energy E x . The gap in excitation energyreflects the discontinuity of the energy at the shape transi-tion as is seen in Fig. 6. The dashed-dotted line is the ap-proximation in which the δE correction is neglected in thesaddle-point energy Eq. (19) and in the approximate canon-ical entropy of Eq. (20). The inset shows an expanded scaleat low excitation energies. imate canonical entropy and energy include the δE cor-rection in Eqs. (20) and (19) respectively), and compareit with the SMMC state density (solid circles). The kink in the HF density at E ≈
31 MeV signifies the shapetransition from a deformed to spherical shape. At lowerexcitation energies, the HF state density underestimatesthe SMMC values; the SMMC density includes a con-tribution from rotational bands that are built on topof intrinsic K states, and are not captured in the HFapproximation. Above the shape transition energy, theequilibrium shape is spherical and no longer supports ro-tational bands. The HF density is then very close to theSMMC density.We can try to repair the HF approximation by recog-nizing that the each of the deformed HF configurationsrepresents a band [21]. The angular momentum J z corre-sponds to the K -quantum number of the band. Assum-ing that K is Gaussian distributed, the K -dependent HFstate density ρ K can be expressed ρ K ( E ) = P K ρ HF ( E ) . (40)where P K = 1(2 πσ ) / exp (cid:18) − K σ (cid:19) . (41)and σ = h J z i . (42)For each positive K there will be a rotational band with J = K, K + 1 , ... . For K = 0 the sequence may skip oddor even J values. For a complete treatment of the bandmodel one next introduces a moment of inertia of theband to calculate the J -dependent level density. How-ever, we do not wish to introduce new parameters thattake us beyond the HF theory so we assume that all thelevels in a band are degenerate. The level density is then ρ J π ( E ) ≈ J X K =0 ρ K ( E ) , (43)treating the K = 0 bands the same as the others. Thefactor of 1 / E = 8 . Dy is reported in Table III. It un-derestimates the SMMC value by a factor of ∼
5. Thisis a substantial disagreement; in our view uncovers a se-rious problem with the HF theory of level densities indeformed nuclei.
6. Independent-particle model
A common approximation in the calculation of statedensities is to take the HF ground-state mean field,and assume the excited states can be calculated in theindependent-particle model (IPM) with single-particleenergies given by that potential field. For an axiallydeformed nucleus such as
Dy, the single-particle HFlevels come in doubly degenerate time-reversed pairs and1for an even number of particles the ground state is non-degenerate, so that the T = 0 entropy is zero. It israther easy to carry out the exact paritcle-number pro-jection in the IPM [22], so we will use it in the compari-son. In Fig. 11, we compare the IPM state density with
0 2 4 6 8 10 12 14 ρ ( M e V - ) E x (MeV) FIG. 11: The particle-projected IPM density of
Dy (dashed line) based on the zero temperature HF single-particle levels is compared with the HF density (solid line fromFig. 10) as a function of excitation energy E x . the FTHF density determined by using Eq. (20) for S c .They agree very well at low excitation energy. At theneutron separation energy, S n ≈ . V. THE FINITE-TEMPERATURE HFBAPPROXIMATION
The HFB is the preferred mean-field approximationfor nuclei exhibiting strong pairing correlations. Like theFTHF, the FTHFB is based on the grand canonical en-semble. However, unlike the FTHF, the simple approxi-mate particle-number projection onto a canonical ensem-ble is not expected to be a good approximation at lowtemperatures. This will be evident as we go through thesteps to calculate the level density of
Sm, starting fromthe HFB energy function E HF B ( β ). To simplify the no-tation, we consider only one type of particles as we didin the FTHF. The HFB thermal energy is expressed in terms of the normal and anomalous densities ̺, κ as E HF B = tr( t̺ ) + 12 tr( ̺v̺ ) + 14 tr( κ † vκ ) . (44)The HFB entropy S HF B ( β ) is given by S HF B = − X k f k ln f k − X k (1 − f k ) ln(1 − f k ) , (45)where f k = 11 + e βE k (46)are the quasi-particle occupations expressed in terms ofthe quasi-particle energies E k .The HFB partition function satisfies relations similarto Eqs. (36), and thus the entropy can be computed fromthe energy function by an integral similar to Eq. (2).All the expressions regarding the level density in Sec. IVcarry over to the HFB approximation, except that theparticle-number variance in Eq. (37) is now calculatedusing all three Wick contraction terms in the expectationvalue h a † k a k a † l a l i . This leads to an additional contribu-tion from the anomalous density κ h (∆ ˆ N ) i = X k ̺ kk − X kl ̺ kl + X kl | κ kl | . (47) A. Application to Sm The mean-field ground state of
Sm is spherical andhas a substantial pairing condensate. Thus FTHFB isthe proper mean-field theory for this nucleus. The pair-ing correlation energy, defined as the difference betweenthe HF and HFB ground-state energies is 1.82 MeV (seeTable II). Also shown in Table II is the correlation energyof the SMMC ground state with respect to the HFB so-lution (i.e., the difference between the HFB and SMMCground-state energies). This correlation energy is labeled“missing” in the table and is about 3 MeV for
Sm.The pairing transitions for protons and neutrons occurat β = 2 . − ( E x = 5 . β = 2 . − ( E x = 3 . ζ , the factor used to convert the grandcanonical quantities to canonical in the HFB. Naively,we may try using the HFB particle-number fluctuationin Eq. (47) in Eq. (25). The resulting ζ is shown as thedashed line in Fig. 12. This approximation is bound tofail at low temperatures because the HFB ground statehas a nonphysical particle-number fluctuation. We havealso calculated ζ from Eq. (21) using the full ∂N i /∂α j matrix suggested by the 3-D saddle-point approximation.This gives a better result for β > − , as may beseen by the solid line in Fig. 12. We will therefore use thismethod for the canonical quantities calculated below.2 l n ζ β l n ζ β FIG. 12: ln ζ vs. β in FTHFB for Sm. Lines are as inFig. 5. E x ( M e V ) β (MeV -1 ) E x ( M e V ) β (MeV -1 ) FIG. 13: Excitation energy of
Sm as a function of inversetemperature β for 0 < β < − , comparing the grandcanonical HFB energy in Eq. (44) (dashed line) with the ap-proximate canonical energy in Eq. (19) (solid line) and theSMMC results (open squares). The inset shows an expandedenergy scale.
1. Excitation energies
The
Sm thermal excitation energy function E x vs. β in FTHFB is compared with the SMMC results in Fig. 13.The HFB energy is shown by the dashed line. The exci-tation energies at β = 0 are nearly equal, differing onlyby the missing correlation energy. The two kinks in theHFB curve, visible in the inset, are associated with theneutron and proton pairing phase transitions. The HFBexcitation energy is higher than that of the SMMC, as tobe expected from the higher limiting entropy of the grandcanonical 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 notethat for β < − (i.e., for temperatures above thepairing transitions) , the absolute HF and HFB ener-gies coincide, so that the HF excitation energy is lowerthan the HFB excitation energy by exactly the amountof pairing correlation energy in the ground state.
2. Entropies
The grand canonical HFB entropy (dashed line) andSMMC entropy (open squares) functions in
Sm areshown in Fig. 14 vs. β . Their absolute values are set byintegrating from the β = 0 point, where their respectivevalues are given in Table I. Both entropy functions ap-proaches zero at large β , confirming the sum rule (6).The neutron and proton pairing transitions are also vis-ible as kinks in the HFB entropy curve [23]. S β (MeV -1 ) -4 0 4 8 12 4 8 12 16 20 S β (MeV -1 ) FIG. 14: Entropy functions of
Sm: the grand canoni-cal HFB entropy (dashed line) and the approximate canoni-cal HFB entropy in Eq. (20) with ζ given by Eq. (21) (solidline) are compared with the SMMC entropy (open squares).The dashed-dotted line is the entropy associated with the 3-D saddle-point point approximation, i.e., omitting the βδE term in Eq. (20). The inset shows an expanded entropy scaleat large β values. We also show in Fig. 14 the approximate canonicalHFB entropy in Eq. (20) where ζ is given by Eq. (21) inthe discrete Gaussian model. Since the particle-numberfluctuations in FTHFB remain relatively large even atlow temperatures, similar results for ζ are found in thesaddle-point approximation Eq. (16).This approximate canonical entropy coincides with theSMMC at low values of β and overestimates the SMMCentropy around β ∼ − , i.e., in the vicinity ofthe proton pairing transition. At larger values of β , forwhich a nonzero pairing condensate exists, the approxi-mate canonical entropy is in overall agreement with the3SMMC entropy up to β ∼ . − but at lower tem-peratures it becomes negative with a value of about -2 at β ∼ − when the system reaches its HFB groundstate. We note that if ζ were to be calculated using theHFB particle-number fluctuations, the large- β entropywould have been even more negative at ∼ − ζ − of having the proper proton and neutron num-bers for Sm at β = 7 MeV − is only ∼ .
17, hence theunphysical negative entropy at low temperatures.To summarize the results of this section, we have notfound a simple acceptable procedure to project the HFBonto a canonical ensemble if we require the correct en-tropy at high temperature and an error of less that oneunit at low temperatures. We will comment further onthis situation in the conclusion.
3. Angular momentum fluctuations
Eq. (39) which describes the angular momentum fluc-tuations in the FTHF approximation, has an additionalcontribution in FTHFB from a contraction in Wick’s the-orem that involves the anomalous density κ h (∆ J q ) i = h J q i − h J q i = tr[ j q (1 − ̺ ) j q ̺ ] − tr[ j q κj ∗ q κ ∗ ] . (48)The additional contribution is negative, leading to a re-duction in the mean-square moment of the angular mo-mentum. This is just what one would expect as an effectof the pairing correlations.In Fig. 15 we show the angular momentum fluctua-tions in Sm as a function of β . The HFB solutionsin Sm is spherical at all temperatures so all directionsare equal and we only need examine one of them. Wecompare h J z i for HFB (solid line) with its SMMC (opensquares) values. Below the pairing transition tempera-ture, the HFB values are strongly suppressed comparedto the spherical HF solution (dashed line), a known ef-fect of pairing correlations. The SMMC values are fur-ther suppressed compared with HFB, in particular in thevicinity of the pairing transition. We observed substan-tial suppression also above the pairing transition temper-ature, indicating the persistence of pairing correlations inthe SMMC results, even when the mean-field condensateno longer exists.While the difficiency of the HFB around the phasetransition is interesting to see, the magnitude of the erroris not large enough to be of concern in calculating leveldensities. < J z > β (MeV -1 ) FIG. 15: Mean square angular momentum h J z i in Sm,comparing the HFB results (solid line) with the SMMC results(open squares). The dashed-double dotted line correspondsto the spherical HF solution.
4. State densities and level spacing
In Fig. 16 we show the HFB density of
Sm (solidline) in comparison with the SMMC state density (opensquares). The HFB density was calculated with thecanonical saddle-point approximation (4), taking thecanonical entropy from Eq. (20) and ζ from Eq. (21).The result is practically indistinguishable when the δE term in Eq. (20) is omitted. We observe good agreementbetween the HFB and SMMC densities for excitation en-ergies above the pairing transitions E x ≥
0 5 10 15 20 25 30 35 40 ρ ( M e V - ) E x (MeV)
0 4 8 ρ ( M e V - ) E x (MeV) FIG. 16: State densities in
Sm. The HFB density (solidline) calculated from Eqs. (4), (20) and (21) is compared withthe SMMC state density (open squares). The dotted line isthe HFB density with ζ calculated from the particle-numberfluctuations, Eq. (25). The inset is an expanded scale at lowexcitation energies. for the SMMC, except that now ρ ( E ) is taken to be theHFB density. This is justified since the HFB solutionis spherical. The spin cutoff factor σ is taken from theHFB variance of J z (see Fig. 15). As discussed previously,these fluctuations are larger in the HFB than in SMMC.However, this differences will affect the level density byless than a factor of two. We find in HFB a neutronresonance spacing at the neutron threshold of 4.1 eV, invery good agreement with the SMMC value of 3.7 eV (seeTable III). VI. CONCLUSION AND OUTLOOK
Our benchmarking of the finite-temperature HF andHFB approximations for level densities in heavy nucleiprovides a quantitative assessment of the limitations ofthese mean-field theories, but also justifies their use un-der fairly broad sets of conditions. We have emphasizedthe relation between the grand canonical and canonicalstatistics because the mean-field theories are formulatedin the grand canonical ensemble, but the actual level den-sities and the SMMC theory used for benchmarking arecanonical. In the FTHF (which provides an appropriatemean-field approximation for a nucleus with weak pairingcorrelations), we found a simpler way to approximate theprojection from the grand canonical to the canonical en-semble when the particle-number fluctuation is small andthe standard saddle-point approximation breaks down.The corresponding formulas, Eqs. (20) and (25), are accu-rate to much less than a unit of entropy, which is entirelyacceptable in view of other sources of error. However, wefound no simple way to project the HFB to the canonical ensemble with accurate entropies at temperatures belowthe pairing transition. We will come back to this later.A fundamental problem is how to treat broken sym-metries. As has been long known, both the pairing andshape transitions are rather smooth in the finite-size nu-cleus, in contrast to the sharp phase changes that arepresent in the mean-field theory. Otherwise, the issuesthat arise in the context of level densities are rather dif-ferent for the shape and the pairing transitions. Thedeformed shape of the nucleus
Dy is quite robust, per-sisting to excitation energy of ∼
20 MeV, so the changeof shape with temperature can more or less be ignoredin the statistical regime accessible with low-energy neu-trons. The HF state density of
Dy is too low by aboutan order of magnitude at low excitation energy, but thephysics is easy to understand: the HF is describing bandheads and not the individual rotational states in eachband. This may be seen in the calculated canonical en-tropy (see, e.g., inset of Fig. 8). The SMMC entropy is ∼ − β ∼ −
10 MeV − due to the rotationalband contribution. In contrast, the HF entropy is closeto zero in that region. We attempted to take rotationalband physics into account for level densities at small J bytreating each HF state as a rotational band head, and ne-glecting the rotational energies as in Ref. 21. This gavea level density at neutron threshold that is ∼ Dy and
Sm. In fact, we followed the prescriptiondescribed in Eqs. (8-10) of Ref. 21 for extracting the leveldensities from band head densities, and obtained a leveldensity that is too high. Also, the benchmark entropies5of both nuclei are very similar for β in the range ∼ . − − . We conclude that, except for the very lowest ex-citation energies, the deformation is much less importantthan commonly assumed.In a theory of the statistical properties of nuclei, animportant source of error is calculating the baseline forthe excitation energy of the thermal ensemble. Any cor-relations in the ground state beyond mean-field theorywill lower the base ground-state energy and thus tend toincrease the excitation energy of a thermal ensemble. Onthe other hand, if similar correlations are also present atthe excitation energies of interest, their effect will par-tially cancel the shift of the baseline. The more problem-atic correlations are those that are present in the groundstate but are suppressed in the thermal ensemble at theexcitation energies of interest. A possible explanation ofthe too-high level density calculated for Dy might bethat the rotational energy is in this category. On theother hand, we saw that the missing correlation energyof
Sm is similar to that of
Dy (see Table II). Itcannot change very much in
Sm and still keep goodagreement with the SMMC.We conclude with some remarks on the possibilitiesof improving and extending the mean-field treatmentsdiscussed here, short of doing the full SMMC sampling.The mean-field theory is the first approximation in thesystematic many-body perturbation theory. The nextapproximation adds the second-order corrections to theenergy or the grand potential Ω. In the infinite Fermigas, these corrections increase the level density irrespec-tive of the sign of the interaction; it would be interestingto determine if this is the case for the non-collective cor-relations within the shell model Hamiltonian spaces inuse.The next systematic correction to mean-field theoryis the random-phase approximation (RPA). It provides apowerful method to treat correlation energies, even in thepresence of degeneracies that are associated with the bro-ken symmetries. With some exceptions [4, 30], the RPAhas mostly been applied in the framework of the staticpath approximation (SPA) [3]. In early studies it wasfound that the ground-state energy obtained in the SPAwas not accurate enough to be useful for setting the exci-tation energy scale for the level densities. However, latermodel studies that included the RPA corrections to pair-ing have been quite successful [18, 31] and the method hasbeen applied to physical systems [19]. We note also thatthe SPA+RPA with the inclusion of number-parity pro-jection describes well thermodynamic properties of super-conducting nanoscale metallic grains [32]. In the nuclearcontext, up to now there have not been any systematicstudy of the accuracy of the SPA+RPA in the frameworkof the shell model Hamiltonian such as the one used inthe 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 thanthe grand canonical ensemble [33]. This method com-bines 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, butit was still possible the carry out a systematic survey ofthe level densities of deformed nuclei that extends up tothe actinides. We found that the IPM based on a de-formed HF ground state was accurate to within one unitof entropy up to excitation energy of ∼
10 MeV in
Dy.A final issue is the lack of a systematic theory simplerthan the SMMC that applies equally well to deformedand paired spherical nuclei. Ref. 33 limits itself to de-formed nuclei only; the global studies reported in Ref. 34and successor papers uses different formulas for sphericaland deformed nuclei. In principle, the SPA could be abasis for more systematic theory. However, that wouldrequire keeping explicit integrations over the two pairingfields (for protons and neutrons) and the two intrinsicquadrupole fields. Furthermore, the problem of settingthe excitation energy baseline remains an obstacle to us-ing the SPA+RPA as a global theory.We hope that the present availability of realistic Hamil-tonians and accurate SMMC calculations of their thermaland statistical properties will provide guidance to a re-newed search for better methodologies for approximatetheories 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 researchpresented here used resources of the National Energy Re-search Scientific Computing Center, which is supportedby the Office of Science of the U.S. Department of En-ergy under Contract No. DE-AC02-05CH11231. It alsoused resources provided by the facilities of the Yale Uni-versity Faculty of Arts and Sciences High PerformanceComputing Center.
Appendix I: accuracy of the saddle-pointapproximation
Here we use a simple model to assess the accuracy ofthe saddle-point approximations for the state density inSec. II A and the approximate particle-number projec-tions in Secs. II C and II D for the entropy. The model de-scribes independent fermions populating equidistant en-ergy levels. The only parameter in the model (in the limitwhen the number of single-particle levels is large) is thespacing of the single-particle states δ . The Hamiltonianof this system is H = δ Ω − X i =0 ( i + 1 / a † i a i (49)6where δ is the single-particle level spacing and Ω is thetotal number of single-particle levels.We first show that the factorization (26) of the grandcanonical partition function is essentially exact in thismodel when T /δ , the temperature in units of the single-particle mean-level spacing, is much smaller than boththe number of particles N and Ω − N . The key observa-tion is that under these conditions the canonical partitionfunction Z c ( β ), calculated with respect to the ground-state energy of the N particles (i.e., in terms of excita-tion energy) is independent of N . Changing N shifts theFermi energy but since the single-particle spectrum is in-variant under such a shift, the particle-hole excitationsremain unchanged. A typical particle-hole excitation en-ergy is of order T so the number of excited particles istypically smaller than N (under the above conditions).The canonical partition function of a system withground-state energy E is given by e − βE Z c ( β ) where Z c ( β ) is the partition calculated using the excitation en-ergies. The N particle ground-state energy of the Hamil-tonian (49) is given by E = N δ/
2, and thus the N particle canonical partition is Z c ( β, N ) = e − βN δ Z c ( β ) . (50)We expand the grand canonical partition function Z gc ( β, α ) for a value α = α that gives an average num-ber of particles N α = βN δ , (51)and use Eq. (50) to find Z gc ( β, α ) ≈ X N e α N − βN δ Z c ( β ) . (52)The quasi-equality “ ≈ ” is a reminder that the formula isvalid in for T /δ ≪ N, Ω − N . Using (51) we have Z gc ( β, α ) ≈ X N e − βδ ( N − N e βN δ Z c ( β ) . (53)With the help of Eq. (50) for N = N , we can rewritethe last relation in the form Z gc ( β, α ) e − α N = Z c ( β, N ) X N e − βδ ( N − N ! (54)This relation describes the factorization (26) of thegrand canonical partition. The quantity in parenthesis isthe partition function ζ of the discrete Gaussian model[see Eq. (25)], provided that the particle-number varianceis ( βδ ) − . Indeed h (∆ N ) i = Ω X m =0 f m (1 − f m ) ≈ Z Ω0 dm e βmδ − α ( e βmδ − α + 1) = 1 βδ , (55) where we have used βN δ ≫ β (Ω − N ) δ ≫ , N ) = (40 , S c ( β ) (obtained by usingexact particle-number projection) is shown in Fig. 17 bythe solid line. It starts at S (0) = ln (cid:0) (cid:1) = 25 .
65 andapproaches zero at large β .We next turn to the grand canonical ensemble, inwhich we fix the chemical potential at each value of β to get the desired particle number in the ensemble aver-age. The grand canonical entropy for the (40 ,
20) modelis shown as the dashed curve in Fig. 17. Here the startingentropy from Eq. (9) is 27.73, larger than the canonicalentropy by 2.08 units. The approximate reduction to thecanonical ensemble is carried out by Eq. (20). The resultis shown as the dashed curve in Fig. 17. It is accurateto within 0.1 of a unit over entire range of β . We alsoshow for comparison the entropy calculated with Eq. (20)without the δE correction (dotted line). It is much lessaccurate, differing from the exact value by more than ahalf-unit for βδ > S βδ FIG. 17: Entropy of the (Ω , N ) = (40 ,
20) model as a functionof β . See text for explanation. The entropies of Fig. 17 are reploted in Fig. 18 as afunction of excitation energy (in units of δ ) E x /δ . Theentropy shown by the dashed line is the approximatecanonical entropy of Eq. (20) that includes the contri-bution from δE ; in this contribution is omitted in theentropy shown as the dotted line. The two are very closeexcept at low excitation energies (see inset). The differ-ence of the approximate canonical entropy from the truecanonical entropy (solid line) can hardly be seen in thefigure.We next turn to the state density itself. The excitationenergy spectrum is nδ ( n = 0 , , , . . . ), and the state7 S E x / δ S E x / δ FIG. 18: Entropy of the (Ω , N ) = (40 ,
20) model as a functionof excitation energy E x /δ . See text for explanation.
0 5 10 15 20 25 30 35 40 ρ E x / δ
0 1 2 3 4 5 ρ E x / δ FIG. 19: State density for the Hamiltonian (49) with (Ω , N ) =(40 , density for N particles is given by ρ N ( E x ) = X n =0 δ ( E x − nδ ) a ( n ) (56)for n < min( N, Ω − N ). Here a ( n ) = 1 , , , , , , , ... is the well-known partition of the integer n [35, 36]. Thesaddle-point approximation for the state density is alsowell-known [36] and it has been shown to be accurateenough for our purposes for n > δ (i.e., thenumbers a ( n )). The ground state and first excited stateare unique, and then there is an increasing density up tohalf the maximal energy. For comparison we also show bysolid line the level density calculated from (4) using the canonical energy and entropy obtained by exact particle-number projection.The state density in the standard saddle-point approx-imation is shown by the dotted line, while the state den-sity of Eq. (4), in which the saddle-point prefactor iscalculated from dE ζ /dβ (see Eq. (19)) gives the dashed-dotted line. This curve is hardly distinguishable fromthe solid line and improves the agreement with the exactresult at low excitation energies. We observe that theimproved saddle-point expression (20) and (19) is accu-rate to better than 10% to energies as low as 2 δ . Sincethe low-lying region of the spectrum would be calculatedby explicit methods anyway, we conclude that the im-proved saddle-point approximation is entirely adequatefor statistical purposes. Appendix II: thermodynamical consistency of theHF and HFB approximations
The key consistency condition of a finite-temperaturetheory in the grand canonical ensemble is the relationbetween S gc and E gc given by the equation analogous toEq. (2) for the canonical ensemble. It can be easily de-rived from the relations (14) satisfied by the first logarith-mic derivatives of Z gc = Tr exp( − β ˆ H + α ˆ N ). However,these derivatives are more subtle in the case of a mean-field Hamiltonian because the effective Hamiltonian inthe density operator depends on the temperature andchemical potentials.Here we prove that similar relations Eqs. (36) are in-deed valid in the FTHF approximation. The proof fol-lows from the fact that the theory is derived from thevariational principle for the grand potential Ω [1, 2]. Wewrite expression for the HF grand potential Ω HF in theform − ln Z HF = β Ω HF = βE HF − S HF − αN HF . (57)where E HF , S HF were defined in Eqs. (33) and (34).The derivative of Eq. (57) with respect to β has a con-tribution E HF from the explicit dependence on β . How-ever, there is in principle also a contribution from theimplicit dependence on β in E HF , S HF and N HF . To seethat they vanish we go back to the many-body uncorre-lated density matrix ˆ D HF that was the trial density ofthe variational principle. Taking that as the fundamentalvariable, the relevant derivative is − ∂ ln Z HF ∂β (cid:12)(cid:12)(cid:12) α = ∂ ( β Ω HF ) ∂β (cid:12)(cid:12)(cid:12) α = E HF + δ ( β Ω HF ) δ ˆ D HF (cid:12)(cid:12)(cid:12) β,α ∂ ˆ D HF ∂β . (58)Since ˆ D HF is a variational solution at fixed β and α , itfollows that δ ( β Ω HF ) δ ˆ D HF (cid:12)(cid:12)(cid:12)(cid:12) β,α = 0 , (59)8and the second contribution on the right-hand side ofEq. (58) vanishes. The thermodynamic consistency of the finite-temperature HFB can be proven in the sameway. [1] A.L. Goodman, Nucl. Phys. A352
30 (1981).[2] K. Tanabe, K. Sugawara-Tanabe and H.J. Mang, Nucl.Phys.
A357
20 (1981).[3] A.K. Kerman and S. Levit, Phys. Rev. C , 1029 (1981);A.K. Kermin, S Levit and T. Troudet, Ann. Phys.