Variational Wigner-Kirkwood approach to relativistic mean field theory

Semiclassical methods are widely used to deal with pr erties of global character of different types of Fermi syste like atoms, nuclei, helium, or metallic clusters ~see Refs.@1, 2# for comprehensive reviews !. Concerning the nuclea ground-state energy, the simplest approach is the semiem ical mass formula based on the liquid drop model @3#. This approximation gives the smooth part of the energy and produces reasonably well the experimental values. The cess of the mass formula is due to the fact that the fluctua quantal correction to the nuclear energy ~shell effects! is small as compared with the average part. This allows a turbative treatment of the shell effects that is justified from theoretical point of view by the Strutinsky energy theore @4#. Actually, to obtain the correct semiclassical averaged ergy one should solve the problem using the Strutin smoothing. However, this is in general more difficult handle than the quantal problem itself if realistic nuclear tentials are used@5#. The search for alternative methods therefore an interesting and yet open problem. One poss ity is to use density functional theory ~DFT! techniques, such as the Thomas-Fermi ~TF! method and its extensions ~ETF!. These methods can be based on the Wigner-Kirkwood ~WK! \ expansion of the density matrix @2,6#. The WK theory provides an expansion of the particle and kinetic energy d sities in gradients of the single-particle potential up to zero second, or fourth order in \. These\ corrections come from the fact that in the Hamiltonian the momentum operator d not commute with the potential. In DFT the original W expansion of the density r is inverted to recast the kineti energy density as a functional t@r# of the density and its gradients. If the potential part of the interaction is al known as a functional of r, as happens with Skyrme force then minimization of the DFT energy leads to a variation equation forr alone@2,6#. Recently, the variational content of the WK expansion the energy has been studied for a set of nonrelativistic mions submitted to external~Woods-Saxon ! or selfconsistent~Skyrme! one-body potentials@7–9#. It has been shown that the variational solution for the particle dens that minimizes the semiclassical energy at each order in \ expansion is just the WK expansion of r at the same orde


I. INTRODUCTION
Semiclassical methods are widely used to deal with properties of global character of different types of Fermi systems like atoms, nuclei, helium, or metallic clusters ͑see Refs.͓1, 2͔ for comprehensive reviews͒.Concerning the nuclear ground-state energy, the simplest approach is the semiempirical mass formula based on the liquid drop model ͓3͔.This approximation gives the smooth part of the energy and reproduces reasonably well the experimental values.The success of the mass formula is due to the fact that the fluctuating quantal correction to the nuclear energy ͑shell effects͒ is small as compared with the average part.This allows a perturbative treatment of the shell effects that is justified from a theoretical point of view by the Strutinsky energy theorem ͓4͔.
Actually, to obtain the correct semiclassical averaged energy one should solve the problem using the Strutinsky smoothing.However, this is in general more difficult to handle than the quantal problem itself if realistic nuclear potentials are used ͓5͔.The search for alternative methods is therefore an interesting and yet open problem.One possibility is to use density functional theory ͑DFT͒ techniques, such as the Thomas-Fermi ͑TF͒ method and its extensions ͑ETF͒.These methods can be based on the Wigner-Kirkwood ͑WK͒ ប expansion of the density matrix ͓2,6͔.The WK theory provides an expansion of the particle and kinetic energy densities in gradients of the single-particle potential up to zeroth, second, or fourth order in ប.These ប corrections come from the fact that in the Hamiltonian the momentum operator does not commute with the potential.In DFT the original WK expansion of the density is inverted to recast the kinetic energy density as a functional ͓͔ of the density and its gradients.If the potential part of the interaction is also known as a functional of , as happens with Skyrme forces, then minimization of the DFT energy leads to a variational equation for alone ͓2,6͔.
Recently, the variational content of the WK expansion of the energy has been studied for a set of nonrelativistic fermions submitted to external ͑Woods-Saxon͒ or selfconsistent ͑Skyrme͒ one-body potentials ͓7-9͔.It has been shown that the variational solution for the particle density that minimizes the semiclassical energy at each order in the ប expansion is just the WK expansion of at the same order in ប.The method for solving this variational problem was presented in Ref. ͓7͔ and called variational Wigner-Kirkwood ͑VWK͒ theory.
The VWK theory has been mainly applied to the calculation of the surface energy coefficient in semi-infinite nuclear matter ͑SINM͒.It has been found to reproduce the quantal value nicely ͓7-9͔.From a comparison of VWK and DFT calculations ͓2,9,10͔, one can see that the agreement with the quantal results is worse in DFT ͑even at order ប 4 ͒ than in VWK theory ͓7-9͔.
In recent years the investigation of nuclear systems by means of a relativistic approach has attracted a growing interest.Especially successful has been the phenomenological relativistic mean field theory ͑RMFT͒ ͓11͔.The semiclassical approach to the relativistic theory has been recently worked out and ប corrections to the earlier established relativistic TF model have been derived ͓12-14͔, both at the WK and at the relativistic DFT ͑RDFT͒ level.It is the purpose of this paper to formulate, to second order in ប, the relativistic variational Wigner-Kirkwood ͑RVWK͒ theory, i.e., the VWK approach to RMFT.We will do that based on the nonlinearmodel ͓11,15͔.We do not attempt here to investigate the quality of the relativistic model.Rather, we want to compare the new semiclassical theory with the usual RDFT, adopting the quantal results as a standard, and try to find a good alternative to the Strutinsky averaging procedure.As an application of the theory we will calculate the surface energy coefficient in SINM.

