J Dobaczewski

B G Carlsson

M Kortelainen

The Negele-Vautherin density matrix expansion applied to the Gogny force

Keywords: numbers: 21.60.Jz, 21.10.Dr, 21.30.-x

We use the Negele-Vautherin density matrix expansion to derive a quasi-local density functional for the description of systems of fermions interacting with shortranged interactions composed of arbitrary finite-range central, spin-orbit, and tensor components. Terms that are absent in the original Negele-Vautherin approach owing to the angle averaging of the density matrix are fixed by employing a gauge invariance condition. We obtain the Kohn-Sham interaction energies in all spin-isospin channels, including the exchange terms, expressed as functions of the local densities and their derivatives up to second (next to leading) order. We illustrate the method by determining the coupling constants of the Skyrme functional or Skyrme force that correspond to the finite-range Gogny central force. The resulting self-consistent solutions reproduce the Gogny-force binding energies and radii within the precision of 1-2%.

Introduction

The search for a universal energy density functional (EDF) [1], that would be able to provide a spectroscopic-quality [2] description of atomic nuclei, is at the focus of the present-day studies in nuclear structure. Recently, several advanced phenomenological analyses were aimed at improving the standard relativistic [3] or nonrelativistic local Confidential: not for distribution. Submitted to IOP Publishing for peer review 22 February 2010 [4,2,5,6,7] or nonlocal functionals [8]. A significant ongoing effort is also devoted to microscopic derivations of the functionals (see, for example, Refs. [9,10,11].) One of the central points of the current EDF studies is the question: to what extent can the finite-range effective interactions be approximated as quasi-local density functionals?

The framework to build the quasi-local theory was set up in the seminal paper by Negele and Vautherin (NV) [12], which introduced the so-called density matrix expansion (DME) method. Later, other methods to achieve the same goal, like the semiclassical expansion [13,14] were also proposed and studied. The original NV expansion, for the scalar density, allowed for treating only the even-order terms in relative coordinates, and thus was applicable only to even-power gradient densities.

In the present study, we revisit the NV expansion by adding the odd-power gradient densities through the local gauge-invariance condition. There are two other important differences of the present approach with respect to the NV original, namely, (i) we treat the spin (vector) densities analogously to the scalar densities and (ii) we apply a different DME for the direct terms. These differences are motivated by the effective-theory [START_REF] Lepage | Lectures given at the VIII Jorge Andre Swieca Summer School[END_REF] interpretation of the DME advocated in the present study. The ultimate test of these ideas can only be obtained by analysing microscopic properties of nuclear densities [11].

We present a complete set of expressions in all spin-isospin channels, applicable to arbitrary finite-range central, spin-orbit, and tensor interactions. Our study is focused on applying the NV expansion to the Gogny interaction [START_REF] Dechargé | [END_REF]17]. This allows us to look at the correct scale of the interaction range that properly characterizes nuclear lowenergy phenomena. By performing expansion up to second order (or next-to-leading order, NLO), one obtains the local Skyrme functional [18,19]. In this way, we establish a firm link between two different, local and nonlocal, very successful functionals.

The paper is organized as follows. In Section 2, we introduce reformulation of the NV expansion in terms of particles without spin. Discussion of this example allows us to present details of the approach without complications that otherwise could have obscured the main ideas. In Sections 3 and 4, we present results with spin and isospin degrees of freedom reintroduced, with Section 4 containing applications of the formalism to the Gogny force. Conclusions are given in Section 5, and Appendix A contains the discussion of spin and isospin polarized nuclear matter. Preliminary results of the present study were published in Ref. [START_REF] Dobaczewski | Trends in Field Theory Research[END_REF].

Local energy density for spinless particles of one kind

In this section, we consider the simplest (and academic) case of fermions with no spin and no isospin. First we recall that for an arbitrary non-local finite-range interaction V (r 1 , r 2 ; r 1 , r 2 ), the Kohn-Sham interaction energy [START_REF] Kohn | [END_REF] has the form

E int = 1 2 d 3 r 1 d 3 r 2 d 3 r 1 d 3 r 2 V (
r 1 , r 2 ; r 1 , r 2 ) × (ρ(r 1 , r 1 )ρ(r 2 , r 2 ) -ρ(r 2 , r 1 )ρ(r 1 , r 2 )), (1) whereas for a local interaction, V (r 1 , r 2 ; r 1 , r 2 ) = δ(r 1 -r 1 )δ(r 2 -r 2 )V (r 1 , r 2 ), (2) the interaction energy reduces to:

E int = 1 2 d 3 r 1 d 3 r 2 V (r 1 , r 2 ) ρ(r 1 )ρ(r 2 ) -ρ(r 2 , r 1 )ρ(r 1 , r 2 ) , (3)
where ρ(r 1 ) ≡ ρ(r 1 , r 1 ) and ρ(r 2 ) ≡ ρ(r 2 , r 2 ) are local densities. As is well known, the first term in Eq. ( 3) (the direct term) depends only on local densities, whereas the second one (the exchange term) depends on the modulus squared of the non-local density. This markedly different structure of the two terms requires separate treatment, as discussed in the following two subsections.

Direct interaction energy

In nuclei, the range of interaction is significantly smaller than the typical scale of the distance at which the local density varies. Therefore, we may expand the local densities ρ(r 1 ) and ρ(r 2 ) around their average position, and use this expansion to calculate the direct term in Eq. ( 3). Denoting the standard total (R) and relative (r) coordinates and derivatives as

R = 1 2 (r 1 + r 2 ), r = r 1 -r 2 , (4)
∇ = ∂ ∂R = ∂ ∂r 1 + ∂ ∂r 2 , ∂ = ∂ ∂r = 1 2 ∂ ∂r 1 - ∂ ∂r 2 , (5)
we have the expansion of local densities to second order,

ρ(r 1 ) = ρ(R + 1 2 r) = ρ(R) + 1 2 r a ∇ a ρ(R) + 1 8 r a r b ∇ a ∇ b ρ(R) + . . . , (6)
ρ(r 2 ) = ρ(R -1 2 r) = ρ(R) -1 2 r a ∇ a ρ(R) + 1 8 r a r b ∇ a ∇ b ρ(R) + . . . , (7)
and hence

ρ(r 1 )ρ(r 2 ) = ρ 2 (R)+ 1 4 r a r b ρ(R)∇ a ∇ b ρ(R)-[∇ a ρ(R)][∇ b ρ(R)] +. . . , (8
) where we implicitly assumed the summation over the repeated Cartesian indices a and b.

Assuming that the local potential V (r 1 , r 2 ) depends only on the distance between the interacting particles, V (r 1 , r 2 )=V (|r 1 -r 2 |)=V (r), the direct interaction energy is given by the integral of a local energy density H int dir (R),

E int dir = d 3 R H int dir (R), (9)
where up to second order,

H int dir (R) = 1 2 V 0 ρ 2 + 1 12 V 2 ρ∆ρ -(∇ρ) 2 + . . . , (10)
where the coupling constants, V 0 and V 2 , are given by the lowest two moments of the interaction,

V n = d 3 r r n V (r) = 4π dr r n+2 V (r). ( 11
)
After integrating by parts, the two second-order terms in Eq. ( 10) are identical, so we can equally well use:

H int dir (R) = 1 2 V 0 ρ 2 + 1 6 V 2 ρ∆ρ + . . . , (12)
We see that the separation of scales between the range of interaction and the rate of change of the local density leads to a dramatic collapse of information that is transferred from the interaction potential to the interaction energy. Namely, the two constants, V 0 and V 2 , completely characterize the interaction in the direct term, and the detailed form of the potential V (r) becomes irrelevant. Moreover, it can be easily checked that in this approximation, the direct energy density is exactly equal to that corresponding to the contact force corrected by the second-order gradient pseudopotential,

Ṽ (r) = V 0 δ(r) -2 3 V 2 ∂ • δ(r)∂. ( 13
)

Exchange interaction energy

In the exchange term of Eq. ( 3), the non-zero range of the interaction probes the nonlocal space dependence of the density matrix. For short-range interactions, one can expand ρ(R, r) to second order with respect to the variable r, which gives

ρ(r 1 , r 2 ) = ρ(R, r) = ρ(R) + r a ∂ a ρ(R, r) + 1 2 r a r b ∂ a ∂ b ρ(R, r) + . . . = ρ(R) + ir a j a (R) + 1 2 r a r b 1 4 ∇ a ∇ b ρ(R) -τ ab (R) + . . . , (14)
where derivatives ∂ i are always calculated at r a =0, and therefore, the result can be expressed in terms of the standard current and kinetic densities [22]:

j a (R) = 1 i ∂ a ρ(R, r) r=0 , τ ab (R) = ∇ (1) a ∇ (2) b ρ(r 1 , r 2 ) r 1 =r 2 . ( 15
)
This parabolic approximation does not ensure that ρ(r 1 , r 2 )-→0 for large r = |r| = |r 1 -r 2 |. In the spirit of the DME [12], one can improve it by introducing three functions of r, π 0 (r), π 1 (r), and π 2 (r) [START_REF] Dobaczewski | Trends in Field Theory Research[END_REF] that vanish at large r, i.e., we define the quasi-local approximation of the density matrix by:

ρ(r 1 , r 2 ) = π 0 (r)ρ(R) + π 1 (r)r a ∂ a ρ(R, r) + 1 2 π 2 (r)r a r b ∂ a ∂ b ρ(R, r) + . . . = π 0 (r)ρ(R) + iπ 1 (r)r a j a ρ(R) + 1 2 π 2 (r)r a r b 1 4 ∇ a ∇ b ρ(R) -τ ab (R) + . . . . (16)
Such a postulate has to be compatible with the Taylor expansion of Eq. ( 14), which requires that

π 0 (0) = π 1 (0) = π 2 (0) = 1 and π 0 (0) = π 1 (0) = π 0 (0) = 0. ( 17
)
Of course, for π 0 (r)=π 1 (r)=π 2 (r)=1, one reverts to the parabolic approximation (14). The product of nonlocal densities in the exchange integral of Eq. (3) to second order reads (18) where we have introduced the supplementary condition:

ρ(r 1 , r 2 )ρ(r 2 , r 1 ) = π 2 0 (r)ρ 2 (R) + π 0 (r)π 2 (r)r a r b ρ(R)∂ a ∂ b ρ(R, r) -[∂ a ρ(R, r)][∂ b ρ(R, r)] + . . . = π 2 0 (r)ρ 2 (R) + π 0 (r)π 2 (r)r a r b 1 4 ρ(R)∇ a ∇ b ρ(R) -ρ(R)τ ab (R) + j a (R)j b (R) + . . . ,
π 2 1 (r) = π 0 (r)π 2 (r). ( 19
)
This condition ensures that the quasi-local approximation of Eq. ( 16) is compatible with the local gauge invariance [23]. Indeed, the left-hand side of Eq. ( 18) is manifestly invariant with respect to the local gauge transformation,

ρ (r 1 , r 2 ) = e iφ(r 1 )-iφ(r 2 ) ρ(r 1 , r 2 ), (20)
and only the difference of terms in the curly brackets in Eq. ( 18) is invariant with respect to the same transformation [22,19]. Functions π 0 (r), π 1 (r), and π 2 (r) also depend on the parameters defining the approximation [START_REF] Dechargé | [END_REF]. In particular, when the infinite matter is used to define functions π 0 (r), π 1 (r), and π 2 (r), like in the DME, they parametrically depend on the Fermi momentum k F . By associating the local density ρ(R) with k F , functions π 0 (r), π 1 (r), and π 2 (r) become dependent on ρ(R), and hence the damping of the density matrix in the non-local direction r can be different in different local points R. However, in order to keep the notation simple, we do not explicitly indicate this possible dependence on density.

Within the quasi-local approximation, one obtains the exchange interaction energy,

