Semiclassical evaluation of average nuclear one and two body matrix elements

Thomas-Fermi theory is developed to evaluate nuclear matrix elements averaged on the energy shell, on the basis of independent particle Hamiltonians. One- and two-body matrix elements are compared with the quantal results and it is demonstrated that the semiclassical matrix elements, as function of energy, well pass through the average of the scattered quantum values. For the one-body matrix elements it is shown how the Thomas-Fermi approach can be projected on good parity and also on good angular momentum. For the two-body case the pairing matrix elements are considered explicitly.


I. INTRODUCTION
The solution of the nuclear many-body problem presents a formidable challenge. Not only bare and effective nucleon-nucleon forces are not completely known, but still for those given as granted one has to solve the many-body problem of a highly quantal, strongly interacting, self-bound, and therefore inhomogeneous Fermi system. Over the years semiclassical techniques have helped to solve this problem, for instance in regard to the latter aspect. In practice it is mainly the Thomas-Fermi (TF) method and its extensions for the description of nuclear ground-state properties which has been considered (see [1] and references therein). The nuclear density and kinetic energy density are the main ingredients of this approach.
The semiclassical approximation often gives a direct physical insight, yielding the shell average of the quantities under consideration and providing their main trend (which, in certain cases, may be obscured by strong shell fluctuations). A known example is e.g. the nuclear binding energy which coincides with the liquid drop part in the semiclassical approach. Another quantity of longstanding interest is the average single-particle level density. It is well known that the TF approximation to the level density (includingh corrections) coincides analytically with the Strutinsky averaged quantal level density for the harmonic oscillator (HO) potential [2]. First performing the quantal calculation and then the average is more cumbersome than calculating the shell average directly via the TF method. The technical advantage of the latter becomes significant in the deformed case [2] or, for instance, when one wants to go beyond the independent particle description to include correlations [3]. The TF approach is also very helpful for the calculation of surface and curvature energies. Actually, the latter quantity can only be correctly extracted in a semiclassical procedure [4].
In general, however, we believe that the true virtue of the TF method shows up not only in calculating average properties in the independent particle approximation, where it can replace the results obtained through the more cumbersome Strutinsky method [5], but rather in many-body applications going beyond the mean field or independent particle picture where a straightforward quantum solution for finite systems may reach its limits. A case where we treated correlation effects in TF approximation was, as already mentioned, the level density parameter [3]. Pairing correlations in finite nuclei have also already successfully been treated in the past [6]. This is one of the aspects which we shall consider again in this work in more detail.
In this work we want to dwell on an aspect of Thomas-Fermi theory which in the past has been exploited only very little. This concerns the evaluation of matrix elements averaged over a certain energy interval which may be typically of the order ofhω, i.e., the separation of major shells. It must be pointed out that this is, to our knowledge, the first attempt to evaluate not only one-body but also two-body matrix elements in the TF approximation. This semiclassical calculation provides the smoothly varying part of the matrix elements dropping the shell effects according to the idea of the Strutinsky averaging method [7]. We will describe several tests of the accuracy of the TF method for on-shell densities. In a first part we will develop our approach for the matrix elements of single-particle operators, for given parity and angular momentum. This goes along similar lines already developed in the domain of systems with chaotic behaviour [8][9][10][11]. In a second part we will address our main objective, which is to show that the method also works for two-body matrix elements. Some preliminary results have been published previously in Ref. [12]. As a specific example we will treat the pairing matrix elements.
Let us give a short summary of the approach we are going to develop. Consider for example the expectation value of a single particle operatorÔ in some shell model state |ν : Instead of knowing O ν quantum state by quantum state it may some times be advantageous and instructive to only know how the matrix element (1) changes as a function of energy. We therefore introduce a single-particle matrix element averaged over the energy shell: where we callρ E the density matrix on the energy shell. It is related with the so-called spectral density matrix and will be defined immediately below. The spectral density matrix δ(E −Ĥ) has the characteristic discontinuous behaviour due to the quantization of the eigenvalues of the single-particle HamiltonianĤ. It can be written, however, as a sum of a smooth partδ(E −Ĥ) and of a strongly oscillating part, i.e., δ(E −Ĥ) =δ(E −Ĥ) + δ osc (E −Ĥ) [1,13,14]. Analogously, the single-particle level density g(E) = T r[δ(E −Ĥ)] = ν δ(E − ε ν ) is obtained as a sum of two terms: g(E) =g(E) + g osc (E), whereg and g osc stand for the smooth and the rapidly fluctuating contributions, respectively. Using the smoothδ(E −Ĥ) andg(E), we define the density matrix averaged on the energy shell aŝ It is therefore a smooth function of E. The smeared level densityg(E) (per spin and isospin in this paper) in the denominator of expression (3) ensures the right normalization ofρ E , The smooth quantities entering eq.(3) are to be evaluated in some continuum limit [13,14]: this is the case for example when one introduces the Strutinsky averaging procedure [5,7] or, alternatively, and this is the approach we adopt in this work, it can be done by replacinĝ H, the independent-particle (mean field) Hamiltonian, by its classical counterpart H cl which corresponds to the TF approximation [1]. Such an approximation has been used very early by Migdal [15] and later, as already mentioned, in the context of chaotic motion dynamics [8][9][10][11]. Recently, we have employed it to describe Bose condensates in traps [16], but we are not aware of any systematic use in the context of nuclear physics. The approach is not limited to the evaluation of expectation values of single particle operators. Also the average behavior of two-body matrix elements can be calculated. For instance the semiclassical evaluation of the average pairing matrix element where |Φ(ν,ν) is an antisymmetric normalized two-body state constructed out of a state |ν and its time-reversed state |ν , can be of great practical interest and shall be considered in this work. As it is known [1,2], the Strutinsky method averages the density matrix over an energy interval corresponding roughly to the distance between two major shells. Implicitly the same holds if the equivalent Wigner-Kirkwood expansion (TF approximation at lowest order) is used forρ E . In detail our paper is organized as follows. In the next section the Wigner function on the energy shell is introduced and applied to the evaluation of some semiclassical one-body and two-body matrix elements of physical interest. Our conclusions are given in the last section.

