Average ground-state energy of finite Fermi systems

Semiclassical theories like the Thomas-Fermi and Wigner-Kirkwood methods give a good description of the smooth average part of the total energy of a Fermi gas in some external potential when the chemical potential is varied. However, in systems with a fixed number of particles N, these methods overbind the actual average of the quantum energy as N is varied. We describe a theory that accounts for this effect. Numerical illustrations are discussed for fermions trapped in a harmonic oscillator potential and in a hard wall cavity, and for self-consistent calculations of atomic nuclei. In the latter case, the influence of deformations on the average behavior of the energy is also considered.


I. INTRODUCTION
A basic problem in the physics of finite fermion systems such as, e.g., atoms, nuclei, helium clusters, metal clusters, or semiconductor quantum dots, is the determination of the ground-state energy E. A standard decomposition, deeply rooted in the connection of classical and quantum physics, is to write E as the sum of an average energyĒ and a fluctuating part E [1,2,3]: The largest contribution,Ē, is a smooth function of the number N of fermions. The shell correction E has a pure quantal origin and displays, instead, an oscillatory behavior as a function of N . Equation (1) underlies the usefulness of the so-called mass formulae, like the liquid drop model for nuclei or for metal clusters, of which the oldest example is the wellknown Bethe-Von Weizsäcker mass formula for the binding energy of nuclei. The decomposition (1) is also at the basis of semiclassical and statistical techniques that are used to investigate how the properties of global character of fermion systems vary with the particle number N . Such is the case for instance of the celebrated Thomas-Fermi and Wigner-Kirkwood theories [1,2]. These methods often provide deep physical insights which may be otherwise obscured behind a full quantum calculation.
It is recognized, however, that the semiclassical calculations ofĒ(N ) for fermion systems in either external potentials or self-consistent mean fields show systematic deviations with respect to the actual average of the exact quantum results [1,2,4,5,6,7,8]. For example, in spherically symmetric calculations one finds that, as a function of the number N of particles, the difference E(N )−Ē(N ) between the (fluctuating) exact value E(N ) and the (smooth) semiclassical averageĒ(N ) does not oscillate around zero. In general, for fermions in a fixed external potential, semiclassical calculations overbind the true average of the quantum energy. One of our purposes in the present work is to explain the origin of this effect, and to derive an explicit formula that allows to compute the correct average behavior of E(N ) in fermion systems. Related studies are the works of Refs. [9,10] where a particle number conserving shell correction method has been pursued.
Additional contributions to the average part of the ground-state energy come in fact from a careful analysis of the oscillatory term E(N ). Because this fluctuating function is evaluated at discrete values of the chemical potential (that correspond to integer values of the particle number), its average value is generically non-zero and therefore contributes to the average part of E(N ). This phenomenon is related to the different physical descriptions of quantum mechanical systems obtained from different thermodynamic ensembles, the grand canonical and the canonical in the present context. This subtle topic has played, in recent years, a crucial role in understanding the physics of, e.g., persistent currents in mesoscopic metallic rings [11], or in trapped Bose-Einstein condensates [12].
Our results are illustrated with two schematic models. Namely, we study the average of E(N ) for fermions in a harmonic oscillator (HO) potential, via the semiclassical Wigner-Kirkwood (WK) theory [13], and for fermions in a spherical cavity with sharp boundaries, via the Weyl expansion [14]. In the former case, analytical expressions are available. Finally, and in contrast to the previous examples where the confining potential is fixed, we consider the influence of deformations and self-consistency on the average behavior of E(N ), as well as other related topics. We find that for self-consistent potentials with deformation degrees of freedom the behavior of the average energy is qualitatively different.