E int exc = d 3 R H int exc (R), (21)
where up to second order,

H int exc (R) = -1 2 V 00 π0 ρ 2 + 1 3 V 02 π2 1 4 ρ∆ρ -(ρτ -j 2 ) + . . . , (22)
and where τ = τ aa . The coupling constants, V 00 π0 and V 02 π2 , are given by the following moments of the interaction,

V ij πn = d 3 r r n π i (r)π j (r)V (r) = 4π dr r n+2 π i (r)π j (r)V (r). ( 23
)
Unlike the coupling constants defining the direct term (11), these in Eq. ( 23) have to be understood as the running coupling constants, which in the DME depend on the scale of the Fermi momentum k F or density ρ(R).

Again, the separation of scales between the range of interaction and the rate of change of the density matrix in the non-local direction results in the dependence of the local energy density on two coupling constants only, and not on the details of the interaction. For the parabolic approximation of Eq. ( 14), the coupling constants that define the direct and exchange energies are identical, i.e., V 00 π0 =V 0 and V 02 π2 =V 2 ; however, for the quasi-local approximation of Eq. ( 16) they are different. This important observation is discussed in Sec. 3.2 in more detail.

In nuclei, the separation of scales discussed above is not very well pronounced. The characteristic parameter, which defines these relative scales, is equal to k F a, where a stands for the range of the interaction. For example, within the Gogny interaction, which we analyze in detail in Sec. 4, there are two components with the ranges of a = 0.7 and 1.2 fm, whereas k F = 1.35 fm -1 , which gives values of k F a = 0.95 and 1.62 that are dangerously close to 1. Therefore, in the expansion of the local energy density (22), one cannot really count on the moments V ij πn decreasing with the increasing order n. Instead, as demonstrated by Negele and Vautherin [12], one can hope for tremendously improving the convergence by using at each order the proper counterterms, which make each order vanish in the infinite matter. Without repeating the original NV construction, here we only note that the net result consists of adding and subtracting in Eq. ( 16) the infinite-matter term, that is,

ρ(r 1 , r 2 ) = ν 0 (r)ρ(R) + iν 1 (r)r a j a ρ(R) + 1 2 ν 2 (r)r a r b 1 4 ∇ a ∇ b ρ(R) -τ ab (R) + 1 5 δ ab k 2 F ρ(R) + . . . , (24)
where we have defined functions ν i (r) such that:

ν 0 (r) = π 0 (r) -1 10 (k F r) 2 π 2 , ν 1 (r) = π 1 (r), ν 2 (r) = π 2 (r). ( 25
)
By neglecting the term quadratic in ν 2 , we can now use the approximation (24) to calculate the product of densities in Eq. ( 18), which gives

ρ(r 1 , r 2 )ρ(r 2 , r 1 ) = ν 2 0 (r)ρ 2 (R) + ν 0 (r)ν 2 (r)r a r b 1 4 ρ(R)∇ a ∇ b ρ(R) -ρ(R)τ ab (R) + 1 5 δ ab k 2 F ρ 2 (R) + j a (R)j b (R) + . . . , (26)
where again the gauge invariance requires that

ν 2 1 (r) = ν 0 (r)ν 2 (r). (27)
We note that the gauge-invariance conditions for the functions π i (r) and ν i (r), Eqs. ( 19) and (27), are compatible with one another only up to the term π 2 2 (r), which was shifted to higher orders. Finally, approximation (26) gives the energy density analogous to (22),

H int exc (R) = -1 2 V 00 ν0 ρ 2 + 1 3 V 02 ν2 1 4 ρ∆ρ -(ρτ -j 2 ) + 3 5 k 2 F ρ 2 + . . . , (28)
where the coupling constants, V 00 ν0 and V 02 ν2 , are given by the moments of the interaction calculated for functions ν i (r), namely,

V ij νn = d 3 r r n ν i (r)ν j (r)V (r) = 4π dr r n+2 ν i (r)ν j (r)V (r). ( 29
)
We see that the two sets of auxiliary functions, π i (r) and ν i (r), are suitable for discussing the approximate forms of the nonlocal density ρ(r 1 , r 2 ) and exchange energy density H int exc (R), respectively. Although for the nonlocal density they correspond to a simple reshuffling of terms, which gives relations (25) between π i (r) and ν i (r), for the exchange energy density they constitute entirely different approximations, given in Eqs. ( 22) and ( 28), with expansion (28) having a larger potential for faster convergence. Only this latter expansion is further discussed.

Determination of functions π i (r)

For a given nonlocal density ρ(r 1 , r 2 ), the auxiliary functions π i (r) or ν i (r), which define its quasi-local approximation, can be calculated as their best possible approximations in terms of local densities. However, the usefulness of the expansion relies on the assumption that generic forms of these functions can be estimated and then applied to all many-body systems of a given kind.

The standard Slater approximation [24,25], which is routinely used to evaluate the Coulomb exchange energy (cf. Refs. [26,27]), corresponds to

π 0 (r) = ν 0 (r) = 3j 1 (k F r) k F r and π 2 (r) = 0. ( 30
)
The NV expansion [12] gives a second-order estimate by making the momentum expansion around the Fermi momentum k F of an infinite system. This gives:

π 0 (r) = 6j 1 (k F r) + 21j 3 (k F r) 2k F r 1 - (k F r) 4 504 , (31)
π 2 (r) = 105j 3 (k F r) (k F r) 3 1 - (k F r) 2 18 + (k F r) 4 792 , (32)
ν 0 (r) = 3j 1 (k F r) k F r 1 - (k F r) 2 10 + (k F r) 4 280 , (33)
where j n (k F r) are the spherical Bessel functions.

The NV functions π 0 (r), π 2 (r), and ν 0 (r) are plotted in Fig. 1 with solid, dashed, and dotted lines, respectively. One can see that for large k F r, functions π 0 (r) and π 2 (r) have zeros close to one another and the same signs. Indeed, asymptotically both behave like cos(k F r)/k F r. Therefore, the gauge-invariance condition ( 19) can be satisfied almost everywhere. On the other hand, functions ν 0 (r) and ν 2 (r) = π 2 (r), have asymptotically opposite signs, and the corresponding gauge-invariance condition ( 27) can be almost nowhere satisfied. Nevertheless, as discussed in Sec. 2.2, what really matters are the moments of interaction ( 23) and (29), where functions π i (r) and ν i (r) are probed only within the range of the interaction, that is, up to k F r 1-2. Therefore, for the NV expansion, one can safely use approximations:

π 1 (r) = + |π 0 (r)π 2 (r)|, (34)
ν 1 (r) = + |ν 0 (r)ν 2 (r)|, (35)
which are valid up to the first zero of j 3 or j 1 , respectively, that is, up to k F r 7.0 and 4.5.

By the same token, we can replace in Eqs. ( 31)-(33) the Bessel functions by Gaussians having the same leading-order dependence on k F r, namely,

π 0 (r) = exp - (k F r) 4 504 , (36)
π 2 (r) = exp - (k F r) 2 18 - (k F r) 4 1100 , (37)
ν 0 (r) = exp - (k F r) 4 504 - (k F r) 2 10 exp - (k F r) 2 18 - (k F r) 4 1100 . ( 38
)
As seen in the bottom panel of Fig. 1, in this way, in the region of small k F r, one obtains a very good reproduction of the NV functions π i (r) and ν i (r).

Local energy density for particles with spin and isospin

Density matrix expansion with spin and isospin

For nucleons, the density matrix ρ(r 1 σ 1 τ 1 , r 2 σ 2 τ 2 ) depends not only on positions r 1 and r 2 but also on spin σ 1 , σ 2 = ±1 and isospin τ 1 , τ 2 = ±1 coordinates. Since the strong two-body interaction is assumed to be isospin and rotationally invariant, it is convenient to represent the standard density matrix ρ(r 1 σ 1 τ 1 , r 2 σ 2 τ 2 ) through nonlocal densities ρ µk (r 1 , r 2 ) as:

ρ(r 1 σ 1 τ 1 , r 2 σ 2 τ 2 ) = 1 4 µ=0,x,y,z 3 k=0 ρ µk (r 1 , r 2 ) σ 1 |σ µ |σ 2 τ 1 |τ k |τ 2 , (39)
where σ 0 (σ x,y,z ) and τ 0 (τ 1,2,3 ) are the unity (Pauli) matrices in the spin and isospin coordinates, respectively. For µ = 0 and µ = x, y, z, the densities are scalars and vectors, respectively, and for k = 0 and k = 1, 2, 3, they are isoscalars and isovectors, so altogether the density matrix is split into the four standard spin-isospin channels.

For the direct term, we can proceed as in Sec. 2.1, by making the Taylor expansions of local densities (at r 1 = r 2 ) in each spin-isospin channel; that is, similarly as in Eqs. ( 6) and ( 7), we have

ρ µk (R ± 1 2 r) = ρ µk (R) ± 1 2 r a ∇ a ρ µk (R) + 1 8 r a r b ∇ a ∇ b ρ µk (R) + . . . . (40)
For the exchange term, Sec. 2.2, the analogous Taylor expansions of nonlocal densities, similarly as in Eq. ( 14), read

ρ µk (R, ±r) = ρ µk (R) ± ir a j µak (R) + 1 2 r a r b 1 4 ∇ a ∇ b ρ µk (R) -τ µabk (R) + . . . , (41
) where the current (j µak (R)) and kinetic (τ µabk (R)) densities are defined in each channel as in Eqs. [START_REF] Lepage | Lectures given at the VIII Jorge Andre Swieca Summer School[END_REF], namely,

j µak (R) = 1 i ∂ a ρ µk (R, r) r=0 , τ µabk (R) = ∇ (1) a ∇ (2) b ρ µk (r 1 , r 2 ) r 1 =r 2 . (42)
The local density approximation of densities in all channels, analogous to Eqs. ( 16) and (24), is now postulated as

ρ µk (R, ±r) = π 0 (r)ρ µk (R) ± iπ 1 (r)r a j µak (R) + 1 2 π 2 (r)r a r b 1 4 ∇ a ∇ b ρ µk (R) -τ µabk (R) + . . . , (43)
and

ρ µk (R, ±r) = ν 0 (r)ρ µk (R) ± iν 1 (r)r a j µak ρ(R) + 1 2 ν 2 (r)r a r b 1 4 ∇ a ∇ b ρ µk (R) -τ µabk (R) + 1 5 δ ab k 2 F ρ µk (R) + . . . . (44
) At this point, we have assumed that functions π i (r) and ν i (r) are channelindependent, that is, that they are scalar-isoscalar functions. In Appendix A we discuss this point in more detail, and we show that the postulate of simply channel-dependent functions π i (r) and ν i (r) is incompatible with properties of infinite matter, whereas the proper treatment of the problem leads immediately to the channel mixing and the energy density, which is not invariant with respect to rotational and isospin symmetries. This question certainly requires further study, whereas at the moment, a consistent approach can only be obtained by assuming the scalar-isoscalar functions π i (r) and ν i (r).

We can now apply derivations presented in Secs. 2.1 and 2.2 to the general case of an arbitrary finite-range local nuclear interaction composed of the standard central, spin-orbit, and tensor terms:

V (r 1 , r 2 ) = W (r) + B(r)P σ -H(r)P τ -M(r)P σ P τ + P (r) + Q(r)P τ L • S + R(r) + S(r)P τ S 12 , (45)
where r=|r|=|r 1 -r 2 |, and

P σ = 1 2 (1 + σ 1 • σ 2 ), P τ = 1 2 (1 + τ 1 • τ 2 ), (46)
L = -ihr×∂, S = h 2 (σ 1 +σ 2 ), S 12 = 3 r 2 (σ 1 •r)(σ 2 •r)-σ 1 •σ 2 .( 47
)
After straightforward but lengthy calculations, one obtains the interaction energy in the form of a local integral, analogous to that for the Skyrme interaction [28,22,19],