II. WIGNER FUNCTION ON THE ENERGY SHELL
As stated in eq. (3) of the introduction, we are interested in the density matrixρ E on the energy shell. Consider the Wigner transform [1] of the density matrix (3), namely where are the centroid and relative coordinates respectively. In order to obtain the TF approximation plush corrections to the Wigner function on the energy shell (5), it is convenient to differentiate with respect to E the Wigner-Kirkwood expansion of the full single-particle one-body density matrixρ = Θ(E −Ĥ), which is amply given in the literature [1]. Up to orderh 2 the result is One should realize thatg(E) also containsh corrections and that, strictly speaking, in order to get a consistent expansion of f E one should also take into account theh expansion of g(E) and then correctly sort out relation (6) to orderh 2 (see also comments at the end of Section II.A). The first term in eq. (6) represents evidently the pure TF approximation which is of lowest order inh. In a first attempt and to assess the accuracy of our approximation we will content ourselves with the TF approximation. Integration over the momenta yields the local density on the energy shell: where is the local momentum at the energy E in the potential V (R) and the level densityg(E) is given byg For the following it is important to first elaborate on the meaning and accuracy of this density on the energy shell. For demonstration purposes we will take as an example the spherical HO potential but later we will see that our method works equally well for a Woods-Saxon (WS) potential.
In fig. 1 we display the quantal (solid line) and TF (dash-dotted line) densities of the N = 4 and N = 5 HO shells withhω = 41A −1/3 and A = 224. For the TF densities we have taken the quantal energies. We see that in both cases the TF result passes accurately through the average, terminating at the classical turning point defined by The features of the TF densities on the energy shell are quite analogous to the ones already known for the full TF density [1]. However, on the quantal side one remarks that there is a strong difference between odd parity (N = 5) and even parity (N = 4) shells. The former shows a pronounced hole at the origin whereas the second ones shows, on the contrary, an enhancement. Both features can obviously be related to the absence or presence of the swave contributions in the corresponding HO shell, respectively. One may try to recover this even-odd parity effect in projecting the TF density matrix on good parity. This is easily done as follows. We calculate the inverse Wigner transform of f T F E (R, p). This yields where j 0 is the zeroth-order spherical Bessel function. Now the even/odd parity density on the energy shell is obtained as We have drawn this expression in fig. 1 (dashed lines) as well. The bump (hole) structure exhibited by the quantal density is now well reproduced in the interior. The agreement only deteriorates near the classical turning point. One should mention, however, that in spite of the seemingly rather spectacular improvement of formula (12) over (7), the former presents some small problems. This concerns the behavior of (12) around the turning point. The presence of the second term in (12) can induce a slightly negative value of the density around the turning point. Also the second term is not naturally limited to r values inside the turning point and thus is oscillating around zero due to the Bessel function. This leads to ambiguities in evaluating matrix elements such as (2) which, however, numerically are rather unimportant. Thus we advocate to use (7) instead of (12), excepting for some problems where the even/odd bump structure may be particularly important. The latter may for example be the case for the evaluation of matrix elements of operators more concentrated at the nuclear interior like e.g. 1/r 2 , etc.

