Microscopic-Macroscopic Approach for Binding Energies with Wigner-Kirkwood Method

The semi-classical Wigner-Kirkwood $\hbar$ expansion method is used to calculate shell corrections for spherical and deformed nuclei. The expansion is carried out up to fourth order in $\hbar$. A systematic study of Wigner-Kirkwood averaged energies is presented as a function of the deformation degrees of freedom. The shell corrections, along with the pairing energies obtained by using the Lipkin-Nogami scheme, are used in the microscopic-macroscopic approach to calculate binding energies. The macroscopic part is obtained from a liquid drop formula with six adjustable parameters. Considering a set of 367 spherical nuclei, the liquid drop parameters are adjusted to reproduce the experimental binding energies, which yields a {\it rms} deviation of 630 keV. It is shown that the proposed approach is indeed promising for the prediction of nuclear masses.


I. INTRODUCTION
Production and study of loosely bound exotic nuclei using Radioactive Ion Beam facilities is of current interest [1,2]. These experiments have given rise to a number of interesting and important discoveries in nuclear physics, like neutron and proton halos, thick skins, disappearance of magicity at the conventional numbers and appearance of new magic numbers, etc. Further, advances in detector systems, and in particular, the development of radioactive beam facilities like Spiral, REX-Isolde, FAIR, and the future FRIB may allow to investigate new features of atomic nuclei in a novel manner.
The study of nuclear masses and the systematics thereof is of immense importance in nuclear physics. With the advent of mass spectrometry, it is possible to measure masses of some of the short lived nuclei spanning almost the entire periodic table [3,4]. For example, the ISOL (isotope separator online) based mass analyzer for superheavy atoms (MASHA) [5,6] coming up at JINR-Dubna will be able to directly measure the masses of separated atoms in the range 112 ≤ Z ≤ 120. The limitation on measurements is set by the shortest measurable half-life, T 1/2 ∼ 1.0 s [5]. The JYFLTRAP [7] developed at the University of Jyväskylä, on the other hand, enables to measure masses of stable as well as highly neutron deficient nuclei (for masses up to A = 120) with very high precision (∼50 keV) [7].
On the theoretical front as well, considerable progress has already been achieved in the accurate prediction of the nuclear masses, and it is still being pursued vigorously by a number of groups around the globe. This is of great importance, since an accurate knowledge of the nuclear masses plays a decisive role in a reliable description of processes like the astrophysical r-process (see, for example, [3]). There are primarily two distinct approaches to calculate masses: a) the microscopic nuclear models based on density functional theory like, Skyrme [8,9] and Gogny [10] Hartree-Fock-Bogoliubov or Relativistic Mean Field (RMF) models [11]), b) microscopic-macroscopic (Mic-Mac) models [12,13,14,15] The Mic-Mac models are based on the well-known Strutinsky theorem. According to this, the nuclear binding energy, hence the mass can be written as sum of a smooth part, and an oscillatory part which has its origins in the quantum mechanical shell effects. The latter consists of the shell correction energy and the pairing correlation energy which in the Mic-Mac models are evaluated in an external potential well. The smooth part is normally taken from the liquid drop models of different degrees of sophistication. The largest uncertainties arise in the calculation of shell corrections. The shell correction is calculated by taking the difference between the total quantum mechanical energy of the given nucleus, and the corresponding 'averaged' energy. Usually, the averaging is achieved by the well-established Strutinsky scheme [16,17]. This technique of calculating the averaged energies runs into practical difficulties for finite potentials, since for carrying out the Strutinsky averaging, one requires the discrete single-particle spectrum, with cut-off well above (at least 3 ω 0 , ω 0 being the major shell spacing) the Fermi energy. For a realistic potential, this condition is not met, since continuum may start within ∼ ω 0 of the Fermi energy. Standard practice is to discretise the continuum by diagonalising the Hamiltonian in a basis of optimum size. A number of Mic-Mac calculations with varying degree of success are available in the literature (see, for example, [12,13,14,15]). The Mic-Mac models typically yield better than ∼0.7 MeV rms deviation in the masses. All these models agree reasonably well with each other and with experiment, but deviate widely among themselves in the regions far away from the valley of stability.
The semi-classical Wigner-Kirkwood (WK) approach [18,19,20,21,22,23,24,25], on the other hand, makes no explicit reference to the single-particle spectrum, and achieves an accurate averaging of the given one-body Hamiltonian. Thus, the WK approach is a good alternative to the conventional Strutinsky smoothing scheme. The quantum mechanical energy is calculated by diagonalising the one-body Hamiltonian in the axially symmetric deformed harmonic oscillator basis with 14 shells. The difference between the total quantum mechanical energy and the WK energy in the external potential well yields the value of the shell correction for a given system. In the present work, we propose to carry out a reliable microscopic-macroscopic calculation of the nuclear binding energies (and hence the masses), employing the semi-classical Wigner-Kirkwood (WK) expansion [18,19,20,21,22,23,24,25] for the calculation of shell corrections instead of the Strutinsky scheme. An exploratory study of using the WK method to compute the smooth part of the energy has been reported earlier to test the validity of the Strutinsky scheme, especially near the driplines [27].
It is known that the WK level density (g W K (ε)) with the 2 correction term exhibits a ε −1/2 divergence as ε → 0, for potentials which vanish at large distances as for instance Woods-Saxon potentials (see, for example, Ref. [26]). The Strutinsky level density, on the contrary, exhibits only a prominent peak as ε → 0. It was therefore concluded in Ref. [28] that the divergence of the WK level density as ε → 0 is unphysical, and the Strutinsky smoothed level density should be preferred. It should however be noted that the WK level densities, energy densities, etc., have to be understood in the mathematical sense of distributions and, consequently, only integrated quantities are meaningful. In fact, it has been shown [25] that the integrated quantities such as the accumulated level densities are perfectly well behaved, even for ε → 0.
Pairing correlations are important for open shell nuclei. In the present work, these are taken into account in the approximate particle number projected Lipkin-Nogami scheme [29,30,31]. Odd-even and odd-odd nuclei are treated in an entirely microscopic fashion (odd nucleon blocking method in the uniform filling approximation), allowing an improved determination of odd-even mass differences, see e.g. the discussion in [32]. The majority of nuclei in the nuclear chart are deformed. In particular, it is well known that inclusion of deformation is important for reliable predictions of nuclear masses. Therefore, here we incorporate in all three deformation degrees of freedom (β 2 , β 4 , γ). To our knowledge, no such detailed and extensive calculation based on the WK method is available in the literature.
The paper is organised as follows. We review the WK expansion in Section 2. The choice of the nuclear, spin-orbit, and Coulomb potentials forms the subject matter of Section 3.
Details of the WK calculations are discussed in Section 4. A systematic study of the WK energies for neutrons and protons as a function of the deformation degrees of freedom is presented in Section 5. The shell corrections for the chains of Gd, Dy and Pb isotopes obtained by using our formalism are reported, and are compared with those calculated employing the traditional Strutinsky averaging technique, in Section 6. Section 7 contains a brief discussion on the Lipkin-Nogami pairing scheme. As an illustrative example, the calculation of the binding energies for selected 367 spherical nuclei is presented and discussed in Section 8. Section 9 contains our summary and future outlook. Supplementary material can be found in appendices A and B.