II. SMOOTH BEHAVIOR: GRAND CANONICAL VERSUS CANONICAL ENSEMBLES
The usual computation of the different terms in Eq. (1) is as follows. The single-particle level density g(ε) = Tr[δ(ε −Ĥ)] of a quantum fermion system can be expressed as [1,2,13] in terms of the phase space Wigner function f ε ( r, p). We have included a factor 2 to account for spin degeneracy. Then, for a set of fermions in a potential well filled up to an energy µ, the number of states (accumulated level density) and the ground-state energy are obtained from g(ε) through Inserting the Wigner-Kirkwood expansion of the Wigner function f ε ( r, p) in powers of in Eq. (2) produces a smooth functionḡ(ε), where the leading order gives rise to the Thomas-Fermi term. This procedure is well documented in the literature [1,2,15,16,17,18]. Inserting the latter series for g(ε) ≈ḡ(ε) into Eq. (3) yields the semiclassical expansions forN (µ) andĒ(µ). Alternatively, for a Fermi gas contained in a hard wall cavity, one inserts in Eqs. (3) the corresponding Weyl expansion [14] of the average single-particle level densityḡ(ε). In both cases, Eqs. (3) produce a series in decreasing powers of µ whose coefficients depend on the shape of the potential. These expressions provide in general accurate descriptions of the average behavior of g(ε), N (µ), and E(µ). For instance, for an isotropic three dimensional HO potential of frequency ω one obtains the well-known WK expressions [1,2,15,16,17,18] The last term in Eq. (4) contains the derivative of the delta function δ(ε). This term and the last term in Eq. (6) stem from the corrections of order 4 toḡ andĒ, respectively. In the HO potential the 4 contribution tō N vanishes. Figure 1 displays the comparison between the exact quantum mechanical quantities and the smooth approximations (5)- (6). The upper panel shows the accumulated level density N (µ) for a set of fermions in a spherical HO potential, as a function of µ/ ω. The quantum result exhibits discontinuities at each major shell (N = 2, 8, 20, 40, 70, 112 in the present case) and is represented by a staircase function which fluctuates around the smooth WK curve provided by Eq. (5). The oscillatory part of N (µ) (dashed curve) contains the fluctuations due to shell effects. They are seen to oscillate around zero, with a vanishing net average, as µ is varied. The lower panel of Fig. 1 displays the ground state energy E(µ)/ ω for the same potential [19]. Again, the smooth WK curve excellently averages the quantum result and the shell energy fluctuates around zero.
The fact that the average behavior of the remaining shell corrections is zero for E(µ) and N (µ) can be explicitly checked. The general semiclassical theory expresses the fluctuating parts E(µ) = E(µ) −Ē(µ) and N (µ) = N (µ) −N (µ) as sums over the classical periodic orbits of the system at energy µ [1,20,21,22]. Each term in E(µ) and N (µ) is an oscillatory function of the chemical potential (through the action of the corresponding orbit), whose average over a chemical potential window is zero. In the particular case of the HO potential the semiclassical approximation turns out to be exact (see e.g. [1,23]), and takes the form where M r ≡ 2πM , µ r ≡ M r µ/ ω, andN (µ) andĒ(µ) are given by Eqs. (5) and (6), respectively. These expressions illustrate explicitly that the average of the fluctuating parts E(µ) and N (µ) over a chemical potential window is zero, E(µ) µ = N (µ) µ = 0. In real physical Fermi gases with a well defined number of particles the various quantities, like masses or manybody level densities, are not studied as a function of the chemical potential µ but rather as a function of the particle number N . For instance, the ground-state energy E(N ) of the system consists in the sum of the singleparticle energies of the N lowest single-particle states (taking into account spin-degeneracy). Thus Eqs. (7) and (8) are related to the grand canonical ensemble. The qualitative behavior of the function E(N ) as a function of N is in general quite different from the behavior of E(µ).
Based on e.g. the Wigner-Kirkwood method, the usual way to calculate the function E(N ) is as follows. Having determined, in that approximation, the energyĒ(µ) and the accumulated level densityN (µ), one first fixes the chemical potential (or Fermi energy) in terms of the particle-number N by inverting the function To be consistent in the notation, we useμ N to denote the value of the chemical potential for a given N determined from the WK (or Weyl) approximation, to stress that, in this approximation, µ is computed by inverting the smooth partN and not the exact counting function N . Finally, one obtains the smooth Wigner-Kirkwood or Weyl termĒ(N ) replacing µ byμ N inĒ(µ), For example, applying this procedure to the isotropic HO potential, from the leading terms of Eqs. (5) and (6) one straightforwardly obtains which is the leading-order Thomas-Fermi result. It shows that in a HO the leading dependence of the average energy per particle, in units of ω, is ∝ N 1/3 . The full Wigner-Kirkwood functionĒ(N )/N computed for the HO potential including up to the fourthorder contributions in is plotted in Fig. 2 (dashed line) as a function of the particle number N , in units of ω. It is compared to the exact quantum result (solid line). To better visualize the quantum oscillations with changing N , we have subtracted the dominant N 1/3 dependence [recall Eq. (11)] from both the quantum and the WK curves. The upper panel of Fig. 2 displays the results for the isotropic HO potential. The lower panel is for a strongly deformed potential and it will be discussed later on. Focusing on the isotropic HO, one sees that, as expected, the general trend of the smooth WK result turns out to be quite correct in comparison with the global particle number dependence of the quantum energies. There is, however, a systematic deviation in the sense that the WK curve does not pass as a function of N through the average of the quantum values. This is clearly seen from the large asymmetry of the shaded regions above and below the WK curve in the upper panel of Fig. 2. One notices that the WK result overbinds with respect to the true average of the quantum values when N is varied in the spherical symmetry. The same situation prevails in other problems of atomic and nuclear physics as well as in self-consistent mean field calculations [1,2,4,5,6,7,8].