A. One-body matrix elements
We now proceed to calculate in TF approximation as a function of the energy the rms radius of a nucleon confined in a WS potential with V 0 = −44 MeV, a= 0.67 fm and R = 1.27A 1/3 fm with A=224 nucleons. We choose the rms radius for demonstration purposes but we could have taken as well any other smoothly varying single-particle operator. We use the TF approximation (11) and show the results (dashed line) together with their quantum mechanical counterparts, represented by dots, in fig. 2. We see that the TF calculation very nicely passes through the average of the scattered quantal values, with the exception of the lowest s-state. This is a first confirmation of the accuracy of our approach.
In a next step we want to project the TF density matrix on different partial waves and calculate matrix elements as a function of the energy for different l-values. One way to project on partial waves has been elaborated by Hasse [17]. There one pre-multiplies the Wigner function with the semiclassical projectors on the orbital angular momentum and its z-component, i.e., Then one can calculate single particle matrix elements as For local operators it is sufficient to know the density, which can be obtained in integrating (13) over momenta. Assuming spherical symmetry we can also sum over the m quantum numbers and after some algebra one finally finds [17]: where the level densityg E,l is chosen in such a way that the integral of (15) over the available R space is normalized to unity. Let us mention that we have replaced l(l + 1) by (l + 1 2 ) 2 as it is done in the WKB method to recover the right asymptotic behaviour of the wave function in the free (V (R) = 0) case [18].
We recently employed, however, a different way to do the l-projection which for some purposes may be more convenient [19]. For this we first perform the inverse Wigner-transform of the TF part of eq. (6): Expanding the plane waves in spherical harmonics, we can read off the l-projected density matrix. For the local density we then obtain whereg E,l,m is again chosen to ensure the right normalization of (18), i.e., dRρ E,l,m (R) = 1. We can calculate for example the mean square radius as a function of E for different l values: Notice that eq.(19) becomes independent of m after the angular integration.
In table I we show various moments R n 1/n (in fm) obtained using the TF local densities on the energy shell provided by eqs. (15) and (18), as compared with the corresponding quantal values for the aforementioned WS potential with A = 224 nucleons. From these tables one can see that the quantal rms radius (n = 2) of each nl state is, on average, well reproduced using the semiclassical l-projected on-shell densities given by the TF approaches eqs. (15) and (18), where E has been replaced by the quantal eigenvalue corresponding to the nl state (for consistency we should perform a WKB quantization from which we refrain for simplicity as it should not affect the result by much). A more quantitative analysis shows that the quantal rms radii are reproduced by eq.(15) within 2% in average over the range of energies considered, while using eq.(18) the average relative error is around 4%. Higher moments obtained with the TF densities on the energy shell also reproduce reasonably well the results of the full quantal calculation. For the highest moment considered here (n = 10), the average relative error with respect to the quantal calculation is around 4.3% for both semiclassical approximations. One should point out that the TF local densities on the energy shell are free of the shell effects that are present in the quantal calculation. Actually, the TF results represent the shell averaged values of the moments and their difference with the quantal calculations provides an estimate of the shell correction for the considered state. Of course, a precise calculation of the shell correction in the considered moments would require a Strutinsky calculation which is not an easy task for a WS potential. For a related discussion we refer the reader to [20], where the moments ranging from n = 1 to n = 10 of the total density of A = 224 particles in a HO potential were evaluated using several semiclassical approaches and the Strutinsky averaging method. The fact that TF works for moments as high as n = 10 makes us believe that one can certainly consider in our approach all single-particle operators which are low-order polynomials in the phase-space variables.
It is also instructive to directly compare the densities. For this we again take E equal to its quantal value in eqs. (15) and (18). The comparison between the quantal and TF on-shell densities is shown in figs. 3-5. In each one of these figures the quantal on-shell density (solid line) and the semiclassical densities provided by eq.(15) (dash-dotted line) and eq.(18) (dashed line) are displayed for all the bound s and p states of a WS potential with A = 224. From these figures one can see that the quantal on-shell densities are rather well reproduced by the TF approach eq.(18). In particular the quantal 1l on-shell densities are well reproduced by eq.(18) in the interior. However, as it happens in fig. 1, there are discrepancies between the quantal and TF densities in the outer part around the classical turning point. For nl states (n > 1), the inner n − 1 bumps are well reproduced by eq. (18) and the agreement deteriorates in the outer bumps due to the classical turning point where the TF densities vanish. On the other hand the TF on-shell densities obtained with eq.(15) do not reproduce the quantal density profiles at all. As it can be seen from eq.(15), this density is defined in the region in between the two roots (turning points) of In this approach the three-dimensional problem has been reduced to an equivalent onedimensional problem for R with an effective potential that in addition to V (R) contains the centrifugal barrier, as it happens in the WKB method. Thus, in this case we find two different turning points. The largest root of (20) is very close to the classical turning point of eq.(18) given by k E = 0 [see eq.(8)], while the smallest root gives the inner turning point due to the centrifugal barrier. Since the TF on-shell density (15) has square root singularities at the two turning points, its integral as well as the corresponding expectation values converge. We arrive at the, at first sight, paradoxical result that the densities (15) which have no detailed resemblance with the quantal ones reproduce the rms values (and very likely most of other expectation values of smoothly varying operators) better than the densities given in (18), which show quite reasonable overall behaviour in comparison with the quantal results. We here find an illustrating example that the Thomas-Fermi and Wigner-Kirkwood local densities are to be regarded as mathematical distributions, in the sense that in spite of their possible divergences they yield finite, accurate results when used to compute a restricted class of expectation values by integration [13,20]. In this respect we again refer the reader to tables I and II where he finds confirmation of our statement. A similar situation is found in computing the kinetic energy for a bosonic system in Ref. [16]: the semiclassical and quantal kinetic energy densities are clearly different but give similar values of the integrated kinetic energy. On the other hand, one should note that sorting out correctly the various orders inh is very important to achieve optimal results as shown on other occasions [20,21]. In expression (18) there remain someh corrections to all orders in the form of the spherical harmonics which are quantal wave functions, i.e., solutions of a Schrödinger equation. This mixing of resummation inh on the one hand and of lowest order inh on the other hand [in the form of the δ-function in (16)] finally makes the on-shell density (18) slightly less accurate than (15), which represents the correcth → 0 limit as shown in Ref. [17].