II. SEMI-CLASSICAL WIGNER-KIRKWOOD EXPANSION
Following Ref. [20], we consider a system of N non-interacting fermions at zero temperature. Suppose that these fermions are moving in a given one-body potential including the spin-orbit interaction. To determine the smooth part of the energy of such a system, we start with the quantal partition function for the system: Here,Ĥ is the Hamiltonian of the system, given by: where V ( r) is the one-body central potential andV LS ( r) is the spin-orbit interaction.
In order to average out shell effects, the simplest one could do is replace the partition function in the above expression by the classical partition function. In this work, we shall use the WK expansion up to fourth order. For brevity, we represent the potentials and form factors without mentioning the dependence on the position vector.
Ignoring the spin-orbit interaction, the WK expansion of the partition function, correct up to fourth order is given by [20]: The spin-orbit interaction, in general, can be written as: whereσ is the unit Pauli matrix, κ is the strength of spin-orbit interaction, and f is the spin-orbit form factor. With the inclusion of such spin-orbit interaction, the WK expansion for the full partition function splits up into two parts: Here, Z (4) (β) is given by Eq. (3), and the spin-orbit contribution to the partition function, correct up to fourth order in , reads [20]: where The level density g W K , particle number N, and energy E can be calculated directly from the WK partition function by Laplace inversion: and where λ is the chemical potential, fixed by demanding the right particle number, and L −1 denotes the Laplace inversion. Using the identity and noting that, in order to get inverse Laplace transforms in convergent form, one obtains the level density for each kind of nucleons assuming spin degeneracy: the particle number: and the energy: It should be noted that we have explicitly assumed that all the derivatives of the potential V and the spin-orbit form factor f exist. The expansion defined here is therefore not valid for potentials with sharp surfaces. This automatically puts a restriction on the choice of the Coulomb potential: the conventional uniform distribution approximation for the charge distribution cannot be used in the present case. We shall discuss this point at a greater length in the next section. The integrals in the above expressions are cut off at the turning points, defined via the step function. The chemical potential λ appearing in these equations is determined from Eq. (16), separately for neutrons and protons. Further, it is interesting to note that the spin-orbit contribution to the particle number N as well as to the energy E appears only in the second order in . Secondly, the level density and particle number are calculated only up to the order 2 . It can be shown [20] that for the expansion correct up to fourth order in , it is sufficient to take Z W K up to order 2 in Eq. (11) to find the chemical potential (and hence the particle number), whereas one has to take the full partition function Z (4) W K up to order 4 in Eq. (12) to compute the energy in the WK approach.
The divergent terms appearing in Eq. (17) are treated by differentiation with respect to the chemical potential. Explicitly: In practice, the differentiation with respect to chemical potential is carried out after evaluation of the relevant integrals. Numerically, this approach is found to be stable. Its reliability has been checked explicitly by reproducing the values of fourth-order WK corrections quoted in Ref. [20].
The WK expansion thus defined, converges very rapidly for the harmonic oscillator potential: the second-order expansion itself is enough for most practical purposes. The convergence for the Woods-Saxon potential, however, is slower than that for the harmonic oscillator potential, but it is adequate [33]. For example, for ∼ 126 particles, the Thomas-Fermi energy is typically of the order of 10 3 MeV, the second-order ( 2 ) correction contributes a few 10's of MeVs, and the fourth-order ( 4 ) correction yields a contribution of the order of 1 MeV. This point will be discussed in greater details later. It is also important to note that the WK expansion of the density matrix has a variational character and that a variational theory based on a strict expansion of the of has been established [34].
The WK approach presented here should be distinguished from the extended Thomas-Fermi (ETF) approach. Divergence problems at the classical turning points (see the particle number and energy expressions above) can be eliminated by expressing the kinetic energy density as a functional of the local density. This is achieved by eliminating the chemical potential, the local potential, and the derivatives of the local potential (for further details, see Ref. [35]). It cannot be accomplished in closed form, and has to be done iteratively, leading to a functional series for the kinetic energy density. The resulting model is what is often referred to as the ETF approach. The WK approach as presented here, in this sense, is the starting point for ETF approach (further details of ETF can be found in Refs. [22,23,25,36,37,38]). The conventional ETF approach exhibits somewhat slower convergence properties which has been attributed to a non-optimal sorting out of terms of each given power in [25,35].