III. CONTRIBUTION OF THE OSCILLATORY CORRECTIONS
The previous results show that the functionĒ(N ) does not describe appropriately the average behavior of E(N ). We now discuss the origin of the discrepancy, and the way to correct it.
In systems with a well defined number of particles the chemical potential µ takes discrete values. These values do not occur at random. For instance, for an even number of particles and non-degenerate single-particle states, a standard rule is to locate the chemical potential halfway between the last occupied and the first unoccupied single-particle states. Fixing a particular rule to determine the chemical potential at a given number of particles introduces a bias in the sampling of the values of µ (with respect to a uniform, random distribution of µ). Due to this bias, when the oscillatory part of the energy E(µ) is evaluated over the set of discrete points it produces, generically, a function whose average is different from zero. To compute that average we proceed as follows.
The decomposition of the single-particle level density into a smooth part and a fluctuating part, whereḡ(ε) is the WK (or Weyl) smooth part and g(ε) is given by the sum over periodic orbits mentioned above [1], induces a corresponding decomposition of the integrated density [cf Eq. (3)]: For a given number of particles N , the chemical potential is defined by inversion of the exact accumulated level density As the particle number N increases, it is natural to decompose the chemical potential into smooth and fluctuating parts: The average partμ N satisfies Eq. (9). Assuming that µ N ≪μ N , a Taylor expansion of the smooth part in powers of µ N around µ =μ N in Eq. (12) (no Taylor expansion is allowed for the fluctuating term N , because it is not a regular function) yields, to lowest order, where we denoteḡ Similarly, the energy may be decomposed as In a system with a well-defined number of particles, the smooth partĒ(N ) of the exact function E(N ) = E(µ N ) was defined in (10), The fluctuating part is thus defined as In order to compute E(N ), and in particular its average over some particle-number window ∆N around N , it is convenient to express the energy in terms of the grand potential Ω = − µ N (ε)dε using the thermodynamic relation Recalling the definition of µ N andμ N [Eqs. (9) and (13)], Eq. (17) may be written Decomposing Ω(µ N ) into its average and fluctuating parts, expandingΩ(µ N ) aroundμ N to second order in µ N , and using the thermodynamic relations we get (14), this takes the form Equation (18) connects the fluctuations of the grand potential (grand-canonical ensemble) to those of the energy at a fixed number of particles (canonical ensemble). This connection, to lowest order, has been exploited in recent years to analyze nuclear-mass fluctuations [24,25]. One may be tempted to think that the average of E(N ) over some particle-number window ∆N around is thus lowered with respect toĒ(N ) (due to the minus sign in front of N 2 (µ N ) in Eq. (18)). However this is wrong because, for the same reasons as for E(N ), the average Ω(µ N ) N = 0. A detailed calculation (cf the Appendix) shows that where s 2 N N is the variance of the spacing s N = ε N +1 −ε N between two consecutive single-particle levels around µ N . Taking the average with respect to the discrete points µ N in Eq.(18), using Eq. (20) for the average of Ω(µ N ) and expressing, for convenience, the average N 2 (µ N ) N over the discrete points µ N in terms of the average N 2 (µ) µ over the continuous variable µ around µ N (cf the Appendix), we get The final expression for the average value of the energy in a system conserving the number of particles is, according to Eqs. (19) and (22), It follows from Eq. (23) that, with respect to the WK or Weyl smooth terms, the true average of the energy as a function of N is increased by a term proportional to the variance of the accumulated level density. Equation (23) contains all the relevant information on the average, and allows to understand the numerical results presented above. Before making a quantitative comparison, we first discuss the general aspects involved in that equation.
Equation (20) is demonstrated in the Appendix for a system without degeneracies (intrinsic and/or due to spin). However, it is easy to see that it is also valid in the presence of degeneracies. This is because the thermodynamic quantities we are considering are continuous variables of a given set of external parameters λ. Assume that for some λ = λ 0 degeneracies occur, and that for slightly different values λ = λ 0 all the degeneracies are lifted (for instance, some of the components of λ may be associated to a shape deformation, with λ = λ 0 the spherical case, and another component may be a magnetic field that lifts the spin degeneracy, with λ = λ 0 = 0 no magnetic field). Then for λ = λ 0 Eq. (20) is valid. One can therefore consider the case with degeneracies as the limit λ → λ 0 and, by continuity, Eq. (20) remains valid.
The variance N 2 (µ) µ in Eq. (23) depends on the system under consideration. However, its general properties can easily be determined. In systems where the typical size of the fluctuations is important, then the shift of the true average E(N ) N with respect toĒ(N ) will also be important. On the other hand, in systems with small fluctuations, N 2 (µ) µ /2ḡ will be small, and the termĒ(N ) will give not only a good approximation to E(N ) N , but also to E(N ) as well (since fluctuations are small). In general, the more regular and/or symmetric a system is, the greater the fluctuations are, and the larger the correction N 2 (µ) µ /2ḡ will be. As the regularities or symmetries are broken, the typical size of the fluctuations diminishes, and E(N ) will be well approximated byĒ(N ). This point is illustrated in Fig. 2, where the upper panel shows E(N )/N for the isotropic HO, where large fluctuations are observed and clear deviations of the average with respect toĒ(N ) are found. In contrast, the lower panel shows a strongly deformed HO, with frequencies ω x /ω = 0.460, ω y /ω = 1.111, and ω z /ω = 1.954 (ω x ω y ω z = ω 3 ), where small fluctuations are observed, and as a consequence good agreement between E(N )/N andĒ(N )/N is found. Another manifestation of the same phenomenon is provided in Fig. 3, where E(N )/N is shown for N = 92 fermions (with spin degeneracy) in a triaxially deformed HO potential as a function of the deformation parameter d, where ω x /ω = δ −1/2 /σ 1/3 , ω y /ω = δ 1/2 /σ 1/3 , and ω z /ω = σ 2/3 , with σ = 1 + d √ 3 and δ = 1 + |d| √ 2. We see that for most deformations (mid-shell configurations) the quantum and the smooth WK values practically agree, up to small fluctuations. Large deviations are observed, instead, when sphericity is approached, and for other special deformations, e.g., for d ∼ 0.65 where the frequency ratio ω x : ω y : ω z is close to 1 : 2 : 3 (when the three frequencies ω x , ω y , and ω z are integer ratios, the energy levels of the HO are degenerate and the classical trajectories of the Hamiltonian become closed periodic orbits [1]).
We have made a quantitative check of Eqs. (22) and (23) for the case of a Fermi gas in a spherical cavity. shows a nontrivial dependence with the particle number (which reflects, to a large extent, the supershell structure), that is very well reproduced by theory.
In the case of the spherical HO, it is possible to get easily an analytical expression for E(N ) N . The function N (µ) is given by the second term in the right hand side of Eq. (7). Squaring it, the main contribution to N 2 (µ) comes from terms where both indices of the double sum are equal. Hence to leading order in µ, we get: N 2 (µ) µ is calculated by averaging over the rapidly fluctuating factors, given by the sine terms. This yields Since N 2 (µ) µ is a smooth function, we replace µ byμ. Thus we can use the WK expression Eq. (5) to compute the dependence of the variance with the number of particles. Using moreover Eq. (4) we finally get, to leading order in N , A comparison with the numerical average of E(N ), obtained from an isotropic 3D HO, is presented in Fig. 5. The result shows an excellent agreement; compared to the spherical cavity, a much simpler N dependence is observed, due to the absence of supershells. Based on general properties of the single-particle spectrum, it is possible to estimate the typical size of the variance N 2 (µ) µ and of its N dependence for a large class of systems, and therefore to estimate E(N ) N . The relevant classification relies on the type of classical dynamics associated to the confining potential. The two extreme cases that can be treated explicitly are fully regular and fully chaotic dynamics (the case of mixed dynamics is more subtle). Based on this classification, explicit results for the typical size of N 2 (µ) µ were obtained in Ref. [26].