E int = d 3 R 3 k=0 C ρ t ρ 2 k + C ∆ρ t ρ k ∆ρ k + C τ t ρ k τ k -j 2 k + C s t s 2 k + C ∆s t s k • ∆s k + C T t s k • T k -J abk J abk + C F t s k • F k -1 2 J aak J bbk -1 2 J abk J bak + C ∇s t ∇ • s k 2 + C ∇J t ρ k ∇ • J k + s k • (∇ × j k ) , (48)
where ρ k ≡ ρ 0k , τ k ≡ τ 0bbk , j ak ≡ j 0ak , s ak ≡ ρ ak , T ak ≡ τ abbk , F ak ≡ 1 2 (τ babk + τ bbak ), J abk ≡ j abk , and J ak ≡ abc j cbk are the standard local densities. The isoscalar (t = 0) and isovector (t = 1) coupling constants C t correspond to k = 0 and k = 1, 2, 3, respectively.

The coupling constants of the local energy density (48) are related to moments of the interaction in the following way:

8        C ρ 0 C ρ 1 C s 0 C s 1        =        4 2 -2 -1 0 0 -2 -1 0 2 0 -1 0 0 0 -1               W 0 + M 00 ν0 + 1 5 M 02 ν2 k 2 F B 0 + H 00 ν0 + 1 5 H 02 ν2 k 2 F H 0 + B 00 ν0 + 1 5 B 02 ν2 k 2 F M 0 + W 00 ν0 + 1 5 W 02 ν2 k 2 F        , (49) 96 
                              C ∆ρ 0 C ∆ρ 1 C τ 0 C τ 1 C ∆s 0 C ∆s 1 C T 0 C T 1 C F 0 C F 1 C ∇s 0 C ∇s 1                                =                                8 -1 4 -2 -4 2 -2 4
0 0 0 0 0 -1 0 -2 -4 0 -2 0 0 0 0 0 0 4 0 8 0 -8 0 -16 0 0 0 0 0 4 0 8 0 0 0 0 0 0 0 0 0 -1 4 0 0 2 -2 0 -4 1 -2 2 0 -1 0 0 0 0 -2 0 0 1-2 0 0 4 0 0 0 -8 0 0 0 -4 0 -8 0 4 0 0 0 0 0 0 0 -4 0 0 0 0 0 0 0 0 0 0 0 12 0 24 0 0 0 0 0 0 0 0 0 12 0 0 0 0 0 0 0 0 0 0 -12 3 -6 6 0 0 0 0 0 0 0 0 0 3 -6 0

  C ∇J 0 C ∇J 1   =   2 1 0 1     P 2 + Q 01 ν2 Q 2 + P 01 ν2   . ( 51
)
All the coupling constants of the local energy density (48) depend linearly on the following moments of potentials:

X n = d 3 r r n X(r) = 4π d r r n+2 X(r), (52)
X ij νn = d 3 r r n ν i (r)ν j (r)X(r), = 4π d r r n+2 ν i (r)ν j (r)X(r), (53)
where X stands for W , B, H, M, P , Q, R, or S. Again we see that whenever expansions of density matrices, Eqs. ( 40) and (44), are sufficiently accurate within the ranges of interactions, information about these interactions collapses to a few lowest moments. Short-range details of these interactions are, therefore, entirely irrelevant for low-energy characteristics of nuclear states. This is typical of all physical situations, where scales of interaction and observation are well separated, as specified in the effective field theories. The energy density characterizing the low-energy effects is local and depends on local densities and their derivatives up to second order, whereas the dynamic information is contained in a few coupling constants.

Moreover, the detailed large-r dependence of auxiliary functions ν i (r) on position r is also irrelevant, because all that matters are moments (53) which define the coupling constants (49)-(51) describing the exchange energy, and these are influenced only by the small-r properties of functions ν i (r). Finally, the most important feature is the k F or density dependence of ν i (r), which determines the density dependence of the coupling constants.

Local energy density corresponding to the Skyrme force

In general, the number of moments entering Eqs. ( 49)-( 51) is higher than the number of coupling constants, and all the coupling constants are independent. However, it is extremely instructive to check what happens in the vacuum limit of k F = 0. This situation is obtained by setting ν i (r) = 1, which gives the direct and exchange moments equal to one another, namely, X ij νn =X n , and the coupling constants of Eqs. ( 49)-(51) collapse to:

8        C ρ 0 C ρ 1 C s 0 C s 1        =        3 0 -1 -2 -1 2 -1 0          W 0 + M 0 B 0 + H 0   , (54) 96 
                              C ∆ρ 0 C ∆ρ 1 C τ 0 C τ 1 C ∆s 0 C ∆s 1 C T 0 C T 1 C F 0 C F 1 C ∇s 0 C ∇s 1 24   C ∇J 0 C ∇J 1   =   3 1   P 2 + Q 2 . ( 56
)
These coupling constants correspond exactly to those obtained for the Skyrme force (see Ref. [19] for the notations and conventions used), namely,

t 0 = W 0 + M 0 , t 0 x 0 = B 0 + H 0 , (57)
t 1 = -1 3 (W 2 + M 2 ) , t 1 x 1 = -1 3 (B 2 + H 2 ), (58)
t 2 = 1 3 (W 2 -M 2 ) , t 2 x 2 = 1 3 (B 2 -H 2 ), (59)
t e = 1 15 (S 2 -R 2 ) , t o = 1 15 (S 2 + R 2 ), (60)
W = -1 6 (P 2 + Q 2 ) . (61)
The same relations are also obtained by using in the exchange term the pure Taylor expansions (41); that is, by setting π(r) = 1, which gives X ij πn =X n , and by using the classification of terms as in Eqs. ( 21)-( 23). This second way of obtaining the approximate coupling constants leads to results independent of k F , which are, of course, identical to those obtained at k F = 0 above.

Relations (54)-( 56) imply that the coupling constants of the energy functional (48) are dependent of one another, and in fact, half of them determines the other half. This is exactly the situation encountered when the energy density is calculated for the Skyrme interaction. Then one obtains (cf. Ref. [23]):

3   C s 0 C s 1   =   -2 -3 -1 0     C ρ 0 C ρ 1   , (62) 24 
            C ∆s 0 C ∆s 1 C T 0 C T 1 C ∇s 0 C ∇s 1              =              -12 -12 3 9 0 -6 -4 -4 3 -3 -2 4 16 48 -4 12 -8 0 16 -16 4 -12 0 -8 0 0 0 0 0 -18 0 0 0 0 -6 12                           C ∆ρ 0 C ∆ρ 1 C τ 0 C τ 1 C F 0 C F 1              , (63)
C ∇J 0 = 3C ∇J 1 . (64)
It is obvious that the above relations among the coupling constants result from an oversimplified approximation to the exchange energy of the finite-range interaction.

We recall here [23,29] that without the tensor terms, relations (62) and (63) allow us to determine the time-odd coupling constants C s t , C ∆s t , and C T t as functions of the time-even coupling constants C ρ t , C ∆ρ t , and C τ t . Since the time-even coupling constants are usually adjusted solely to the time-even observables, the resulting values of the time-odd coupling constants are simply "fictitious" or "illusory", as noted already in Ref. [START_REF] Negele | Effective Interactions and Operators in Nuclei[END_REF]. In a more realistic case of relations (49) and (50), these constraints are no longer valid, and the time-odd properties of the functional are independent of the timeeven properties. This independence requires breaking the link between the Skyrme force Table 1. The NV coupling constants (49)-(50) calculated for the Gogny interaction D1S [17] for the Fermi momenta of k F = 0 and 1.35fm -1 . First-order coupling constants, C ρ t and C ∆s t , are in units of MeV fm 3 and the second-order coupling constants, C ∆ρ t , C τ t , C ∆s t , and C T t , are in units of MeV fm 5 . 59) calculated for the Gogny interaction D1S [17] for the Fermi momenta of k F = 0 and 1.35 fm -1 . These parameters correspond to the time-even sector of the Skyrme functional. Parameter t 0 is in units of MeV fm 3 ; parameters t 1 and t 2 are in units of MeV fm 5 , and parameters x 0 , x 1 , and x 2 are dimensionless.

t = 0 t = 1 k F = 0 k F = 1.35 k F = 0 k F = 1.
k F = 0 k F = 1.

Application to the Gogny interaction

In this section, we apply the results of Sec. 3 to the finite-range part of the Gogny interaction D1S [17]. This amounts to calculating moments (52) and (53) of the Gaussian functions with the two ranges of 0.7 and 1.2 fm, which constitute the central part of the Gogny interaction. Because this interaction does not contain any finite-range spin-orbit or tensor force in Eqs. ( 50) and (51), the moments P , Q, R, and S are set to zero. On the other hand, the zero-range spin-orbit and density-dependent terms of the Gogny interaction are left unchanged.

In Table 1 we show values of coupling constants (49)-( 51) calculated in the vacuum (k F = 0) and at the saturation density (k F = 1.35 fm -1 ). Similarly, Table 2 shows show values of the isoscalar and isovector coupling constants, respectively. Units are specified in the caption to Table 1.

values of the Skyrme-force parameters (57)-( 59) corresponding to the time-even coupling constants. In Figs. 2 and3 we plot the coupling constants and Skyrme-force parameters, respectively, as functions of the Fermi momentum.

The most important observation resulting from values shown in Tables 1 and2 and Figs. 2 and3 pertains to a significant density (or k F ) dependence of the coupling constants and Skyrme-force parameters. The strongest dependence is obtained for the isoscalar tensor coupling constant C T 0 . (Note that the central finite-range Gogny interaction induces significant values of the tensor coupling constants C T t , even if this force does not contain any explicit tensor term.) Also the kinetic coupling constants C τ t exhibit a strong density dependence, going almost to zero at k F ∼ 2 fm -1 . We note that the obtained density dependencies do not, in general, follow any power laws. Significantly stronger density dependencies are obtained for the Skyrme-force parameters (Fig. 3). The pole appearing in the parameter x 2 is a consequence of the fact that parameters t 2 and t 2 x 2 , derived in Eq. ( 59), change signs at slightly different values of the Fermi momentum.

Let us now consider the question of whether parameters calculated at any fixed value of the Fermi momentum can provide a reasonable alternative. To analyze this point, we first note that the zero-order coupling constants given in Eq. (49) depend on 59), calculated for the Gogny interaction D1S [17], on the Fermi momentum k F . These parameters correspond to the time-even sector of the Skyrme functional. Units are specified in the caption of Table 2.

k F (i) implicitly through the moments of Eq. ( 53) and (ii) explicitly through the k 2 F term, that is,

C ρ t (k F ) = C ρ t,0 (k F ) + C ρ t,2/3 (k F )ρ 2/3 , (65)
C s t (k F ) = C s t,0 (k F ) + C s t,2/3 (k F )ρ 2/3 , (66)
where we employed the standard association of the Fermi energy with density, namely, ρ = 2k 3 F /3π 2 . Therefore, we can consider fixed values of the coupling constants calculated at a given Fermi momentum k 0 F as

C ρ t (k 0 F ) and C s t (k 0 F ), or (67)
C ρ t,0 (k 0 F ), C ρ t,2/3 (k 0 F ), C s t,0 (k 0 F ), and C s t,2/3 (k 0 F ). (68)
The second option gives the coupling constants that still depend on the density as ρ 2/3 , in analogy with the standard density-dependent term, which for the Gogny force depends on the density as ρ 1/3 , and which in the present study is alway kept untouched.

In the upper and lower panels of Fig. 4, thin lines show the nuclear matter equations of state (energy per particle in function of density or Fermi momentum) obtained for the coupling constants fixed according to prescriptions (67) and (68), respectively. For comparison, thick lines show the NV results obtained for density-dependent coupling constants. Since in the nuclear matter, the factor multiplying the term ν 0 (r)ν 2 (r) in Eq. ( 18) vanishes exactly (by construction), the thick lines correspond to the exact Gogny force results.