A. Form of the Nuclear Potential
The spherically symmetric nuclear mean field is well represented by the Woods-Saxon (WS) form [39], given by: where V 0 is the strength of the potential, R 0 is the half-density radius, and a is the diffuseness parameter. The WS form factor defined here, can be easily generalised to take the deformation effects into account. Note that the distance function l(r) = r − R 0 appearing in Eq. (20) can be interpreted as the minimum distance of a given point to the nuclear surface, defined by r = R 0 . One might thus generalise it to the case of deformed surfaces as well.
Using the standard expansion in terms of spherical harmonics, a general deformed surface may be defined by the relation r = r s , where Here, the Y λ,µ functions are the usual spherical harmonics and the constant C is the volume conservation factor (the volume enclosed by the deformed surface should be equal to the volume enclosed by an equivalent spherical surface of radius R 0 ): The distance function to be used in the WS potential would be the minimum distance of a given point to the nuclear surface defined by r = r s . Such definition has been used quite extensively in the literature, with good success (see, for example, Refs. [40,41,42,43,44]).
However, in the present case, this definition is not convenient, since it should be noted that the calculation of this distance function involves the minimisation of a segment from the given point to the nuclear surface. This in turn implies that each calculation of the distance function (for given r, θ, and φ coordinates: we are assuming a spherical polar coordinate system here) involves the calculation of two surface angles θ s and φ s , and these are implicit functions of r, θ, and φ. See Fig. (8) in Appendix A for details. Since the WK calculations involve differentiation of the WS function, one also needs to differentiate θ s and φ s , which are implicit functions of r, θ, and φ.
Alternatively, the distance function for the deformed Woods-Saxon potential can be written down by demanding that the rate of change of the potential calculated normal to the nuclear surface and evaluated at the nuclear surface should be a constant [45] which, indeed, is the case for the spherical Woods-Saxon form factor. Thus, wheren is the unit vector normal to the surface (r = r s ) and is given bŷ In fact, the above condition (23) is related to the observation that the second derivative of the spherical Woods-Saxon form factor vanishes at the nuclear surface, defined by r = R 0 .
The resulting distance function is given by [46]: where r s is as defined in Eq. (21). The denominator is evaluated at r = r s . Writing the θ and φ derivatives of r s as A and B respectively, we get: with In the present work, we use the distance function as defined in Eq. (25). The WS potential thus reads It is straightforward to check that the Woods-Saxon potential defined with the distance function as given by Eq. (25) satisfies the condition (23). Substituting this Woods-Saxon potential inn · ∇V ( r), we get Here, f ( r) = [1 + exp (l( r)/a)] −1 is the Woods-Saxon form factor. Clearly, at the surface defined by r = r s , the quantityn · ∇V ( r) is constant.