IV. DEFORMED AND SELF-BOUND SYSTEMS
Up to now we have considered fermion systems confined by an external potential. This may be applicable to e.g. quantum dot systems or magnetically trapped atomic gases where the self-consistent mean field part plays a minor role with respect to the external confining potential. However, many relevant systems are self bound and then the mean field potential is essentially given by the solution of the self-consistent Hartree-Fock equations which are obtained by minimising the energy of the system. The mean field in these situations may turn out to be spherical, but in many cases rotational invariance will be broken and the mean field becomes deformed. We will see that in these cases the results can show interesting differences with regard to the scenario found in the upper panel of Fig. 2 Fig. 5 for the harmonic oscillator, and in Fig. 4 for the hard wall cavity, where the potential was kept spherical. We want to investigate such cases now.

or in
First, to illustrate the situation, we again consider the HO potential. In contrast to the previous section, for each particle number we now minimise the groundstate energy of the quantum solution with respect to deformation, i.e., with respect to free variation of ω x , ω y , and ω z , under the constraint of volume conservation (ω x ω y ω z = ω 3 ). This must be done in carefully checking simultaneously the optimal choice of the occupancies n x , n y , n z . The semiclassical energiesĒ(N ) always have their absolute minimum at sphericity as given by Eq. (6). As a particular example of self-bound system, we consider the case of the atomic nuclei. We mimic the saturation properties of nuclear forces by including the standard particle-number dependence of the HO frequency ω = 41A −1/3 MeV [2] with A = 2N (i.e., A here represents the mass number of a hypothetical uncharged nucleus with N protons and N neutrons).
In Fig. 6 we show the difference δE between the fully minimised quantum energies and the corresponding isotropic semiclassical expressionĒ(N ) obtained from Eqs. (5) and (6). For comparison, we include in the same figure the fluctuating part E(N ) for the spherical HO (same curve as in the upper part of Fig. 5). We observe that with respect to the purely spherical case, the situation changes very much. Now, contrary to the spherical case, practically all values of δE from the minimised quantum solutions are negative, meaning that the minimised quantum energies stay below the semiclassical curveĒ(N ). The minimised quantum energies coincide with the spherical ones only in a small neighborhood around the shell closures, whereas away from the latter the system is axially deformed or even, around the middle of the shells, a slight triaxiality can appear (in the case of axial symmetry, typical deformations show an axis ratio of 2:3).
It seems natural that the deformed quantum energies are more bound than the approximate energies obtained from the semiclassical theory, in spite of the fact that to our knowledge no upper bound theorem like the Rayleigh-Ritz principle exists for the semiclassical approach. We wish to point out that in Fig. 6 for most values of N the system is actually rather well deformed and that, with the exception of a couple of particle numbers around closed shells, the energy differences are in most cases very close to the zero line, i.e., to the semiclassical WK values. This is consistent with the results obtained in the previous section, where it was shown that for deformed systems where degeneracies are lifted the energy E(N ) is expected to be well approximated by the WK theory.
The magnitude of the difference δE of the minimised quantum solutions to the semiclassical values slightly increases with increasing particle number, as the average curve shows in Fig. 6. However, the magnitude of the same quantity δE divided by the particle number decreases as a function of increasing N , and the minimised deformed quantum energies per particle are extremely close to the semiclassical ones. Notice the opposite trend in Fig. 6 of the two average curves, with and without energy minimisation. In contrast to the latter case, for which an explicit formula for the average behavior was developed and successfully checked in the previous section, we do not yet have an equivalent result for a selfbound system.
We are interested in checking whether this simplified harmonic oscillator scenario remains valid in realistic Hartree-Fock-type mean field calculations. In Ref. [8] self-consistent calculations of the ground-state binding energy of atomic nuclei were carried out using the variational Wigner-Kirkwood method [27]. The nuclear interaction was described by the relativistic mean field (RMF) meson exchange model [28]. Quantum calculations for the RMF model are available in the literature. In particular, a mass table of deformed (axially symmetric) quantum calculations for nuclei with an accurately calibrated RMF nuclear interaction is published in Ref. [29]. From this table we took for each value of the mass number A the most bound (in general deformed) isotope and traced E/A as a function of A [30]. The quantum values together with the RMF semiclassical results, computed following Ref. [8], are shown in Fig. 7 for nuclei with an even number of protons and neutrons. Most of the WK energies lie on top of the deformed quantum energies on the scale of the figure. We plot in addition the RMF quantum values constrained to sphericity. The typical arch structure found in Fig. 2 for the spherical HO potential is then recovered. These arches in nuclei take place between the so-called magic numbers, i.e., the proton or neutron numbers where effects analogous to the shell closures of the HO or of the electron shells in atoms occur. The fact that for nuclei above iron E/A raises whereas in Fig. 2 it keeps decreasing is a trivial effect due to the Coulomb repulsion among protons in the atomic nucleus.
In Fig. 8 we display for the self-consistent RMF the difference δE between the quantum energies (that are, as mentioned, minimised with respect to deformation) and the corresponding semiclassical values (that attain their absolute minima at sphericity). For reference, also the values of the energy difference δE obtained from experimental data [31] instead of the quantum results are displayed. The similarity of Fig. 8 with Fig. 6 for the HO is striking. Systems with the largest deformations are again located mostly around mid shells. When the system approaches the spherical shape, δE becomes increasingly negative and displays the downward peaks seen in Fig. 8 on reaching neutron or proton magic numbers.
It is clear that in the self-consistent case as in the schematic case of the HO with optimized shapes considered above, the average of δE as a function of particle number is, at least for the heavier systems, negative. In self-bound systems we again find that in between shells  , as obtained from self-consistent relativistic mean field calculations. The deformed calculations are from Ref. [29]. For each value of the mass number A we plot the most bound isotope according to the tabulation of Ref. [29] (excepting for A = 40 [30]).
the quantum energies are closer to the semiclassical WK values. In between shells the system deforms in search of minimum energy and avoids the large positive shell corrections to the energy that occur if a spherical shape is kept. As the deformation increases, symmetries are broken and the amplitude of the shell corrections diminishes. This is in agreement with the basic ideas underlying Eq. (23), that imply an energy E(N ) which is well approximated by the WK theory.