Equations of state calculated with prescription (68) for fixed values of k 0 F from 1.35 to 1.98 fm -1 with the step of 0.07 fm -1 completely miss the saturation point. This means F . Upper and lower panels show results for fixed coupling constants (67) and for coupling constants (68) that depend on density as ρ 2/3 . The full square and circles mark the minima of curves.

that Skyrme forces with constant parameters (68) derived by the NV expansion cannot be equivalent to the finite-range Gogny interaction. On the other hand, prescription (67), for fixed value of k 0 F = 1.35 fm -1 , reproduces the equation of state fairly well, with some deviations seen only at densities beyond the saturation point. For comparison, we also show results obtained with k 0 F = 1 and 1.7 fm -1 , which fit the equation of state at low and high densities, respectively. These results show that the explicit dependence on the Fermi momentum, which appears in the NV expansion, cannot be used to define the density dependence of the coupling constants.

The full density dependence of all Skyrme-force parameters was recently implemented in a spherical self-consistent code [31]. Here we use this implementation to test the results of the NV expansion against the full-fledged solutions known for the Gogny D1S force [32,33]. In Table 4, we show results obtained for three sets of Skyrme-force parameters:

• S1Sa: fixed parameters given in Table 2 for k 0 F = 1.35 fm -1 , used together with the standard D1S parameters [17] of the zero-range spin-orbit (W = 130 MeV fm 5 ) and density-dependent (t 3 = 1390.60 MeV fm 4 ) terms.

• S1Sb: density-dependent parameters shown in Fig. 3, implemented for ρ = 2k 3 F /3π 2 and used with the standard values of W and t 3 .

• S1Sc: fixed parameters identical to those of S1Sa, but with the value of t 3 = 1385. 35 MeV fm 4 slightly changed, so as to bring the ground-state energy of 208 Pb exactly to the true D1S value [33].

Table 4 shows the ground-state energies E of seven doubly magic spherical nuclei, calculated by using the Skyrme-force parameters S1Sa, S1Sb, and S1Sc, and compared with the Gogny-force energies E G . To facilitate the comparison, we also show relative differences ∆E = (E -E G )/|E G | in percent along with the corresponding RMS deviations (the last row in Table 4). On can see that for the Skyrme-force parameters S1Sa and S1Sb, which are directly derived from the Gogny force by using the NV expansion, one obtains the nuclear binding energies smaller by 1-2% as compared to those given by the original Gogny force. This should be considered a very good result, although it cannot compete in precision of describing experimental data with the original Gogny or Skyrme forces, which have parameters directly fitted to experimental binding energies. A simple rescaling of the parameter t 3 brings the RMS deviation to below 0.4%, and makes the Skyrme force S1Sc competitive with most other standard Skyrme parameterizations. At the same time, the analogous RMS deviations obtained for the neutron and proton radii are 0.20 and 0.30% (S1Sa), 1.01 and 0.94% (S1Sb), and 0.26 and 0.44% (S1Sc), respectively.

We note here in passing that in this study we defined the S1Sb parameter set by considering the density-dependent coupling constants C ∆ρ t (Fig. 3) that multiply densities ρ t ∆ρ t . An attempt of using the same coupling constants along with densities -(∇ρ t ) 2 gives, in fact, the RMS deviations of binding energies (5.34%), neutron radii (2.44%), and proton radii (2.27%), which are significantly worse than those of S1Sb (Table 4). This shows that the prescription to replace in EDFs integrated by parts (see Refs. [29,31]) the Fermi momentum by ρ = 2k 3 F /3π 2 may lead to significantly different results. Another recipy is to associate k F with density before taking products of the two density matrices which means both of the above terms will be active. However in order to have a better correspondence with the traditional Skyrme functionals we have made this association after taking the product.

It is not our purpose here to propose that any of the Skyrme-force parameterizations introduced in the present work are better solutions to the problem of finding the best agreement with data. It is already known that within the standard second-order Skyrme-force parameterizations, a spectroscopic-quality [2] force cannot be found [34]. Nevertheless, it is gratifying to see that the NV expansion allows us to bridge the gap between the non-local and quasi-local EDFs, or between the finite-range and zero-range effective forces. A more quantitative discussion of the accuracy of the NV expansion will be possible by considering higher-order NV expansions [35].

Finally, in Fig. 5 we compare the time-odd coupling constants calculated by using Eqs. ( 49) and (50) with those corresponding to the Skyrme-force parameters; that is, calculated by using Eqs. ( 62) and (63). Similarly, Table 3 lists the numerical values of the Skyrme-force time-odd coupling constants. As one can see, differences between both sets of the time-odd coupling constants, shown in Fig. 5 with open and full symbols, are quite substantial. These results illustrate the fact that the NV expansion of the Gogny force leads to the Skyrme functional and not to the Skyrme force.

We conclude this section by noting that functions π i (r) approximated by Gaussians (see Eqs. (36)-(38) and Fig. 1) lead to the coupling constants and Skyrme-force parameters, which, when plotted in the scales of Figs. 2,3, and 5, are indistinguishable from those presented in these figures.

Conclusions

In the present study, we derived a set of compact expressions giving the local energy density corresponding to the Kohn-Sham potential energy for an arbitrary local finiterange central, spin-orbit, and tensor interactions. The method is based on the Negele-Vautherin density matrix expansion augmented by the odd-power gradient terms fulfilling the gauge-invariance condition. The coupling constants of the local energy density depend on a set of moments of the interaction conforming to the ideas of the effective theory. The expansion is based on the separation of scales between the range of the force and space characteristics of the one-body density matrix. It leads to a representation of dynamical properties of the system in terms of a set of numbers, whereby complicated short-range characteristics of effective interactions remain unresolved. We pointed out the fact that to correctly describe the exchange properties of the functional, proper treatment of the density matrix in the nonlocal direction is essential. This immediately leads to the local energy density that does not correspond to an averaged zero-range pseudopotential. Therefore, the Negele-Vautherin expansion performed up to NLO leads to the Skyrme functional and not to the Skyrme force. Within this formalism, the only way to define the Skyrme force is to match it to the time-even properties of the non-local functional, disregarding those pertaining to the time-odd channel.

We applied the general NLO expressions to the case of the central finite-range part of the Gogny interaction. It turns out that the obtained coupling constants of the local Skyrme functional quite strongly depend on the Fermi momentum or on the density. Nevertheless, the equation of state obtained for fixed coupling constants, calculated at the saturation point of k F = 1.35 fm -1 , fairly well reproduces the exact Gogny-force result. On the other hand, partial density dependence, inferred from the explicit dependence of the coupling constants on k F , gives very unsatisfactory results.

By solving the self-consistent equations with the Skyrme-force parameters derived from the Gogny force, one obtains an excellent agreement (up to 1-2%) of binding energies and radii with those corresponding to the true Gogny force. This shows that the ideas of the effective theory, whereby the finite-range nuclear forces are sufficiently short-range to be replaced by contact quasi-potentials, are applicable to low-energy nuclear observables.

We also discussed properties of the one-body density matrix for the spin and isospin polarized infinite nuclear matter. In this case, one obtains the nonlocal densities with mixed spin-isospin channels. As a result, the local energy density is not invariant but covariant with respect to rotational and isospin symmetries; that is, it does not have the form of the standard Skyrme functional. It means that the standard Negele-Vautherin expansion can only be performed for unpolarized densities. the relative position vector r is given by four scalar functions πστ (r) [compare Eq. ( 30)], πστ (r) = 3j 1 (k F,στ r) k F,στ r . (A.3)

We note that the spin-isospin indices στ pertain here to the preselected quantization axis defined by the chosen directions of H and λ, which define an "intrinsic" reference frame. In this reference frame, densities are marked with a bar symbol. The ground state of the system is obtained by filling the four Fermi spheres up to the common Fermi energy F ,

F = F,στ = h2 k 2 F,στ 2m
-H z σ -λ 3 τ, (A.4) which defines the four Fermi momenta k F,στ . Finally, by varying F , one obtains systems with different total densities ρ = στ ρστ .

It is, of course, clear that for the asymmetric and polarized infinite nuclear matter, the density matrix of Eq. ( 39) is diagonal in spin and isospin, ρ(r 1 σ 1 τ 1 , r 2 σ 2 τ 2 ) = ρσ 1 τ 1 (r 1 , r 2 )δ σ 1 σ 2 δ τ 1 τ 2 (A.5)

and thus the nonlocal densities ρµk (r 1 , r 2 ) have non-zero components only for µ = 0 or z and k = 0 or 3, that is, where functions πµk (r) are defined similarly as in Eq. (A.6), namely,

       ρ00 (r 1 , r 2 ) ρ03 (r 1 , r 2 ) ρz0 (r 1 , r 2 ) ρz3 (r 1 , r 2 )        =        1 1 1 1 1 -1 1 -1 1 1 -1 -1 1 -1 -1 1               ρ++ (
       π00 (r) π03 (r) πz0 (r) πz3 (r)        = 1 4        1 1 1 1 1 -1 1 -1 1 1 -1 -1 1 -1 -1 1               π++ (r) π+-(r) π-+ (r) π--(r)        . (A.8)
Already here we see the main problem: for the spin and isospin polarized systems, the spin-isospin channels of nonlocal densities ρµk (r 1 , r 2 ) in Eq. (A.7) are linear combinations of the spin-isospin channels of local densities ρµk ; that is, the spin-isospin channels become mixed.

To make the preceding result even more clear, we note that the spin-isospin directions of the Lagrange multipliers H and λ can be arbitrarily varied and the spinisospin directions of the nonlocal densities ρ µk (r 1 , r 2 ), local densities ρ µk , and functions

Figure 1 .

 1 Figure 1. Dependence of the functions π 0 (r) (solid lines), π 2 (r) (dashed lines), and ν 0 (r) (dotted lines) on k F r. The top and middle panels show the NV functions of Eqs. (31)-(33), with the top panel plotted in expanded scale to better show the details at large k F r. The bottom panel shows the functions π 0 (r), π 2 (r), and ν 0 (r) approximated by Gaussians as in Eqs. (36)-(38).

Figure 2 .

 2 Figure 2. Dependence of the NV coupling constants (49)-(50), calculated for the Gogny interaction D1S[17], on the Fermi momentum k F . Full and open symbols (lower and upper panels) show values of the time-even and time-odd coupling constants, respectively. Solid and dashed lines (left and right panels) show values of the isoscalar and isovector coupling constants, respectively. Units are specified in the caption to Table1.

2 Figure 3 .

 23 Figure 3. Dependence of the standard Skyrme-force parameters (57)-(59), calculated for the Gogny interaction D1S[17], on the Fermi momentum k F . These parameters correspond to the time-even sector of the Skyrme functional. Units are specified in the caption of Table2.

Figure 4 .

 4 Figure 4. Dependence of the infinite-matter energy per particle on the density (left panels) or Fermi momentum (right panels). Thick lines show the NV results for density-dependent coupling constants. Thin lines show the NV results for fixed coupling constants calculated at the given values of k 0 F . Upper and lower panels show results for fixed coupling constants (67) and for coupling constants (68) that depend on density as ρ 2/3 . The full square and circles mark the minima of curves.

Figure 5 .

 5 Figure 5. Similar to Fig. 2 except for the time-odd NV coupling constants only. The open symbols show exact results (49)-(50) as plotted in Fig.2, whereas the full symbols show results inferred from the time-even sector by using Eqs. (62) and (63). Solid and dashed lines (left and right panels) show values of the isoscalar and isovector coupling constants, respectively.

r 1 , 1 ,

 11 r 2 ) ρ+-(r 1 , r 2 ) ρ-+ (r 1 , r 2 ) ρ--(r 1 , r 2 ) abbreviated the indices of στ just to their signs. After expressing the right-hand side of this equation in terms nonlocal densities (A.2), one obtains: r 2 ) ρ03 (r 1 , r 2 ) ρz0 (r 1 , r 2 ) ρz3 (r 1 , r 2 ) π03 (r) πz0 (r) πz3 (r) π03 (r) π00 (r) πz3 (r) πz0 (r) πz0 (r) πz3 (r) π00 (r) π03 (r) πz3 (r) πz0 (r) π03 (r) π00 (r)