II. RELATIVISTIC VARIATIONAL WIGNER-KIRKWOOD THEORY
Our starting point in setting the RVWK theory is the constrained energy E c of a finite nucleus: where E tot is the total energy density and the chemical potential is the Lagrange multiplier that ensures the correct par-ticle number A. Assuming for simplicity the nucleus to be symmetric and uncharged, in the relativistic nonlinearmodel E tot reads1 ͓11,15͔ In this equation V and are the vector and scalar fields, respectively, and the energy density E stands for ϩmϪg s s , with the relativistic kinetic energy density and s the scalar density.Quantally, we have ϭ͚ ␣ ␣ † ␣ , s ϭ ͚ ␣ ␣ † ␤ ␣ , and ϭ͚ ␣ ␣ † (Ϫi␣•"ϩ␤mϪm) ␣ on a single-particle basis ␣ .The free couplings and meson masses of the relativistic energy functional are usually fixed by adjusting them to reproduce nuclear matter properties or experimental data on finite nuclei within the Hartree approximation.
To second order in the WK expansion one has The well-known relativistic TF expressions are with the definition of a local Fermi momentum a Dirac effective mass m*ϭmϪg s , and ⑀ F ϭͱk F 2 ϩm* 2 .The second-order WK corrections 2 , s,2 , and E 2 contain squared gradients and Laplacians of the fields V and .Their expressions are considerably lengthier and will not be reproduced here as they can be found in Refs.͓12-14͔.On account of Eqs.͑3͒-͑7͒, the WK constrained energy density E tot Ϫ becomes a functional of the vector and scalar fields only.This is in accordance with the fact that in the relativistic framework there exist two different densities, namely, the baryon and the scalar densities, in contrast to the nonrelativistic situation where the basic quantity is the ground-state density, related to the one-body potential through the Fermi momentum.
With Eq. ͑3͒, however, the constrained energy ͑1͒ does not represent yet a complete separation into ប 0 and ប 2 parts.The vector and scalar fields and the chemical potential also have to be split into zeroth-and second-order terms: VϭV 0 ϩV 2 , ϭ 0 ϩ 2 , and ϭ 0 ϩ 2 .For instance, In what follows, the WK functionals 0 , 2 ,E 0 ,E 2 , etc., are to be evaluated at the V 0 , 0 , and 0 values.Introducing ϭ 0 ϩm, the separation of E c into pure TF (ប 0 ) and second-order (ប 2 ) contributions will then read

͑11͒
The gist of the derivation of the variational equations from the constrained energy E c ϭE c (0) ϩE c (2) ϩO(ប 4 ) in the RVWK theory is the fact that the minimization must be performed for each order in ប separately ͓7͔ ͑i.e., ␦E c With the aid of Eqs.͑12͒, at lowest order one gets from Eq. ͑9͒ the usual Euler-Lagrange equations for the TF fields V 0 and 0 ͓11,15͔: The Fermi momentum at the TF level is easily found from Eq. ͑7͒: It should be noticed that the density is already normalized at the TF level ͓7͔ ͑see the Appendix for details͒: Aϭ ͵ dr 0 .

͑15͒
Consequently, the full ប 2 contribution to A must vanish: ͑see below for the correction k 2 to the Fermi momentum͒.
The second-order Euler-Lagrange equations are derived from Eq. ͑11͒ and, using Eqs.͑12͒, turn out to be In practice, the procedure to solve these variational equations for a finite nucleus with A particles is the following.First, one solves the TF equations ͑13͒-͑15͒ in the usual manner ͓11͔ to find the fields V 0 and 0 and the chemical potential 0 .Inserting in Eqs.͑17͒ and ͑18͒ the TF solutions, one gets two linear differential equations to compute the so-far unknown corrections V 2 and 2 .These quantities will depend on 2 , which in turn has to be determined so that the normalization condition ͑16͒ is satisfied.Equations ͑16͒-͑18͒ will be iterated until consistency is reached.We will show below, however, that to calculate the correction of order ប 2 to the energy it is not necessary to obtain V 2 , 2 , and 2 ; only the TF solutions V 0 , 0 , and 0 are needed.
Once the solutions for the fields and the chemical potential are known, it is immediate to obtain the Fermi momentum to second order by direct expansion of Eq. ͑7͒: In principle, due to divergences at the classical turning point, in RVWK theory the densities and potentials beyond the TF order must be considered as distributions.However, they are very efficient to compute expectation values by integrals over the space ͓6,12,16͔.
It is worthwhile noting that the same variational equations ͑13͒, ͑14͒ and ͑17͒, ͑18͒ can be obtained if we start by directly varying Eq. ͑1͒ with respect to V and , now employing expressions ͑4͒ for and ͑6͒ for E, and only afterwards do we perform in the new Euler-Lagrange equations the splitting of V, , and into their ប 0 and ប 2 parts.Thus, expansion and variation can be interchanged.For brevity, we prove this statement for the vector field equation only.The variation with respect to V of Eq. ͑1͒ before expansion results in with given by Eq. ͑4͒.Now one makes VϭV 0 ϩV 2 and It is simple to see that Using these results and separating Eq. ͑20͒ into each order allows one to recover the previous variational equations ͑13͒ and ͑17͒.
On the other hand, one can follow an alternative route that consists in working with k F as if it were an independent variable of V and , i.e., without replacing it by Eq. ͑7͒.In this case, after expanding to order ប 2 (k F ϭk 0 ϩk 2 ), minimization with respect to k 0 (␦E c (0) /␦k 0 ϭ␦E c (2) /␦k 0 ϭ0) yields two variational equations which are just the zeroth-and second-order contributions to k F given by Eq. ͑19͒.Therefore, at each order in ប, the relation ͑7͒ between k F and the fields V and is the variational Euler-Lagrange equation for k F , a basic test of consistency.It also can be shown that when k F is kept as an independent variable, the variational equations for the fields reduce to Eqs. ͑13͒, ͑14͒ and ͑17͒, ͑18͒.
On the basis of Eqs.͑13͒ and ͑14͒ it is easy to verify that several contributions to Eq. ͑11͒ cancel out so that, taking into account the condition ͑15͒, the correction of order ប 2 to the total energy EϭE c ϩA finally is

͑23͒
Since the quantities E 2 and 2 are to be evaluated using the V 0 and 0 values, we arrive at the remarkable result that to calculate the energy to second order only the lowest-order TF equations ͑13͒ and ͑14͒ need to be solved.In this sense we have a perturbational approach, as the calculation of the energy to a given order in ប 2 requires knowledge of the solution only to the next lower order.A similar idea of including ប 2 corrections perturbatively, but in a rather more heuristic way, has been carried out in Ref. ͓17͔.In particular, this calculation differs from the present method in that the chemical potential is not split into ប 0 and ប 2 parts.
Though we have derived the VWK equations for the relativistic energy functional ͑2͒, it is easy to realize that the method is more general and that it can be applied to more realistic functionals ͑e.g., with an isovector meson, other types of scalar couplings, or self-interactions of the vector field͒.To conclude this section, we would like to mention again that the semiclassical calculation provides only the average part of the quantal energy.However, a semiclassical method can be useful to replace some complicated full quantal calculations ͓6͔.In such cases the remaining shell correc-tion should be added perturbatively to the semiclassical quantity, as indicated by the Strutinsky procedure or the expectation value method ͓2,5,18͔.

III. COMPARISON WITH RELATIVISTIC DENSITY FUNCTIONAL THEORY
The normalization condition ͑15͒ brings us to the discussion of an important difference between the RVWK theory we have just introduced and the standard RDFT ͓12,19,20͔.While in RVWK theory the result of the integral ͑16͒ vanishes, RDFT is more restrictive and reinforces this condition by imposing that the integrand vanishes locally.Then, in RDFT one has that, at each point, We write k ˜0 and k ˜2 to distinguish them from k 0 and k 2 of the RVWK theory.Owing to the ansatz ͑24͒, in RDFT the functional 0 (k ˜0) equals the exact density , i.e., ϭ2k ˜0 3 /3 2 .Notice that at the TF level, RVWK theory and RDFT are equivalent.
Next we briefly recall the derivation of the RDFT variational equations.In RDFT the scalar and vector fields and the chemical potential are not explicitly split into ប 0 and ប 2 parts ͓12,19,20͔; only k F is. Thus, expanding Eq. ͑2͒ into k ˜0 and k ˜2 and utilizing Eq. ͑24͒, the RDFT constrained energy to second order is where and ˜s,2 indicates that the gradients of V have been inverted as in E ˜2 and ˜2 .Note that having varied the energy ͑25͒ as a whole, and not independently for each order in ប, the final RDFT solution of Eqs.͑26͒-͑28͒ mixes different powers of ប.Nevertheless, the RDFT functionals are free from divergence problems at the classical turning point and generally provide a good description of the local density profile ͓2,7,12͔.

IV. NUMERICAL APPLICATION
To exemplify all the above on a concrete case, we shall present numerical calculations of the surface energy coefficient in SINM for several parameter sets of the relativistic interaction.The semi-infinite system corresponds to a onedimensional geometry where half the space is filled with nuclear matter at saturation and the other half is empty.The particle density then varies only along one axis, e.g., the z axis, and develops a surface around zϭ0.Following Ref. ͓3͔, the surface energy coefficient E s in SINM is written as where r ϱ and a v refer to the radius and energy per particle in saturated nuclear matter.In the self-consistent problem for SINM, a v equals the chemical potential , owing to the Hugenholtz-Van Hove theorem.Consequently, the surface energy is stationary with respect to variations of the density.
Our previous formulation can be fully applied to the onedimensional semi-infinite geometry ͓compare Eqs.͑1͒ and ͑30͔͒ with the simplification that is fixed and thus 2 ϭ0.First, we will discuss SINM results for the linearmodel (bϭcϭ0).In spite of its simplicity, the linear model allows one to investigate the incidence of the gradient corrections more easily.Second, we will consider the more general nonlinearmodel.We are mainly interested in extracting the average part of the energy associated withinteractions whose parameters have been obtained from a mean field calculation.Within the framework of the RMFT, and taking the quantal Hartree calculations as a standard, the authors of Refs.͓12, 18-21͔ studied extensively the quality of the relativistic TF and RDFT ͑to order ប 2 ͒ approximations, and how the results depend on the parametrization of the effective interaction.It was seen that one must analyze more than just one single parameter set to draw conclusions about the quality of the semiclassical approximations, for it much depends on some properties of the interaction.For our application to the surface of SINM, the mass of the scalar meson, m s , and the effective mass at saturation, m ϱ */m, have the major influence.The discussion that follows is not significantly altered by changing the remaining saturation properties if they lie within ordinary values ͓19-21͔.
Table I collects the results for the surface energy coefficient E s calculated in the linearmodel.Table II  The values in Table I correspond to the discussed semiclassical approaches ͑TF, RVWK, and RDFT͒ and the fully quantal Hartree ͑H͒ calculation.͑See Ref. ͓22͔ for details on the quantal treatment of relativistic SINM.͒In any WK calculation beyond the TF order one is faced with dealing with divergences at the turning point.The main difficulty is to treat them in such a way that the principal part of the energy can be extracted.In practice, we encountered only one divergent term in the WK functional E 2 .To get rid of the divergence we added and subtracted the analytical asymptotic integrand, so that after having isolated the infinity we rejected it.Other possibilities are the method of integrals in the complex plane ͓16͔ or low-temperature expansions ͓2͔.
The saturation properties of infinite nuclear matter are governed by the meson coupling-to-mass ratios g s 2 /m s 2 and g v 2 /m v 2 and by the nonlinear couplings b and c ͓11,15͔.On the contrary, the nuclear surface properties extracted from SINM depend on the meson coupling constants and masses separately.The mass of the vector meson, m v , is given its physical value ͑783 MeV͒.The mass of the scalar meson, m s , should lie somewhere between 400 and 700 MeV, since the nonexperimental particle is interpreted as simulating two-pion exchange contributions.For our purposes, it will be sufficient to look at the region 400 MeVрm s р550 MeV.The scalar mass sets the range of the scalar interaction and, therefore, there is a strong correlation of E s and t with m s .A larger m s determines a shorter range of the attractive potential, leading to a steeper surface and to smaller E s and t, as seen from Tables I and II.
The TF surface energy coefficients in Table I overestimate the quantal ones from ϳ4% for m s ϭ400 MeV to ϳ13% for m s ϭ550 MeV.When the ប 2 gradient corrections are taken into account, one finds that the surface energies calculated in the RVWK approach are larger than the H results, whereas the RDFT energies are smaller than in H calculations.In both cases E s is brought closer to the H value than in the TF calculation.In RVWK theory the deviations lie between 0.5% (m s ϭ400 MeV) and 2.5% (m s ϭ550 MeV), while in RDFT they range from 2% (m s ϭ400 MeV) to 9% (m s ϭ550 MeV).From Table I  show an upward tendency with the scalar mass m s , and that E s TF ϪE s DFT is systematically larger than E s TF ϪE s VWK .These trends can be qualitatively understood looking at the values of the surface thickness t in Table II.The inhomogeneity corrections of RVWK theory and RDFT concentrate at the nuclear surface, where the gradients are more important.A flatter surface ͑small m s , large t͒ results in smaller corrections.In RVWK theory one calculates the gradients with the TF density distributions that have a larger thickness than the RDFT profiles.Therefore, one expects the RVWK modifications to the TF energy to be smaller.
For relativistic harmonic oscillator scalar and vector fields, it has been numerically shown ͓12͔ that the energies calculated using the Strutinsky average and the WK approach to second order almost coincide, a well-known fact in the nonrelativistic frame ͓5͔.It is thus reasonable to identify approximately the difference between H and RVWK calculations in the self-consistent problem with the quantal effects ͑in SINM, Friedel oscillations, and the fluctuating part of the spin-orbit force ͓21,22͔͒.Even though the quantal surface energy coefficient is acceptably reproduced by RDFT in general, the difference with H calculations is larger than in RVWK theory and displays a stronger dependence on the particular value of m s .Altogether RVWK theory appears as more reliable to estimate the quantal effects.This feature, also found in the nonrelativistic context ͓7͔, stems from the following reasons.First, RVWK theory properly sorts out the different orders in ប.Second, the restrictive local condition ͑24͒ for normalization within RDFT is replaced by the more logical global condition ͑16͒ in RVWK theory.
The surface energy is also strongly correlated with the value of the effective mass in nuclear matter m ϱ */m ͓9,19,20͔.To analyze this fact we have considered the nonlinear parameter sets of Ref. ͓20͔.They have a v ϭϪ15.75MeV, ϱ ϭ0.16 fm Ϫ3 , and Kϭ200 MeV, with 0.55рm ϱ */mр0.80 and 400 MeVрm s р550 MeV which covers the range of commonly accepted values for m ϱ */m and m s .Figure 1 illustrates the dependence of the difference between the TF and H surface energy coefficients on m ϱ */m and m s .The discrepancies between the TF and H results exhibit a nearly linear behavior with m ϱ */m.For small values of the effective mass the TF surface energy coefficients are larger than the H ones.They practically agree with the H results for m ϱ */mӍ0.65,and become smaller than the H re- sults for larger m ϱ */m.These trends have been found in a similar fashion for the total energies of finite nuclei in the model ͓18͔ and in nonrelativistic calculations ͓9͔.
Figure 2 displays the difference in the surface energy coefficient between the semiclassical approaches to secondorder and H calculations.Again, for all the analyzed parameter sets, we observe that the RVWK energies are larger than the H ones while the RDFT energies are smaller.Also, the deviations to the H values are an increasing function of the scalar mass m s .Their dependence on m s and on m ϱ */m is weaker in RVWK theory than in RDFT ͑and, in both approaches, much weaker than in the TF method͒.Finally, one realizes that the semiclassical surface energies can show large discrepancies with the quantal value for some parameter sets ͑e.g., m s ϭ550 MeV and m ϱ */mϭ0.55 and 0.60͒.This is due to the fact that in such cases the thickness of the particle density and local effective mass is very small.Their gradients are therefore very large, causing the semiclassical expansions to break down.

V. SUMMARY
In this paper we have developed the relativistic variational Wigner-Kirkwood theory, extending first work in the nonrelativistic domain.This generalization is not trivial because of the presence of two different fields.The RVWK theory has been based on a strict expansion of the energy in powers of ប, together with a global normalization to the particle number.Self-consistency enters at the TF level, whose solution is the input to calculate the higher-order corrections.In obtaining the variational equations we have shown that the steps of variation and expansion can be interchanged.We also have discussed the equivalence of working with the scalar and vector fields as the fundamental variables or with the Fermi momentum as an additional variable.The new theory has been compared with the RDFT.
Semi-infinite nuclear matter calculations in the relativistic problem have shown that the average part of the quantal surface energy is acceptably estimated in both the RVWK and RDFT approaches.However, for a quantitative estimate of the quantal effects, it has been seen that the RVWK theory is preferable.In addition, it must be considered an advantage of RVWK theory that its quality is less dependent on the properties of the effective interaction than in the TF and RDFT approximations.
According to textbook theory of Lagrange multipliers, to locate the extrema of Eq. ͑A1͒ with the restriction ͑A2͒, one has to seek the critical points of the auxiliary functional K.

͑A10͒
Equation ͑A8͒ is just ␦K 0 ␦ 0 ϭ0.͑A11͒ Inserting Eq. ͑A8͒ into Eq.͑A7͒ shows that Eq. ͑A7͒ is equivalent to ␦K 2 ␦ 0 ϭ0.͑A12͒ Equation ͑A10͒ tells us that the lowest order already fulfills the restriction condition ͑A2͒, and in view of Eq. ͑A9͒ the total second-order contribution to A vanishes: Therefore, one can paraphrase the variational problem under consideration as follows.Minimization must be performed for each order in the expansion parameter separately, Eqs.͑A11͒ and ͑A12͒, and the constraint must be satisfied by the zeroth-order solution.Once 0 and 0 are known by solving Eqs.͑A10͒ and ͑A11͒, they are inserted into Eqs.͑A12͒ and ͑A13͒ to calculate the second-order corrections 2 and 2 .Expanding Eq. ͑A1͒,

͑A16͒
This result shows that the extremum of Eq. ͑A1͒ can be computed to second order from the knowledge of 0 and 0 only.We then observe that the whole procedure is consistent with the spirit of perturbation theory, since the lowest-order solution serves as the input to calculate the higher-order corrections.It furthermore guarantees that different powers of the expansion parameter do not mix at each order of the expansion.
shows the corresponding surface thickness t of the density profile ͑standard 90% to 10% falloff distance͒.The nuclear matter properties of the linear model are as follows ͓20͔: volume energy a v ϭϪ15.75MeV, density ϱ ϭ0.193 fm Ϫ3 , incompressibility Kϭ546 MeV, and effective mass m ϱ */mϭ0.56.
one can check that the differences E s TF ϪE s DFT and E s TF ϪE s

TABLE I .
Linearmodel.Surface energy coefficient in the quantal Hartree approach, E s H , and difference of the relativistic TF, DFT, and VKW calculations to E s H .All quantities are in MeV.