B. Two-body matrix elements
As a further interesting application we want to consider the semiclassical evaluation of average two-body matrix elements. An example of particular interest is e.g. the case of matrix elements of the pairing type Φ(ν,ν)|v|Φ(ν ′ ,ν ′ ) [see eq. (4)], which we shall address in this section for identical nucleons. It is straightforward to recast the state-dependent pairing matrix element as The two-particle states |νν on the r.h.s. are product states of |ν and |ν . The states |ν are represented by single-particle wave functions φ ν (r, σ) = φ nlm (r) χ σ (with σ = ± 1 2 ). Assuming spherical symmetry, and considering that the time reversal of |ν involveŝ According to this result, we obtain the following expression for the average pairing matrix elements of eq. (4): where ρ E (r, r ′ ) = r|ρ E |r ′ . In TF approximation the non-local on-shell density matrix ρ E (r, r ′ ) is given by eq. (11). We see that it is a symmetric function in r and r ′ . For the case of a force v 0 δ(r − r ′ ) eq.(23) reduces to (a practically identical expression can be found in Ref. [22]) Using the TF expression (7) for ρ E (r) we can evaluate (24) with a HO potential V (r) = mω 2 0 r 2 /2 and compare with the quantum mechanical matrix elements averaged on each major shell N of energy E = E N = (N + 3 2 )hω. This is done in Table II with v 0 = −345.723 MeV fm 3 andhω 0 = 41A −1/3 MeV. We again see that the semiclassical results agree very well with the averaged quantal values, even for the non-diagonal elements.
With the above positive experience at hand, next we proceed to calculate the average pairing matrix elements v(ε F , ε F ) of the Gogny D1S force [23] which is known to reproduce experimental gap values when used in microscopic Hartree-Fock-Bogolyubov calculations [24]. Writing the diagonal matrix element (4) at the Fermi energy ε F , by using (23), and expressing it through the lowest-order Wigner function in inverting (5), one arrives in TF approximation at (25) Here, H cl = p 2 /2M * + V (R) is the classical Hamiltonian of independent particles with effective mass M * (see below) moving in an external potential well V (R), and v(p − p ′ ) is the Fourier transform of the particle-particle part of the Gogny force which describes the pairing. For the numerical application we use for V (R) a slight variant of the potential given by Shlomo [14]: with R 0 = 1.12A 1/3 + 1 [1 + (πd/R 0 ) 2 ] 1/3 fm, d = 0.70 fm, and In this equation R 0 has to be determined iteratively.
Equation (25) can be reduced to a one dimensional integral over R which can be performed numerically: where R c is the classical turning point defined in (10) and The factors z c correspond to pairing in the S = 0 and T = 1 channel and are written in terms of the parameters of the Gogny force W c , B c , H c , M c and µ c [23]. We have introduced the position-dependent effective mass M * (R) from the Gogny force in order to make the calculation of the pairing matrix element more realistic. It is obtained by evaluating at k = k ε F (R) the position-and momentum-dependent effective mass [25] M M * (R, k) where U(R, k) is the Wigner transform of the single-particle potential obtained from the Gogny interaction assumed spherically symmetric in k. Also the level densityg(ε F ) is calculated in the TF approach using the same potential and effective mass:g In fig. 6 we show Av(ε F , ε F ) as a function of the mass number A. The Coulomb force has been switched off in the present calculation. We see that there is a quite pronounced A dependence which is somewhat in contradiction with the value G ∼ 28/A MeV for the constant pairing matrix element at the Fermi level usually employed in more schematic treatments of the nuclear pairing. On the other hand, if we calculate Av(ε F , ε F ) not with a WS potential but with the HO potential, we obtain a practically constant value. This may indicate that the A −1 dependence of the constant pairing matrix element is better fulfilled in conjunction with a harmonic potential. The difference in the behaviour with A using the HO and WS potentials may come from the absence (HO) or presence (WS) of a surface contribution to v(ε F , ε F ), like it is the case for the level density [1].