V. CONCLUDING REMARKS
In this work we have revisited the old problem of the semiclassical approach to finite fermion systems, based either on the Wigner-Kirkwood expansion with corrections, or the Weyl expansion. We have addressed the nature of one of the most elusive features of the theory, namely the systematic overbinding compared to the quantum average for fermions in fixed potentials.
In the first part, we have shown that this discrepancy is due to the fact that these methods do not incorporate appropriately the conservation of the particle number. There is, generically, a contribution to the average ground-state energy that comes from the fluctuating part, or shell contribution. We derived an explicit formula that takes into account that contribution, and have tested it for different fixed confining potentials. In all cases, a positive correction with respect to the semiclassical result is predicted [cf Eq. (23)], whose magnitude depends on the size of the shell effects. When the confining potential has symmetries, the shell corrections are large, and important deviations between the exact quantum and the WK energies are observed, in agreement  Fig. 7 for the self-consistent relativistic nuclear mean field. The result obtained by taking the difference of the experimental data [31] to the calculated WK energies is also shown for the purpose of illustration. The location of the magic neutron (N ) and proton (Z) numbers is indicated.
with our predictions. In contrast, when symmetries are broken shell effects are smaller, and the exact energies (and not only their average part) are better described by the WK theory.
The description of the behavior of self-bound systems is more difficult and subtle, because at each particle number N the energy is minimised, and hence the shape of the potential is a function of N . In this case, a shell correction which is nearly always negative with respect to the spherical WK result is observed for the HO potential. In between shells, when the system deforms and symmetries are broken, the value of the shell correction is smaller, and the energy is well approximated by the WK theory, in agreement with the general considerations that follow from Eq. (23). Interestingly, all these features have been qualitatively confirmed by a more realistic model based on a mean-field self-consistent calculation of the groundstate energy of atomic nuclei. However, the problem of deriving an explicit formula for the average behavior of the ground-state energy of self-bound systems that correctly takes into account the N -dependence with deformation degrees of freedom is still open.