Table 2 .

 2 The standard Skyrme-force parameters (57)-(

	35

Table 3 .

 3 Similar to Table1except for the time-odd NV coupling constants inferred from the time-even sector by using Eqs. (62) and (63) for a Fermi momentum of 1.35 fm -1 .

		t = 0	t = 1
	C s t	-27.94756	200.2052
	C ∆s t	20.92673	23.21032
	C T t	-26.47083 -44.78631

Table 4 .

 4 Binding energies E of seven doubly magic nuclei calculated by using the Skyrme-force parameters S1Sa, S1Sb, and S1Sc (see text) compared with the Gognyforce energies E G . All energies are in MeV.

		D1S [33]	S1Sa		S1Sb		S1Sc
		E G	E	∆E	E	∆E	E	∆E
	40 Ca	-342.689	-335.312 2.15%	-340.642 0.60%	-339.369	0.97%
	48 Ca	-414.330	-409.118 1.26%	-410.698 0.88%	-414.213	0.03%
	56 Ni	-481.111	-473.497 1.58%	-471.970 1.90%	-479.843	0.26%
	78 Ni	-637.845	-630.447 1.16%	-629.066 1.38%	-638.837 -0.16%
	100 Sn	-828.024	-814.568 1.63%	-814.896 1.59%	-826.453	0.19%
	132 Sn -1101.670 -1086.272 1.40% -1086.867 1.34% -1101.445	0.02%
	208 Pb -1637.291 -1612.634 1.51% -1617.419 1.21% -1637.291	0.00%
	RMS	n.a.	n.a. 1.56%	n.a. 1.33%	n.a.	0.39%

Interesting comments by Thomas Duguet are gratefully acknowledged. This work was supported by the Polish Ministry of Science and Higher Education under Contract No. N N 202 328234, by the Academy of Finland and University of Jyväskylä within the FIDIPRO program, and by the U.S. Department of Energy under Contract Nos. DE-FC02-09ER41583 (UNEDF SciDAC Collaboration) and DE-FG02-96ER40963 (University of Tennessee).

Appendix A. Spin and isospin polarized infinite nuclear matter

Let us consider the infinite nuclear matter with spin and isospin polarizations, which is described by the one-body Hamiltonian,

where σ and H are the space vectors of spin Pauli matrices and spin-polarization Lagrange multipliers, respectively, and τ and λ are the isovectors of the analogous isospin Pauli matrices and isospin-polarization Lagrange multipliers, whereas the dot "•" (circle "•") denotes the scalar (isoscalar) product. Each eigenstate of Hamiltonian (A.1) is a Slater determinant that depends on the orientations of the Lagrange multipliers H and λ in the space and isospace, respectively. However, since the kinetic energy is scalar and isoscalar, we can arbitrarily fix these orientations to H z and λ 3 , which gives the nonlocal densities for spin-up and spin-down (σ = ±1) neutrons and protons (τ = ±1) in the form,

The system simply separates into four independent Fermi spheres for spin-up and spindown neutrons and protons, with four constant densities ρστ , whereas the dependence on π µk (r) are always aligned with those of the Lagrange multipliers. Therefore, we can use the directions of the local densities instead of those pertaining to the Lagrange multipliers. By using the standard densities [19] in the (i) scalar-isoscalar channel (ρ = ρ 00 ), (ii) vector-isoscalar channel (s µ = ρ µ0 , for µ = x, y, z), (iii) scalar-isovector channel ( ρ k = ρ 0k , for k = 1, 2, 3), and (iv) vector-isovector channel ( s µk = ρ µk , for µ = x, y, z and k = 1, 2, 3), we then define functions π(r) in the four channels as,

| s| πz3 (r).

(A.9)

Here, the "intrinsic" functions πµk (r) do not depend on the spin-isospin directions; that is, they are defined by the following Fermi energies,

Finally, definitions (A.9) allow us to present densities in the "laboratory" reference frame as [compare Eq. (A.7)], .14). Therefore, the postulate of using different functions in different channels [11] is not compatible with the results obtained for the polarized nuclear matter. Again we see that the spin-isospin channels of nonlocal densities are mixed, namely, local densities in all channels contribute to every channel in the nonlocal density. As a consequence, the energy density is not invariant but only covariant with respect to the spin-isospin rotations (see the discussion in the Appendix A of Ref. [5]). Therefore, the NV expansion performed in the polarized nuclear matter does not lead to the standard local functional of Eq. (48). On the other hand, derivation in the unpolarized nuclear matter corresponds to all functions ρστ (r) equal to one another, which leads to vanishing functions π(r), π(r), and π(r). Then, in Eqs. (A.11)-(A.14), only the first terms survive and the spin-isospin channels are not mixed. Such a situation corresponds to postulating a channel-independent function π(r), which we employed in Sec. 3.


Alexis Diaz-Torres
email: a.diaztorres@surrey.ac.uk

Modelling incomplete fusion dynamics of weakly-bound nuclei at near-barrier energies

Keywords: numbers: 25.60.Pj, 25.60.Gc, 25.60.-t, 24.10.-i

The classical dynamical model for reactions induced by weakly-bound nuclei at near-barrier energies is developed further. It allows a quantitative study of the role and importance of incomplete fusion dynamics in asymptotic observables, such as the population of high-spin states in reaction products as well as the angular distribution of direct alpha-production. Model calculations indicate that incomplete fusion is an effective mechanism for populating high-spin states, and its contribution to the direct alpha production yield diminishes with decreasing energy towards the Coulomb barrier. It also becomes notably separated in angles from the contribution of no-capture breakup events. This should facilitate the experimental disentanglement of these competing reaction processes.

Introduction

Nuclear physics research has entered a new era with developments of radioactive nuclear beam facilities, where nuclear reactions are the primary probe of the new physics, such as novel structural changes. In those facilities, the low-energy nuclear reactions research is highly focused on understanding astrophysically important reaction rates involving exotic nuclei. These are often weakly-bound with a few-body, cluster structure that can easily be dissociated in their interaction with other nuclei. Understanding the breakup mechanism and its impact on nuclear reaction dynamics is essential. A major consequence of breakup is that a rich scenario of reaction pathways arises, such as events where (i) not all the resulting breakup fragments might be captured by the target, termed incomplete fusion (icf), (ii) the entire projectile is captured by the target, called complete fusion (cf), and (iii) none of the breakup fragments are captured, termed nocapture breakup (ncbu).

Since the availability of intense exotic beams is still limited, extensive experimental research has recently been carried out exploiting intense beams of stable weakly-bound nuclei, such as 6,7 Li and 9 Be [1,2]. Understanding the effect of their breakup on nearbarrier fusion has been a key aspect of these investigations [3]. These have definitively demonstrated that breakup suppresses the above-barrier fusion cross sections. Most recently, experimental activities have been focused on disentangling breakup and competing reaction mechanisms from inclusive and exclusive coincidence measurements [4,5,6,7,8]. One of the challenges is to obtain a complete quantitative understanding of the breakup mechanism and its relationship with near-barrier fusion. This research is guided by complete sub-barrier breakup measurements [9].

Theoretical works have addressed the low-energy reaction dynamics of weaklybound nuclei using quantum mechanical, classical and mixed quantum-classical approaches [10,11,12,13,14,15,16,17]. Among these, the continuum-discretised coupled channels (cdcc) framework has been very successful [10,13,16,17]. However, existing quantum models have limitations [18], as they cannot calculate integrated icf and cf cross sections unambiguously. Neither, after the formation of icf products, can they follow the evolution of the surviving breakup fragment(s) since icf results in depletion of the total few-body wave function.

These difficulties are overcome by the three-dimensional classical dynamical reaction model suggested in Ref. [19]. A crucial input of this model is a stochastically sampled breakup function proposed in Ref. [20], which can be determined from subbarrier breakup measurements [9,20]. This function encodes the effects of the Coulomb and nuclear interactions that cause the projectile breakup. Hence, this approach is not a breakup model, rather it is a quantitative dynamical model for relating the sub-barrier ncbu to the above-barrier icf and cf of weakly-bound nuclei [20]. The model has been successfully applied to interpreting fusion measurements of weakly-bound nuclei [7,9], and isomer ratio measurements [21].

In this paper, I report on further developments of this classical dynamical model. The key new aspect is the time propagation of the surviving breakup fragment and the icf product, allowing the description of their asymptotic angular distribution and dynamical variables. These should be very useful in (i) current experimental activities aimed at disentangling competing reaction mechanisms from asymptotic observables such as alpha-production yields [8], (ii) interpreting particle-γ-coincidence experiments [22], and (iii) applications to γ ray spectroscopy [23,24,25]. The new developments are illustrated with a simplified test case that does not aim at adjusting any measurements. It is worth mentioning that various dynamical models have been proposed for multi-fragmentation and icf in heavy-ion induced reactions at energies well-above the Coulomb barrier ( 10 MeV/nucleon), as recently reviewed in Ref. [22]. In this different context, some of those models (e.g., Refs. [26,27]) use concepts and techniques (e.g., classical trajectory, Monte Carlo sampling) that overlap with some involved in the present model. The model is explained in Section 2. In Section 3, numerical results are discussed, and a summary is given in Section 4.

Model

The main features of the model are as follows:

(i) The target T is initially at rest in the origin of the laboratory frame, whilst the weakly-bound (two-body) projectile P approaches the target (along the z-axis) with incident energy E 0 and orbital angular momentum L 0 . For each L 0 (chosen to be an integer number of ) an ensemble of N incident projectiles is considered.

Including the P -T mutual Coulomb and nuclear forces, classical equations of motion determine an orbit with a definite distance of closest approach R min (E 0 , L 0 ).

(ii) The complexity of the projectile dissociation is empirically encoded in a density of (local breakup) probability P L BU (R), a function of the projectile-target separation R, such that P L BU (R)dR is the probability of breakup in the interval R to R + dR (see Appendix A). A key feature is that for a given projectile-target combination, both measurements [9,20] and cdcc calculations [19] indicate that the integral of this breakup probability density along a given classical orbit is an exponential function of its distance of closest approach, R min (E 0 , L 0 ):

P BU (R min ) = 2 ∞ R min P L BU (R)dR = A exp(-αR min ). (1)
Consequently, P L BU (R) has the same exponential form, P L BU (R) ∝ exp(-αR). (The factor of 2 highlights that breakup may occur along the entrance or exit branch of the trajectory, although the exponential form will clearly place maximum probability of breakup at R min .) This function is sampled to determine the position of breakup in the orbit discussed in (i). In this position, the projectile is instantaneously broken up into fragments F1 and F2. These interact with T , and with each other, through real central two-body potentials having Coulomb barriers