B. Deformation Parameters
In practice, we consider three deformation degrees of freedom, namely, β 2 , β 4 and γ.

C. Woods-Saxon Parameters
The parameters [47] appearing in the Woods-Saxon potential are as defined below: 1. Central potential: a. Strength: with U 0 =53.754 MeV and U 1 =0.791.
c. Diffuseness parameter: assumed to be same for neutrons and protons, and has the value a = 0.637 fm.
b. Half density radius and diffuseness parameter are taken to be the same as those for the central potential.
The parameters have the isospin dependence of the central and spin-orbit potentials "built-in". This potential yields a reasonably good description of charge radii (both magnitude and isospin dependence) as well as of moments of inertia for a wide range of nuclei.
It has been used extensively in the total Routhian surface (TRS) calculations, and it has been quite successful in accurately reproducing energies of single-particle as well as collective states [48].

D. Coulomb potential
The Coulomb potential is calculated by folding the point proton density distribution ρ( r ′ ), assumed to be of Woods-Saxon form. For simplicity, its parameters are assumed to be the same as those for the nuclear potential of protons. The reason for using folded potential here is, as we have indicated in section II, the WK expansion is not valid for potentials with sharp surfaces.
The Coulomb potential for the extended charge distribution is given by: Here, where as explained in Appendix A. It is instructive at this point, to compare the Coulomb potential calculated from the diffuse density with the corresponding potential obtained by using the conventional uniform density (sharp surface) approximation. Such comparison for 208 Pb is plotted in Fig. 1. The radius parameter for the diffuse density approach as well as for the sharp surface approximation is assumed to be equal to 7.11 fm (see the discussion on the choice of the Woods-Saxon parameters in Section 3). It can be seen that in the exterior region, the two potentials agree almost exactly, as expected. In the interior, however, the potential obtained from the diffuse density turns out to be somewhat less repulsive than that from the density with sharp surface.