APPENDIX
We follow here Appendix B of Ref. [32] in order to prove Eqs. (20) and (21).
Let us consider a single-particle spectrum ε j , with j = 1, 2, . . . and ε j ≤ ε j+1 . The accumulated level density N (µ) is discontinuous at each energy level. At the discontinuity, we assign to N (µ) the "intermediate" value N (ε N ) = N − 1/2. For N ≫ 1 and ε N < µ < ε N +1 , writing N (µ) =N (µ) + N (µ), making a Taylor expansion of N (µ) around ε N , and using that N (ε N ) = N − 1/2, we may write Evaluating this relation just before µ = ε N +1 , and taking into account the value of the function at ε N +1 , we have where s N = ε N +1 − ε N is the level spacing (we neglect the dependence ofḡ with energy). Taking the discrete average over N on both sides of Eq. (A2), defined as where ∆N is the number of levels in the window around the N th level, and using that N (ε N +1 ) N = N (ε N ) N , we obtain the (trivial) relation s N N = 1/ḡ. On the other hand, defining the continuous average over a window ∆µ of a function that depends on the chemical potential as The last equality follows because N (µ) µ = 0 by definition. From Eq. (A3) we deduce Squaring and computing the discrete average in both sides of Eq. (A2) it is possible to deduce that N (ε N ) N = 0 after using Eq. (A4).