V ij B at separations R ij B , i, j = 1, 2, T, i = j. (iii)
The instantaneous dynamical variables of the excited projectile at breakup, namely its total internal energy ε 12 , its angular momentum 12 and the separation of the fragments d 12 are all Monte Carlo sampled. The initial separation d 12 between the fragments and its orientation d 12 are determined by sampling the radial and angular probability distributions of the projectile ground-state (g.s.) wave function, respectively. For a two-body projectile with 0 + g.s., a very good approximation for calculating d 12 is through a Gaussian sampling function in the classically allowed region of the fragments, whilst the orientation of d 12 is isotropic. (This will be used in the test case below.) For high 12 excitations, when there is no barrier between F1 and F2, d 12 is equated with their external turning point. The orientation of 12 is chosen randomly from all directions orthogonal to d 12 . 12 is sampled uniformly in the interval [0, max ], whilst for ε 12 an exponentially decreasing function for energies between the top of the barrier (V 12 B ) and a chosen maximum ε max is sampled. Both max and ε max are increased until convergence of the observables occur. (iv) Having fixed the position and dynamical variables of the excited projectile fragments at the moment of breakup, the instantaneous velocity of the particles F1, F2 and T is determined by conservation of energy, linear momentum and angular momentum in the overall center-of-mass frame (see Appendix B). These breakup initial conditions are transformed to the laboratory frame where the three bodies are propagated in time. The calculated trajectories of F1, F2 and T determine the number of icf, cf and ncbu events, fragment Fj being assumed to be captured if the classical trajectories take it within the fragment-target barrier radius R jT B . (v) From the N breakup events sampled for each projectile angular momentum L 0 , the numbers of events N i in which i = 0 (ncbu), 1 (icf), or 2 (cf) fragments are captured determine the relative yields P i = N i /N of these three reaction processes after breakup, with P 0 + P 1 + P 2 = 1. The absolute probabilities P i (E 0 , L 0 ) of these processes are expressed in terms of the relative yields and the integrated breakup probability over the whole trajectory P BU (R min ):

P 0 (E 0 , L 0 ) = P BU (R min ) P 0 , (2)
P 1 (E 0 , L 0 ) = P BU (R min ) P 1 , (3)
P 2 (E 0 , L 0 ) = [1 -P BU (R min )] H(L cr -L 0 ) + P BU (R min ) P 2 , (4)
where H(x) is the Heaviside step function and L cr is the critical partial wave for projectile fusion. The cross sections are calculated using

σ i (E 0 ) = πλ 2 L 0 (2L 0 + 1)P i (E 0 , L 0 ), (5)
where λ 2 = 2 /[2m P E 0 ] and m P is the projectile mass. Beside the absolute cross sections (5), asymptotic observables, such as the angle, kinetic energy and relative energy distributions of the fragments from ncbu events, are calculated by tracking their trajectories to a large distance from the target.

(vi) For the icf events, the time propagation of the icf product and the surviving breakup fragment is now incorporated into this model as follows. The fragment Fj (after overcoming the Coulomb barrier V jT B ) reaches the target radius, forming the icf product, while the other fragment flies away. At this moment, the three-body propagation turns into a two-body propagation, with definite interaction potentials and initial conditions. These are determined by the position and velocity of the three particles, at the moment when the icf product is formed. These also yield the spin and excitation energy distributions of the primary icf product. The asymptotic angular distribution of the icf product and the surviving breakup fragment is calculated in terms of their trajectories.

Numerical results

This model is implemented in the platypus code, as described in Ref. [START_REF] Diaz-Torres | PLATYPUS: a code for fusion and breakup in reactions induced by weaklybound nuclei within a classical trajectory model with stochastic breakup[END_REF]. In order to illustrate the new developments of (vi), calculations are carried out at a laboratory energy of E 0 = 45 and 65 MeV, for the reaction of a pseudo-8 Be projectile P (modeled as a weakly-bound s-state of two α-particles [19]) with the 208 Pb target T .

The breakup function P BU (R min ) in eq. ( 1) has parameters A = 5.98 × 10 3 and α = 0.85 fm -1 , which was deduced in Ref. [19] from mapping fusion measurements for With decreasing energy towards the Coulomb barrier, the ncbu events dominate, separating its centroid substantially from that of the icf events. The total alpha-production distribution changes its shape notably.

the 9 Be + 208 Pb system with this test reaction. This is because the model is limited to solving a three-body problem. However, very recent measurements [9] have shown that prompt 9 Be breakup occurs dominantly through an excited 8 Be nucleus, validating the approximation of a 9 Be projectile by 8 Be.

The nuclear interaction between the alpha particle and the icf product 212 Po is considered to be the Woods-Saxon potential (V , r, a) ≡ (33.98 MeV, 1.48 fm, 0.63 fm) resulted from the global Broglia-Winther parametrization [START_REF] Reisdorf | [END_REF]. (Please note that in the potential the radius parameter is multiplied by A 1/3 T .) The rest of the model parameters are the same as in Ref. [19].

Figure 1 shows the angular distribution of direct alpha-production for two laboratory energies near the P -T s-wave Coulomb barrier (39.9 MeV), namely (a) E 0 = 65, and (b) 45 MeV. The contribution of the icf and ncbu events is represented by thick solid and thick dashed lines, respectively. Their incoherent sum is represented by the thin solid line. Its shape significantly changes as the incident energy decreases. While the contribution of the icf and ncbu events appears to be similar at well-above Here, its centroid significantly separates from the centroid of the icf contribution. Both centroids shift to higher angles as the incident energy decreases, due to the reduction of relative partial waves affecting these reaction processes. However, the ncbu centroid always remains lower than the icf centroid, as expected. This is because higher partial waves contribute to the ncbu process (see Fig. 2). Fig. 2 presents the incident angular momentum distribution of the icf (solid line) and ncbu (dashed line) processes at a laboratory energy of E 0 = 45 MeV. Figure 3 shows the angular momentum distribution of primary icf ( 212 Po) and cf ( 216 Rn) products at (a) E 0 = 65, and (b) 45 MeV. The icf spin distribution is represented in terms of both the 8 Be incident angular momentum L 0 (thick solid line) and the angular momentum brought in by the α-particle into the 208 Pb target (thin solid line). These representations result in two very different icf angular momentum distributions, both providing relevant features of the icf dynamics as explained below.

Comparing the thick solid line with the cf spin distribution (dashed line) it is seen that the localization of the L 0 -window for icf significantly depends on the incident energy E 0 . Their overlap strongly increases as E 0 decreases towards the barrier [panel (b)], indicating that icf and cf are two competing reaction processes at near-barrier energies. However, cf following 8 Be breakup is here a very small component (1.5%) of the cf cross section (149.2 mb). Instead, either one of the α-particles is captured by the target, contributing to the icf cross section (84.7 mb), or the two α-particles survive and contribute to the ncbu cross section (238.8 mb). At E 0 = 65 [panel (a)], cf following 8 Be breakup is very substantial, representing 41% of the cf cross section (981.5 mb). This significantly shifts down the maximum of the cf spin distribution (dashed line) with respect to the critical angular momentum L cr for 8 Be fusion as an inert projectile (thick arrow). Here, L 0 values around L cr determine the icf cross section (273.6 mb) which is similar to the ncbu cross section (259.3 mb). spin distribution is shown in two representations, i.e., in terms of (i) the 8 Be incident angular momentum L 0 (thick solid line), and (ii) the angular momentum brought in by the α-particle into the 208 Pb target (thin solid line). The thick arrow denotes L cr for 8 Be fusion as an inert projectile, whereas the thin arrow is l cr for fusion of a direct beam of α-particles (with half of E 0 ) on 208 Pb. The relevant features are:

(1) the icf mechanism increases the angular momentum brought in by the α-particle into 208 Pb (thin solid line goes beyond the thin arrow), (2) the relative localization of the L 0 -window for cf and icf significantly depends on E 0 (dashed and thick solid lines), strongly overlapping as E 0 decreases, and (3) cf following 8 Be breakup becomes substantial with increasing E 0 , shifting down the maximum of the cf spin distribution (dashed line) with respect to the thick arrow.

The icf spin distributions (thin solid line) in Fig. 3 are the crucial ones in order to access the effectiveness of the icf mechanism for populating high-spin states in 212 Po. It is observed that the tail of these distributions goes well beyond the thin arrow which denotes l cr for fusion, on the 208 Pb target, of a direct beam of α-particles with half of E 0 . Clearly, the extra torque caused by the interaction between the fusing and surviving α-particles enhances the angular momentum of the icf product. This interaction also affects the excitation energy of 212 Po substantially, as presented in Fig. 4. MeV. The arrow denotes the excitation energy of 212 Po, when this compound nucleus is formed by fusion of a direct beam of α-particles, with half of E 0 , on 208 Pb. Following breakup of 8 Be, the fusing α-particle is significantly affected by the interaction with the surviving α-particle. This spreads the 212 Po excitation energy over a range of values. The distribution becomes narrower and its maximum approaches the arrow, as E 0 decreases towards the barrier.

(thin solid line) at (a) E 0 = 65, and (b) 45 MeV. It is caused by the interaction between the fusing and surviving α-particles during the icf process. The arrow corresponds to the pre-determined excitation energy of 212 Po, when this nucleus is formed through cf of a direct beam of α-particles (with half of E 0 ) on 208 Pb. This value remains smaller (by a few MeV) than the excitation energy associated with the maximum of the distribution. The shape of the distribution also changes with E 0 , becoming narrower as E 0 decreases towards the barrier. These distributions along with the icf spin distributions in Fig. 3 (thin solid line) are vital for a reliable prediction of the yield of final icf products. To my knowledge, these aspects have not been yet included in broadly used evaporation codes which mostly describe cf evaporation residues. The present model can be applied to more realistic cases (e.g., 11 Be induced reactions), provided all the necessary information for Monte Carlo sampling is known (e.g., the breakup function and the projectile g.s. wave-function). When one of the breakup fragments is neutral, its capture by the target nucleus can be assumed to occur when its trajectory takes it within the target radius. Although not implemented yet, the model can also predict the elastic angular distribution of the projectile. This distribution can easily be calculated in terms of the projectile-target orbits and the probability of the weakly-bound projectile's survival (see Appendix A), i.e., this probability weighting the contribution of the different projectile-target orbits. Direct reaction processes other than the elastic breakup may significantly contribute to observables associated with ncbu, cf and icf. These are not included in the model yet. However, their contributions could be included using transfer and inelastic-breakup functions. Of course, these would make the calculations more complex, as these direct reaction processes result in additional bifurcation points along the projectile-target orbits. Nevertheless, their inclusion as well as the treatment of more complex projectiles (e.g., 6 He and 11 Li) will be interesting, necessary developments of the present approach.

Summary

The classical dynamical model for reactions of weakly-bound nuclei at near-barrier energies has been developed further. It allows us to quantify the role and importance of icf dynamics in asymptotic reaction observables, such as the angular distribution of the direct alpha-production. The icf contribution to this yield diminishes with decreasing energy towards the barrier. However, the icf and ncbu contributions are clearly separated in angles at near-barrier energies, facilitating the experimental disentanglement of these competing reaction processes. The present developments also provide the spin and excitation energy distributions of primary icf products, which are essential ingredients for calculating the yield of final icf products with present evaporation codes. The icf mechanism appears to be an effective route for producing high-spin states. All these observables may also be affected by other direct processes, such as transfer [30], which are not included in the model yet. Nevertheless, the present classical dynamical model is a powerful tool for interpreting fusion measurements involving radioactive nuclei and in applications to γ ray spectroscopy. The development of a unified quantum dynamical description of relevant reaction processes (icf, cf, ncbu and transfer) remains a great theoretical challenge. One possibility of tackling this issue could be through a time-dependent density-matrix approach incorporating the concept of quantum decoherence [18,31].

Figure 1 .

 1 Figure 1. Angular distribution of direct alpha-production for 8 Be + 208 Pb at two laboratory energies E 0 : (a) 65 MeV, and (b) 45 MeV.With decreasing energy towards the Coulomb barrier, the ncbu events dominate, separating its centroid substantially from that of the icf events. The total alpha-production distribution changes its shape notably.

Figure 2 .

 2 Figure 2. icf and ncbu cross sections as a function of the relative angular momenta L 0 for 8 Be + 208 Pb at E 0 = 45 MeV. The contribution of high-partial waves shifts the ncbu distributions in Fig. 1 to smaller angles, with respect to the icf distributions.