III. CONCLUSIONS
In this work we showed how average nuclear one-and two-body matrix elements can very efficiently be evaluated using the Thomas-Fermi approach. The main ingredient is to replace the density matrix for a given quantum stateρ ν = |ν ν| by its counterpart averaged over the energy shellρ with simultaneous application of the Wigner-Kirkwoodh expansion for the smoothly varying spectral densityδ(E −Ĥ) and level densityg(E) = T r[δ(E −Ĥ)]. We calculated one-and two-body matrix elements restricting ourselves, in this exploratory work, to the lowest order, i.e., the pure Thomas-Fermi approximation. We compared quantal and semiclassical values of the matrix elements using harmonic oscillator and Woods-Saxon type of potentials. We did this also for parity projected and for angular momentum projected Thomas-Fermi theory. In all the cases close agreement with the average quantal behaviours was found, showing the accuracy of the method. As in the case of the well tested Wigner-Kirkwood expansion of the full density matrix [1], one expects that some improvement also could be achieved for the matrix elements by inclusion ofh corrections.
With the positive result for the single-particle matrix elements at hand, we also calculated the average pairing matrix elements of some effective nuclear two-body forces. First, we used a delta interaction and compared diagonal and off-diagonal semiclassical elements with the corresponding quantal values. Again the TF values nicely reproduce the quantum results on average. Next we estimated semiclassically the diagonal pairing matrix elements of the Gogny D1S force at the Fermi energy. Since the Gogny force is known to reproduce very well nuclear pairing properties [24], it is interesting to evaluate e.g. the A dependence of its pairing matrix elements around the Fermi energy and to see to which extent the common assumption of a A −1 dependence holds. Using for the mean field a potential of the Woods-Saxon type, it turned out that the fall off of the pairing matrix element is stronger than the A −1 law. For this problem we have no comparison with quantum values available, but the experience with the one-body matrix elements and the pairing matrix elements for the delta force, makes us believe that the values shown in fig. 6 are also reasonably accurate. The stronger than A −1 decrease observed in fig. 6 has its origin very likely in the presence of a surface contribution implicit from the use of a WS potential, whereas the use of a HO potential with the absence of a surface shows agreement with the A −1 law.
It should be emphasized that to our knowledge the TF method to calculate average matrix elements of a two-body force has never been applied before. We think that the conclusive study of this work will allow to use average matrix elements for the calculation of many nuclear quantities where fine shell effects are not needed such as optical potentials, giant resonances and their widths, and many other quantities where the average trend is of interest. In a future publication we will show how the application of these techniques can be used to study in a very transparent way the size dependence of the average pairing gap in finite Fermi systems in an almost analytical way.  TABLES   TABLE I. Moments R n 1/n (in fm) of several nl states for A = 224 particles in the Woods-Saxon potential described in the text. The full quantum mechanical values (QM) are compared with those obtained with the semiclassical approaches eqs. (15) and (18).