IV. DETAILS OF THE WK CALCULATIONS
In the present work, we restrict our calculations to three deformation degrees of freedom, namely, β 2 , β 4 and the angle γ. The inclusion of γ allows to incorporate triaxiality. Thus, the present WK calculation is genuinely three dimensional. In principle it is natural to use a cylindrical coordinate system here. The spherical polar coordinates, however, turn out to be more convenient. The reason is, the cylindrical coordinates involve two length variables, and one angular coordinate which means that the turning points have to be evaluated for two coordinates (ρ and z). This makes the calculations very complicated. On the other hand, the spherical polar coordinates involve only one length variable, and thus the turning points are to be evaluated only for one coordinate (r). The numerical integrals involved are evaluated using Gaussian quadrature.
The first step in the WK calculations is the determination of the chemical potential.
This has to be done iteratively, using Eq. (16). Since the turning points are determined by the chemical potential, they have to be calculated using a suitable numerical technique at each step. Once the values of the chemical potential are known, the WK energies up to second order can be calculated in a straightforward way. The fourth-order calculations are very complicated, since they require higher-order derivatives of nuclear potentials, spinorbit form factors, and the Coulomb potential. The former can be evaluated analytically in the present case. The expressions are extremely lengthy, and we do not present them here. Comparatively, the derivatives of the Coulomb potential look simple; the Laplacian and Laplacian of Laplacian are completely straightforward: the former is proportional to the proton density and the latter is just the Laplacian of the WS form factor. However, the calculations also need terms like Laplacian of the gradient squared of the total potential. In the case of protons, this involves one crossed term: where V C is the Coulomb potential and V N is the nuclear potential. The determination of such objects is tricky. It turns out that if one uses the form of the Coulomb potential defined above, the calculation of expression (42) becomes numerically unstable.
There exists an alternative for of the Coulomb potential: where the notation ∇ 2 r ′ means that the Laplacian is calculated with respect to the variables r ′ , θ ′ , and φ ′ . Eqs. (39) and (43)  It turns out that the WK calculations for protons are very time consuming. This is due to the fact that the calculation of Coulomb potential (Eq. (39)), in general, involves evaluation of three dimensional integral for each point (r, θ, φ). Typically, it takes few tens of minutes to complete one such calculation. This is certainly not desirable, since our aim is to calculate the masses of the nuclei spanning the entire periodic table. To speed up the calculations, we use the well-known technique of interpolation. Since we are using spherical polar coordinates, the turning points are to be evaluated only for the radial coordinate, r.
For the entire WK calculation, the θ and φ mesh points remain the same (over the domains [0, π] and [0, 2π], respectively), whereas the r mesh points change from step to step. This happens in particular during the evaluation of the chemical potential. Once the convergence of the particle number equation (Eq. (16)) is achieved, the r mesh points as well, remain fixed.
Motivated by the above observations, we apply the following procedure: 1. Before entering into the actual WK calculations (determination of chemical potential, etc.), for each pair of θ and φ mesh points, we calculate the Coulomb potential (Eq. 2. Next, for each pair of θ and φ mesh points, we fit a polynomial of degree 9 in the radial coordinate r to the Coulomb potential calculated in the above step. Thus, the fitting procedure is to be repeated N θ ×N φ times, N θ (N φ ) being total number of mesh points for the θ (φ) integration.
Thus, for any given value of radial the coordinate r (and fixed θ and φ), the Coulomb potential can be easily calculated just by evaluating the 9 th degree polynomial in r. It is found that this interpolation procedure is very accurate. The maximum percentage difference between the fitted and the exact Coulomb potentials is 0.4% for a highly deformed nucleus.  Fig. (2). The other two deformation parameters, β 4 and γ, are set to zero in this test case. The partial contributions to the WK energy are plotted separately for protons and neutrons. It is found that all the correction terms vary smoothly as a function of deformation. As expected, the value of the contributions from the 2 and 4 terms to the averaged energy decreases rapidly. It is found that the proton and neutron Thomas-Fermi energies have opposite trends with respect to increasing β 2 . If Coulomb potential is suppressed, it is found that the Thomas-Fermi energies for protons follow the same trend as those for the neutrons. Further, it is interesting to note that comparatively, the variation in the second-order corrections with respect to deformation parameters is stronger than that in the Thomas-Fermi energies (∼ 10% for second-order corrections and ∼ 3% for Thomas-Fermi energies).
Next, the variation of the Thomas-Fermi energy and of the correction terms as a function of the hexadecapole deformation parameter β 4 is plotted in Fig. (3). Here, β 2 is taken to be 0.2 and γ is set to zero. It is seen that again, the different energies vary smoothly as a function of β 4 . The Thomas-Fermi energy for protons is found to have very little variation with respect to the β 4 deformation parameter. In contrast, the corresponding energies for neutrons have a stronger dependence on β 4 . The same behaviour is also observed in the corresponding quantum mechanical energies. It is found that the proton and neutron Thomas-Fermi energies have a very similar behaviour if the Coulomb potential is suppressed.
Further, to check if this conclusion depends on the value of β 2 , the analysis is repeated for β 2 = 0.4, and the same conclusion is found to emerge.
The behaviour of the Thomas-Fermi energies for protons in the above cases (Figs. (2) and (3)) seems to be due to the Coulomb potential. In the case of variation with respect to β 2 , qualitatively it can be expected that with increasing quadrupole deformation, protons are pulled apart and Coulomb repulsion decreases, thereby making the system more bound.
The β 4 deformation also affects the proton distribution, but, as expected, the effect of hexadecapole deformation is less prominent in comparison with that of quadrupole deformation.
Thus, the repulsion among protons does decrease with increasing β 4 , but the decrease is not large enough to make the system more bound with larger β 4 .
By keeping β 2 and β 4 fixed, if the parameter γ is varied, then it is found that the resulting energies are independent of the sign of γ. Moreover, the γ dependence of the WK energies is found to be rather weak. Therefore, here we do not present these result explicitly.
The fourth-order calculation for protons is very time consuming. Typically, it takes tens of minutes to do a complete WK calculation. Most of the run-time being consumed by particle number determination and the fourth-order calculations for protons. Thus, it is necessary to find an accurate approximation scheme for the fourth-order calculation for protons. Since in the nuclear interior, the Coulomb potential has approximately a quadratic nature (see Fig. (1)), it is expected that the Coulomb potential will have small influence on the fourthorder calculations (note that one needs higher-order derivatives in the fourth-order energy calculations). One may therefore drop the Coulomb potential completely from the fourthorder corrections; we shall refer to this approximation as "quadratic approximation". This approximation has been checked explicitly by performing exact fourth-order calculations for protons. The maximum difference between the WK energies obtained by using exact calculation and the quadratic approximation is found to be of the order 100 keV for 82 protons. It turns out that the difference between the quadratic approximation and exact calculation decreases with decreasing charge number. This approximation can be improved by keeping the Laplacian of the Coulomb potential in the fourth-order contribution i.e., the terms of the form (∇ 2 V ) 2 and ∇ 4 V in Eq. (17). This means that for protons, only the term ∇ 2 (∇V ) 2 is dropped from Eq. (17). It is found that with this modification, the value of the fourth-order correction energy for the mean field part for protons almost coincides with the value obtained by taking all of the derivatives of the Coulomb potential into account. This helps in reducing the total runtime further. Thus, effectively, with the interpolation for Coulomb potential as discussed before (see section IV), and the approximations introduced in the fourth-order correction terms for protons in the present section, the runtime reduces from tens of minutes to just about two minutes, without affecting the desired accuracy of the calculations.

VI. WIGNER-KIRKWOOD SHELL CORRECTIONS AND COMPARISON WITH STRUTINSKY CALCULATIONS
Numerically, it has been demonstrated that the WK and Strutinsky shell corrections are close to each other [20]. This is expected, since it has recently been shown [50] that the Strutinsky level density is an approximation to the semi-classical WK level density.
For illustration, we present and discuss the WK and the corresponding Strutinsky shell corrections for the chains of Pb, Gd and Dy isotopes. For the sake of completeness, we first present and discuss the essential features of the Strutinsky smoothing scheme.
According to the Strutinsky smoothing scheme, the smooth level density for a one-body Hamiltonian is given by [49]: where ǫ i are the single-particle energies calculated by diagonalising the Hamiltonian matrix.
The smoothing constant γ is taken to be of the order of ω 0 ( ω 0 = 1.
is the smoothing order, and is assumed to be equal to 6 in the present work; H j are the Hermite polynomials; and S j is a constant, defined as [49]: The Strutinsky shell correction is given by: where N n is the number of nucleons. This, upon substituting the expression for g st , yields [49], The error integral erf(x) is defined as: It should be noted that the Strutinsky procedure described here uses positive energy states generated by diagonalizing the Hamiltonian matrix, and not by taking resonances into account and smoothing them. Further, in practice, the summations defined above do not extend up to infinity, but are cut off at a suitable upper limit. The limit is chosen in such a way that all the states up to ∼ 4 ω 0 are included in the sum. It has been shown that the uncertainty in the Strutinsky shell corrections obtained this way is typically of the order of 0.5 MeV [49]. For lighter nuclei, however, it has been concluded [49] that this uncertainty is larger.
The total WK shell correction for the chain of even even Lead isotopes ( 178−214 Pb) is plotted in Fig. (4), along with the corresponding values obtained by using the Strutinsky smoothing method. It is found that both the WK and Strutinsky results exhibit very similar trends. As expected, there is a prominent minimum observed for 208 Pb, indicating the occurrence of shell closure. The WK and Strutinsky shell corrections slightly differ from each other. The difference is not a constant, and is found to be increasing slowly towards the more neutron deficient Lead isotopes.
Next we plot the calculated (WK) and the corresponding Strutinsky shell corrections for the chains of even even Gd and Dy isotopes, with neutron numbers ranging from 72 to 92. Apart from 144,146,148 Gd and 146,148,150 Dy, the rest of the nuclei considered here are known to be deformed [12]. For this test run, we adopt the deformation parameters from the Möller -Nix compilation [12]. It is seen that the WK and the corresponding Strutinsky shell corrections agree with each other, within few hundred keVs. The prominent minimum at shell closure at neutron number 82 is clearly seen. In these cases as well, the difference between the two calculations is not a constant. It is larger in the neutron deficient region, and becomes smaller as neutron number increases. by determining λ 1 and λ 2 using certain conditions. Here,Ĥ is the pairing Hamiltonian, and N is the particle number operator. Minimisation of the expectation value ofĤ − λ 1N leads to the usual BCS model, with λ 1 determined from the particle number condition. Thus, in Eq. (50) above, the quantity λ 1 is a Lagrange multiplier, but the particle number fluctuation constant λ 2 is not.
In practice, the LN calculation is carried out by assuming a constant pairing matrix element, G. For a given nucleus (assumed to be even-even for simplicity), one considers N h doubly degenerate states below, and N p doubly degenerate states above the Fermi level.
These states contain N nucleons. In practice, one takes N h = N p = N/2 or Z/2, depending on whether it is being applied to neutrons or protons. The occupation probabilities v 2 k , the pairing gap ∆, the chemical potential λ (= λ 1 + 2λ 2 (N + 1), see Ref. [31]), and the constant λ 2 are determined iteratively using the conditions [13,31]: such that and where E k are the single-particle energies and u 2 k = 1 − v 2 k . The particle number fluctuation constant λ 2 is given by: The pairing matrix element G is calculated by the Möller-Nix prescription [13]: Here,ρ L = g W K /2 is the Wigner-Kirkwood averaged level density (see Eq. (15). Factor of 2 appears because each quantal level here has degeneracy of 2. The level density is evaluated at fermi energy.); a 2 = N /2ρ L and a 1 = −N /2ρ L and∆ is the average pairing gap, taken to be 3.3/N 1/2 [13].
The ground-state energy within the LN model is given by: The pairing correlation energy, E pair is obtained by subtracting the ground-state energy in absence of pairing from Eq. (57):

VIII. CALCULATION OF BINDING ENERGIES
As an illustrative example, we now present and discuss the calculated binding energies (in this paper, we take binding energies as negative quantities) for 367 even-even, evenodd, odd-even and odd-odd spherical nuclei. These nuclei are predicted to be spherical or nearly spherical (β 2 < 0.05) in the Möller-Nix calculations [12] and include 38 ). Of course, it is known that the prediction of sphericity does depend to some extent on the details of the density functional employed [52]. Therefore, it may so happen that some of the nuclei assumed to be spherical here, may actually turn out to be slightly deformed when energy minimization is carried out on the grid of deformation parameters.
Our calculation proceeds in the following steps. For each nucleus, the quantum mechanical and WK energies are calculated as described earlier. This then yields values of the shell corrections (δE) for these nuclei. The pairing energies (E pair ) are then calculated using the Lipkin-Nogami scheme [29,30,31] described previously in the same potential well where the shell correction is computed. These two pieces constitute the microscopic part of the binding energy. The macroscopic part of the binding energy (E LDM ) is obtained from the liquid drop formula. Thus, for a given nucleus with Z protons and N neutrons (mass number A = N + Z), the binding energy in the Mic-Mac picture is given by: The liquid drop part of binding energy is chosen to be: where the terms respectively represent: volume energy, surface energy, Coulomb energy and correction to Coulomb energy due to surface diffuseness of charge distribution. The coefficients a v , a s , k v , k s , r 0 and C 4 are free parameters; T z is the third component of isospin, and e is the electronic charge. The free parameters are determined by minimising the χ 2 value in comparison with the experimental energies: where E(N j , Z j ) is the calculated total binding energy for the given nucleus, E (j) expt is the corresponding experimental value [53], and ∆E (j) expt is the uncertainty in E (j) expt . In the present fit, for simplicity, ∆E (j) expt is set to 1 MeV. The minimisation is achieved using the well-known Levenberg-Marquardt algorithm [54,55].   Table I. Clearly, the obtained values of the parameters are reasonable. The detailed table containing the nuclei considered in the present fit, and the corresponding calculated and experimental [53] binding energies may be found in Ref. [51].
To examine the quality of the fit further, first, we plot the difference between the fitted and the corresponding experimental [53] binding energies for the 367 nuclei as a function of the mass number A in Fig. (5) The corresponding differences obtained for the Möller-Nix Next, the difference between the calculated and the corresponding experimental [53] binding energies (denoted by "WK") for Ca, Ti, Sn, and Pb isotopes considered in this fit are presented in Fig. (6). The differences obtained by using the Möller-Nix [13] values of binding energies (denoted by "MN") are also shown there for comparison. It can be seen that the present calculations agree well with the experiment. It is found that the differences vary smoothly as a function of mass number: the exceptions being the doubly closed shell nuclei 48 Ca, 132 Sn, and 208 Pb, where a kink is observed. The overall behaviour of the differences is somewhat smoother than that obtained by using the values of Möller and Nix. To investigate the effect of the parameters of the single-particle potential, we make a refit of the liquid drop parameters, by using the Rost parameters [56] in the microscopic part Thus, overall, this potential is more realistic than the Rost potential. This is reflected in the calculated binding energies as well, showing clearly that the choice of the single-particle potential (or in other words, the parameters) is indeed important for reliable predictions of binding energies (and hence the masses).
Single and two neutron separation energies (S 1n and S 2n ) are crucial observables. They are obtained by calculating binding energy differences between pairs of isotopes differing by one and two neutron numbers, respectively. The single neutron separation energies govern asymptotic behaviour of the neutron density distributions [57]. They exhibit oddeven staggering along an isotopic chain, indicating that the isotopes with even number of neutrons are more bound than the neighbouring isotopes with odd number of neutrons. The systematics of S 2n primarily reveals the shell structure in an isotopic chain. The correct prediction of these separation energies is crucial for determination of the neutron drip lines.
The calculated S 1n and S 2n values for Sc, Sn and Pb isotopes are displayed in Fig. (7). The corresponding experimental values of S 1n and S 2n [53] are also plotted for comparison. The agreement between calculations and experiment is found to be excellent. The odd -even staggering is nicely reproduced. The shell closures at 132 Sn and 208 Pb are clearly visible both in single and two neutron separation energies. At a finer level, however, a marginal underestimation of the shell gap at the neutron number 82 (126) is observed in 132 Sn ( 208 Pb).
Finally, we remark that the calculated single and two proton separation energies are also found to be in close agreement with the experiment.
The results presented in this section indicate that the present calculations of binding energies, indeed, are reliable.

IX. SUMMARY AND FUTURE OUTLOOK
In the present work, we intend to carry out reliable mass calculations for the nuclei spanning the entire periodic table. For this purpose, we employ the 'microscopic-macroscopic' framework. The microscopic component has two ingredients: the shell correction energy and the pairing energy. The pairing energy is calculated by using the well-known Lipkin-Nogami scheme. To average out the given one-body Hamiltonian (and hence find the shell corrections, given the total quantum mechanical energy of the system), we use the semiclassical Wigner-Kirkwood expansion technique. This method does not use the detailed single-particle structure, as in the case of the conventional Strutinsky smoothing method.
In addition to the bound states, the Strutinsky scheme requires the contributions coming in from the continuum as well. Treating the continuum is often tricky, and in most of the practical calculations, the continuum is taken into account rather artificially, by generating positive energy states by means of diagonalisation of the Hamiltonian matrix. For neutronrich and neutron-deficient nuclei, the contribution from the continuum becomes more and more important as the Fermi energy becomes smaller (less negative). Uncertainty in the conventional Strutinsky scheme thus increases as one goes away from the line of stability.
It is therefore expected that the Wigner-Kirkwood method will be a valuable and suitable option especially for nuclei lying far away from the line of stability.
We now summarise our observations and future perspectives: Therefore, before performing the large scale calculations, we intend to make a refit to the Woods-Saxon potential, with the Coulomb potential obtained from folding.

5.
Having established the feasibility of the present approach, we now intend to extend our binding energy calculations to deformed nuclei. For this purpose, we plan to minimise the binding energy on a mesh of deformation parameters to find the absolute minimum in the deformation space. Work along these lines is in progress.

APPENDIX A: GEOMETRY OF DISTANCE FUNCTION
Consider an arbitrary surface, defined by the relation r = r s , where r s is given by Eq.
(21) of the text. Let us fix the origin of the coordinate system at the centre of mass of the object. Let r ≡ (r, θ, φ) define an arbitrary point in space. This point could be inside or outside the surface. Here, for concreteness, we assume that it is within the volume of the object. Our aim is to find the minimum distance of the point r to the surface r = r s .
To achieve this, construct a vector R s from the centre of mass to the surface. To find the minimum distance, one has to minimise the object | r − R s |. Denoting the angle between the r and R s vectors as Ψ, we have: where, from Fig. (A), the cosine of the angle Ψ is given by: cos Ψ = cos θ cos θ s + sin θ sin θ s cos(φ s − φ). The Coulomb potential for an arbitrary charge distribution is given by: Let, for brevity, | r − r ′ | = R. Consider: Here, the symbol ∇ r ′ means that the differentiation is done with respect to the r ′ , θ ′ , φ ′ coordinates. Let us consider the above derivative component-wise. The contribution coming from the first component is: Adding contributions coming from all the three components, one gets: With this, the potential becomes: Here, the term in the curly brackets has been represented as a unit vectorR. Using the identity: one obtains: which, upon integrating by parts and transferring derivatives to density, becomes: q.e.d.

Derivatives of Coulomb Potential
The calculation of the higher-order derivatives of the Coulomb potential (third and above), even with the form defined in Eq. (43), turns out to be numerically unstable.
For this purpose, we employ the Poisson's equation. According to this, the Laplacian of Coulomb potential is proportional to the charge density: The Laplacian of ∇ 2 V C ( r) is simple to compute, for, all one needs to calculate there are the derivatives of density (assumed to be of Woods-Saxon form).
Thus, it is desirable to generate the required higher order-derivatives of the Coulomb potential (see expression (42) in the text) from Poisson's equation. For this purpose, we evaluate the commutators: The results are: With these expressions, the required higher-order derivatives of the Coulomb potential can be generated. These are then used to evaluate the fourth-order WK energy, as we have described in Section 4.