Figure 3 .

 3 Figure 3. Spin distribution of primary fusion products resulting from the 8 Be + 208 Pb reaction at two laboratory energies E 0 : (a) 65 MeV, and (b) 45 MeV. The cf ( 216 Rn) spin distribution is presented by the dashed line, whilst the icf ( 212 Po)spin distribution is shown in two representations, i.e., in terms of (i) the8 Be incident angular momentum L 0 (thick solid line), and (ii) the angular momentum brought in by the α-particle into the 208 Pb target (thin solid line). The thick arrow denotes L cr for 8 Be fusion as an inert projectile, whereas the thin arrow is l cr for fusion of a direct beam of α-particles (with half of E 0 ) on 208 Pb. The relevant features are:(1) the icf mechanism increases the angular momentum brought in by the α-particle into 208 Pb (thin solid line goes beyond the thin arrow), (2) the relative localization of the L 0 -window for cf and icf significantly depends on E 0 (dashed and thick solid lines), strongly overlapping as E 0 decreases, and (3) cf following 8 Be breakup becomes substantial with increasing E 0 , shifting down the maximum of the cf spin distribution (dashed line) with respect to the thick arrow.

Fig. 4

 4 shows the excitation energy distribution of the primary icf product 212 Po

Figure 4 .

 4 Figure 4. Excitation energy distribution of the primary icf product 212 Po resulting from the 8 Be + 208 Pb reaction at two laboratory energies E 0 : (a) 65 MeV, and (b) 45MeV. The arrow denotes the excitation energy of 212 Po, when this compound nucleus is formed by fusion of a direct beam of α-particles, with half of E 0 , on 208 Pb. Following breakup of8 Be, the fusing α-particle is significantly affected by the interaction with the surviving α-particle. This spreads the 212 Po excitation energy over a range of values. The distribution becomes narrower and its maximum approaches the arrow, as E 0 decreases towards the barrier.

Acknowledgments

The author thanks B. Bayman and M. Dasgupta for suggesting Appendix A, and J.A. Tostevin for discussions related to Appendix B. Support from the UK Science and Technology Facilities Council (STFC) Grant No. ST/F012012/1 is acknowledged.

Appendix A. Breakup probability function

Let us define two probabilities: (i) the probability of breakup between R and R + dR, ρ(R)dR [being ρ(R) a density of probability], and (ii) the probability of the weaklybound projectile's survival from ∞ to R, S(R). The survival probability at R + dR, S(R + dR), can be written as follows

Expression (A.1) suggests the following differential equation for the survival probability S(R),

whose solution is [S(∞) = 1]:

From (A.4), identifying ρ(R) with P L BU (R), we obtain expression (1) for the breakup probability integrated along a given classical orbit.

Appendix B. Matching prior-and post-breakup stages

The integrals of motion in the overall center-of-mass (cm) system are the total energy E tot = m T (m T +m P ) E 0 , the total linear momentum P tot = 0, and the total angular momentum L tot = m P b 0 ( v -V CM ) that is orthogonal to the initial reaction plane. m T , m P , b 0 , v, and V CM are the mass of the target and projectile, the impact parameter between the projectile and the target, the velocity of the incident projectile in the laboratory system and the cm velocity, respectively.

Just after breakup, the two-body projectile is excited to a definite state (ε 12 , 12 and d 12 ), as explained in Section 2. The relative vector between P and T ( R P T ) is also known. Thus, the separation between the three bodies is known. The modulus of the velocity between P and T (V P T = P P T /µ P T ) results from the total energy conservation

where U is the interaction potential between the target and the breakup fragments. The total linear momentum P tot = p T + p 1 + p 2 = p T + p P * , where p P * is the momentum of the center of mass of excited P relative to the overall cm. We need the velocities of P and T relative to the overall cm ( v P and v T ) to complete the initial conditions for subsequent propagation in time of the three bodies. These velocities are related to each other by the expressions

2)

where the magnitude of V P T is known through expression (B.1). To know the direction of this velocity the conservation of total angular momentum is applied. The total angular momentum L tot = 12 + L P T , so the angular momentum ( L P T ) associated with relative motion of P and T about cm is known. This vector can be written as

(B.4)

We now write v P in terms of radial and transverse components as follows:

P q, (B.5)

where r = R P T /R P T and q = n × r, being n = L P T /L P T . The transverse component v (q) P = L P T /(m P R P T ), and for the target v (q) T = -L P T /(m T R P T ). The radial component is obtained using expressions (B.2)-(B.3) and knowing the transverse component:

Both positive and negative roots are consistent with the conservation of the integrals of motion. Hence, both roots are uniformly sampled. Finally, the position and velocity vectors of the projectile fragments and the target are transformed to the laboratory system using Galilean transformations.


Igor N Mishustin

Claudio Ebel

Walter Greiner

Strong electric fields induced on a sharp stellar boundary

Keywords: numbers: 97.60.-s, 12.38.Mh, 04.70-s

Due to a first order phase transition, a compact star may have a discontinuous distribution of baryon as well as electric charge densities, as e.g. at the surface of a strange quark star. The induced separation of positive and negative charges may lead to generation of supercritical electric fields in the vicinity of such a discontinuity. We study this effect within a relativistic Thomas-Fermi approximation and demonstrate that the strength of the electric field depends strongly on the degree of sharpness of the surface. The influence of strong electric fields on the stability of compact stars is discussed. It is demonstrated that stable configurations appear only when the counter-pressure of degenerate fermions is taken into consideration.

Physics of compact stars is under intensive investigation for many decades. Modern studies are focused on the equation of state (EOS) of dense baryonic matter including possible phase transitions to the quark matter. One interesting feature of a strong first order phase transition is that it can generate a sharp discontinuity in the baryon density as a function of radial coordinate, see e.g. ref. [START_REF] Battachariya | Deconfinement Phase Transition in Compact Stars: Maxwell vs Gibbs Constraction of the Mixed phase[END_REF]. In the extreme situation when the EOS has a zero pressure point at a finite baryon density ρ c , as is the case e.g. for the MIT bag and NJL models [2,3], the discontinuity in the baryon density occurs at the surface of the star, i.e. ρ(R) = ρ c . Such a situation is expected in compact stars made of Strange Quark Matter (SQM), as first predicted in refs. [4,5]. If the matter is composed of several species with opposite electric charges and different masses, the presence of a sharp discontinuity should lead to a charge separation and generation of an electric field. This effect was first discussed in ref. [6] in context of the necked SQM stars, where the electrons from the bulk SQM matter can penetrate in vacuum through the sharp star's surface. As demonstrated in ref. [7], supercritical electric fields can be generated at the boundary of a nuclear core surrounded by an electronic cloud. In ref. [8] the structure of compact stars with a net charge at the surface was calculated within the General Relativity. The boundary effect may be important also for understanding the structure of a mixed phase, where domains of two phases have opposite electric charges [9,10]. The electrostatic interactions are important for the description of neutron-star crusts where atomic nuclei are embedded in a dense electron gas [11,12].

Our goal in the present paper is to study the conditions leading to the generation of supercritical fields at a star boundary. We solve a simplified problem replacing the spherical star boundary by the planar one. We expect this approximation to work very well for a compact star with radius of about 10 km. We assume further that the net positive charge, associated with protons or quarks, has a smooth-step distribution of the Woods-Saxon type

ρ p (z) = ρ p0 1 + exp z -z 0 a -1 , (1)
where ρ p0 is the positive charge density in the bulk (z → -∞) and a is a diffuseness parameter. The star boundary is located at z = z 0 and ρ p (z 0 ) = ρ p0 /2. At a → 0 this distribution approaches a rectangular step considered in ref. [7].

For a macroscopic object like a star the condition of global charge neutrality must hold to a very high precision, see discussion in ref. [START_REF] Glendenning | Compact Stars: Nuclear Physics, Particle physics and General Relativity[END_REF]. Therefore, the positive charge of Eq.

(1) must be fully neutralized by the negative charge of electrons. However, since electrons are light and interact only via the electromagnetic force, they will penetrate through the boundary and generate a local charge disbalance around the star surface. The induced electrostatic potential φ(z) is determined from the Poisson equation (see e.g. [START_REF] Landau | Quantum Mechanics -Non-relativistic Theory[END_REF])

d 2 φ dz 2 = -e [ρ p (z) -ρ e (z)] ≡ -eρ ch (z) , (2)
where e = √ 4πα=0.3028 is the proton charge and α is the fine-structure constant [26].

For a given proton distribution, Eq. ( 1), the electron charge distribution ρ e (z) should be determined self-consistently. For this purpose we use the Thomas-Fermi approximation [START_REF] Thomas | [END_REF]16], which should work well for an extended object like a heavy nucleus or star. The relativistic version of this method was considered e.g. in refs. [17,18]. In a semi-classical approximation the electron energy at point z can be written as

ǫ(k, z) = k 2 + m 2 e -V (z) , (3)
where k is its 3-momentum and -V (z) = -eφ(z) is the potential energy. At zero temperature all electronic states with k ≤ k F (z) are occupied, where k F (z) is the local Fermi momentum. It is determined from the condition

ǫ(k F (z), z) = k 2 F (z) + m 2 e -V (z) = µ = const , (4)
that gives

k F (z) = [µ + V (z)] 2 -m 2 e . (5)
The local electron density is found by integration over k:

ρ e (z) = 2 k F (z) 0 d 3 k (2π) 3 = (µ + V (z)) 2 -m 2 e 3/2 3π 2 . ( 6
)
By inserting this expression into Eq. ( 2) we obtain a non-linear differential equation for φ(z). To find its solutions we need to specify boundary conditions.

According to our assumption, the boundary is in direct contact with the vacuum and therefore the electron density as well as the electric potential must vanish at z → ∞. This condition can be fulfilled if µ ≤ m e . For µ < m e the electrons are bound to the surface, i. e.

a finite energy is needed to extract an electron from the star. This energy is well known in ordinary metals as the exit work. In this case the electron density should vanish at a certain distance from the surface, z * , which is found from the condition m e -V (z * ) = µ, see Eq. ( 4). The value of µ is determined from the condition of global charge neutrality:

e ∞ -∞ [ρ p (z) -ρ e (z)] dz ≡ σ + + σ -= 0 , (7)
where σ + and σ -are the areal densities of positive and negative charges accumulated at z < z 0 and z > z 0 respectively. Since calculated µ values are very close to m e , we use µ = m e for our estimates.

Another boundary condition follows from the requirement that far inside the star the electron charge density must completely compensate the positive charge density, i.e.

ρ ch (z → -∞) = 0, that means that k F e (z → -∞) ≡ k F 0 = (3π 2 ρ p0 ) 1/3 and V (z → -∞) = k 2 F 0 + m 2 e -µ = (3π 2 ρ p0 ) 2/3 + m 2 e -µ . (8)
These conditions are sufficient to determine the electron density distribution ρ e (z). Unfortunately, no analytical solution can be found without further approximations, and therefore numerical methods must be applied.

We have performed numerical calculations for different values of the diffuseness parameter a ranging from 1 fm, typical for nuclear surfaces, up to 1 nm, typical for solid-state surfaces.

As an example, in Fig. 1 we present results for a=5 fm which is rather close to a "sharp surface" limit. The upper, middle and lower panels show, respectively, the electron density, electric potential and electric field calculated for the given proton density (dashed line in the upper panel). One can clearly see the deviation of the electron density from the proton one in the layer of ±50 fm around the surface. This gives rise to a sharp fall of the potential, and as a consequence, to a very strong electric field generated in this layer. The maximum of the field, about 1.5 MV/fm, is reached at z = z 0 . This value is almost 600 times higher than the critical electric field for spontaneous electron-positron pair production [START_REF] Greiner | Quantum Electrodynamics of Strong Fields[END_REF], E c = 2m 2 c 3 /e = 0.0026 MV/fm=2.6 • 10 16 V/cm. However, in the considered system the pair production does not happen because all electron states with ǫ e (k, z) < m e are occupied and therefore Pauliblocked.

This situation is clarified in Fig. 2, where the electron levels in a strong electric potential, -V (z), are shown schematically. As explained above, the negative potential inside the medium (z < z 0 ) is generated by the charge separation at the surface. In the vacuum such production of negative muons and pions, which may replace the electrons. They become possible in a strong enough field when V 0 > m µ , m π . This opens the possibility of pion condensation induced by strong electric fields [17].

With increasing a the electron density comes closer and closer to the proton density and the charge separation at the surface diminishes. The electric potential difference at z → ±∞ is fixed by the boundary conditions, but its gradient, i. e. the electric field, across the transition region becomes smaller and smaller. We have checked that the electric field falls below the critical value at a ≈ 10 pm. In Table 1 we give the maximum value of the electric field, E max , calculated for several values of the diffuseness parameter a. E max can be calculated by using the formula which follows from the Gauss law:

E max ≡ E(z = z 0 ) = e z 0 -∞ ρ ch (z)dz = σ + , (9)
where σ + is the areal density of the net positive charge at z < z 0 .

For our discussion below we need to estimate the electric field and other related quantities on a qualitative level. This can be done by using the ultrarelativistic limit, V ≫ m e , which is justified everywhere except of the asymptotic region at z → ∞. Taking for simplicity the limit a → 0, we get the exact solution of the Poisson equation (2)[27]

V (z) =      V 0 1 -c 1 exp z -z 0 λ , z < z 0 V 0 c 2 + (z -z 0 ) √ 6λ -1 , z > z 0 (10)
where

λ = π eV 0 ≈ π ek F 0 . ( 11
)
The constants c 1 = 0.2374 and c 2 = 1/(1 -c 1 ) = 1.3113 are obtained from the continuity conditions for V (z) and dV /dz at z = z 0 . It is interesting that the decay length λ is determined by the electron Fermi momentum in the bulk of the medium k F 0 . One can see that the potential changes much faster at z < z 0 than at z > z 0 . The electric field can be easily obtained now by differentiating Eq. ( 10),

E(z) =      c 1 V 0 eλ exp z -z 0 λ , z < z 0 V 0 √ 6eλ c 2 + (z -z 0 ) √ 6λ -2 z > z 0 (12)
It is obvious that the maximum of the electric field is reached at z = z 0 :

E max = c 1 V 0 eλ = c 1 V 2 0 π = c 1 k 2 F 0 π . ( 13
)
The areal charge density at z = z 0 is given by eq. ( 9). Assuming that this charge is distributed over the layer of width λ we can estimate the average charge density in this layer

eρ ch ≈ σ + λ = c 1 ek 3 F 0 π 2 , (14)
i. e. about 70% of the proton density at z → -∞. The energy density of the electric field is E em = E 2 /2, and therefore the total electrostatic energy is easily calculated as

E em = ∞ -∞ E em dV = c 1 (2 -c 1 ) 12
V 2 0 e 2 λ S ≈ c 1 (2 -c 1 ) πS 4e ρ p , ( 15
)
where S is the surface area and ρ p is the bulk proton density. Note that this expression is approximately equal to the energy stored in a planar capacitor with areal charge density σ + and gap λ.

The authors of ref. [7] discuss the possibility of a new family of stable star-like objects (nuclear cores) where the repulsive electromagnetic force, associated with the charged layer at a sharp boundary, is balanced by the gravitational force. Here we consider this possibility taking into account additionally the kinetic pressure of degenerate fermions. First of all we point out that in a macroscopic nuclear matter, besides neutrons and protons, the electrons must be present to neutralize the bulk proton charge. By this reason the β-equilibrium will be shifted to the neutron-rich side, so that the matter will contain mainly neutrons with a small admixture of protons and electrons. According to standard calculations of β-equilibrated npe-matter in neutron stars (see e.g. refs. [11,[START_REF] Glendenning | Compact Stars: Nuclear Physics, Particle physics and General Relativity[END_REF]), the fraction of protons, ξ = ρ p /ρ B should be on the level of a few percent. As well known, such matter is unbound, in contrast to the symmetric nuclear matter with N n ≈ N p without electrons.

For simplicity we replace the baryon, proton and neutron densities by their mean values:

ρ B = 3N B 4πR 3 , ρ p = k 3 F p 3π 2 = ξρ B , ρ n = k 3 F n 3π 2 = (1 -ξ)ρ B , (16)
where R is the star radius. Now let us consider a star-like configuration with sharp boundary at r = R, where we expect formation of a strong electric field. In our qualitative discussion we represent the star energy per baryon as the sum of three contributions:

W (R) = W kin (R) + W em (R) + W grav (R) , (17)
where terms in the r.h.s. represent the kinetic energy of degenerate fermions (npe), the electrostatic energy and the gravitational energy. The kinetic energy term is calculated easily for any fermionic species of mass m i and Fermi momentum k F i :

W (i) kin (x i ) = 3m i 8x 3 i 2x 3 i + x i 1 + x 2 i -arcsinh (x i ) =    m i 1 + 3 10 x 2 i , x i ≪ 1 3 4 x i , x i ≫ 1 (18)
where i=n,p,e and x i = k F i /m i . Below we include in W kin only neutron contribution and express k F n as a function of R using Eq. ( 16),

k F n = 9π 4 1/3 N 1/3 n R . (19)
The electrostatic energy of a spherical dipole layer can be obtained from Eq. ( 15) by replacing ρ p = 3ξN B /4πR 3 and S = 4πR 2 . Then the energy per baryon is

W em (R) = 3π 4 c 1 (2 -c 1 ) ξ eR ≡ c 3 ξ eR , (20)
where c 3 ≈ 0.426. It should be emphasized that this energy is by factor λ R ∼ 1/N 1/3 p smaller than the energy of net charge Q = 4πR 2 σ + distributed in a spherical layer of wdth λ and radius R (see ref. [8]). Finally, the gravitational energy of a baryon at the surface of a star of mass M = N B m B and radius R is easily calculated in Newtonian gravity:

W grav (R) = - GMm B R = - m B m Pl 2 N B R , (21)
where

m Pl = 1/ √ G ≈ 1.22 • 10 19
GeV is the Planck mass and m B = 939 MeV is the nucleon mass.

Putting together all the contributions we can now rewrite Eq. (17) as

W (R) = 3 4 9π 4 1/3 N 1/3 n R + c 3 ξ e - m B m Pl 2 N B 1 R , (22)
where the ultrarelativistic limit of Eq. ( 18) has been used in the kinetic term. One can immediately see that the electrostatic and gravitational contributions become equal at

N min = c 3 ξ e m Pl m B 2 ≈ 6.26 • 10 36 , (23)
where the value ξ = 1/38 was used for the estimation. The corresponding mass M min = N min m B is about 10 13 g, i. e. almost 20 orders of magnitude smaller than the maximum mass of a neutron star. At N B < N min there exists no equilibrium configurations because the kinetic pressure will force the system to expand indefinitely.

At N B > N min the negative gravitational contribution exceeds the electrostatic one.

The latter becomes completely negligible and can be neglected at baryon numbers ∼ 10 57 characteristic for compact stars. As was first demonstrated by Chandraseckhar [22] and Landau [23], a compact star can be stabilized by the pressure of a degenerate fermion gas alone. From ( 22) one can see that the balance between kinetic and gravitational energy can be achieved if the total number of baryons is smaller than

N max = 3 4 3/2 9π 4 1/2 m Pl m n 3 ≈ 3.78 • 10 57 . (24)
As follows from this simple calculations, the corresponding maximum mass of a neutron star, M max = N max m n is about 3.2 M ⊙ . Numerical calculations within General Relativity [24] for the degenerate noninteracting neutron gas predict somewhat lower maximum mass, M max = 0.7M ⊙ . The repulsive neutron-neutron interaction may increase this value up to 4M ⊙ (see e.g. ref. [25]). One can see that stars with N min < N B < N max have a stable equilibrium state. Indeed, Eq. (22) shows that in this case W (R) > 0 and changes as 1/R, so that the star should expand to lower densities. But this expansion will terminate at a large enough R, when the neutron gas becomes non-relativistic. As follows from Eq. (18), in this limit the energy per particle is m n + 0.3k 2 F n /m n , i. e. changes as 1/R 2 at large R. Since the negative gravitational contribution still changes as 1/R, it wins so that the total energy per baryon approaches m n from below. Therefore, a minimum develops in W (R) allowing for a stable configuration. This argument was presented by Landau in ref. [23]. In ref. [25] it was reformulated for the case of General Relativity and for an arbitrary fermion mass.

It is clear now that ignoring the kinetic-energy term in Eq. ( 22), as done in ref. [7], will lead to the erroneous conclusion regarding the star stability. Although the electrostatic and gravitational forces can balance each other at N B = N min , this state corresponds to a maximum of energy. As one can easily see, the system will either expand or collapse depending on its actual baryon number compared to N min . Moreover, the electrostatic force can support only the surface layer of matter where the electric field is nonzero. The inner layers of the star require the kinetic pressure to resist gravity. Therefore, the mechanism discussed in ref. [7] cannot produce any new stable object. Inclusion of the electrostatic effects will only slightly increase the star mass, see e. g. ref. [8].

In conclusion, we have demonstrated that very strong electric fields can be generated on a sharp boundary separating the bulk matter from the vacuum. The strength of the field is determined by the separation of positive and negative charges in the surface layer. For the proton charge distribution with the diffuseness parameter of a nuclear scale the field strength may exceed significantly the critical value for the electron-positron pair production. The field strength diminishes when the charge distribution becomes more and more smooth. We have also considered star-like configurations with a strong electric field on the boundary.

Our simple analysis shows that such electric field alone cannot balance the gravitational force to produce stable objects. The stable compact star configurations appear only when the counter-pressure of degenerate fermions is taken into consideration. For a more accurate description of compact stars with strong electric fields one should perform calculations within General Relativity and include effects of strong interactions in the energy of the nuclear core.

FIG. 1 :FIG. 2 :

 12 FIG.1:The upper, middle and lower panels show, respectively, the electron density (full line), electric potential and electric field, calculated for the given proton density profile (dashed line in the upper panel) as given by Eq.(1) with ρ p0 = 0.01 fm -3 and a=5 fm.

TABLE I :

 I Maximum electric field and areal charge density calculated for several values of the diffuseness parameter a. a E max [ MV /fm] Emax /Ec σ + e/f m 2 1 fm 1.82 • 10 -0 6.89 • 10 +2 1.01 • 10 -1 10 fm 1.16 • 10 -0 4.38 • 10 +2 6.41 • 10 -2 100 fm 2.05 • 10 -1 7.74 • 10 +1 1.13 • 10 -2 1 pm 2.08 • 10 -2 7.86 • 10 0 1.18 • 10 -3 10 pm 2.08 • 10 -3 7.86 • 10 -1 1.18 • 10 -4 100 pm 2.08 • 10 -4 7.86 • 10 -2 1.18 • 10 -5 1 nm 2.08 • 10 -5 7.87 • 10 -3 1.25 • 10 -6

We are grateful to Prof. Remo Ruffini who attracted our attention to the problems discussed in this paper. We acknowledge his kind hospitality during our visit to the International Center for Relativistic Astrophysics (ICRA) at Pescara in March 2009. We also thank Dr. Thomas Bürvenich and members of the ICRA group for fruitful discussions. This work was supported in part by the DFG Grant 436 RUS 113/957/0-1 (Germany), and by the Grants NS-3004.2008.2 and RFBR 09-02-91331 (Russia).