Solution of Hartree-Fock-Bogoliubov equations and fitting procedure using N2LO Skyrme pseudo-potential in spherical symmetry
aa r X i v : . [ nu c l - t h ] J u l Solution of Hartree-Fock-Bogoliubov equations and fitting procedure using N2LOSkyrme pseudo-potential in spherical symmetry.
P. Becker, ∗ D. Davesne, † J. Meyer, J. Navarro, ‡ and A. Pastore § Universit´e de Lyon, Universit´e Lyon 1, 43 Bd. du 11 Novembre 1918, F-69622 Villeurbanne cedex, FranceCNRS-IN2P3, UMR 5822, Institut de Physique Nucl´eaire de Lyon IFIC (CSIC-Universidad de Valencia), Apartado Postal 22085, E-46.071-Valencia, Spain Department of Physics, University of York, Heslington, York, Y010 5DD, United Kingdom (Dated: October 22, 2018)We present the development of the extended Skyrme N2LO pseudo-potential in the case of spher-ical even-even nuclei calculations. The energy density functional is first presented. Then we derivethe mean-field equations and discuss the numerical method used to solve the resulting fourth-orderdifferential equation together with the behaviour of the solutions at the origin. Finally, a fittingprocedure for such a N2LO interaction is discussed and we provide a first parametrization. Typicalground-state observables are calculated and compared against experimental data.
PACS numbers: 21.30.Fe 21.60.Jz
I. INTRODUCTION
The Nuclear Energy Density Functional (NEDF) theory allows us to describe properties of nuclei from light toheavy nuclei and from drip-line to drip-line [1]. Several functionals have been developed in the recent years, but themost widely used [2, 3] are those derived from the non-relativistic zero-range Skyrme interaction [4]. Since its firstapplications to atomic nuclei [5], this interaction has proven to be very well suited to describe nuclear observables atvery reduced computational cost [6].A crucial aspect in building a functional is to determine the values of its coupling constants. Despite its apparentsimplicity, this is a very delicate aspect: a badly determined coupling constant can give rise to unphysical instabili-ties [7–13] and thus to unphysical results. A possibility for avoiding them is to find an adequate set of observables sothat all coupling constants are properly constrained during the optimization procedure [14, 15]. In Ref. [16], we havepresented an alternative solution to avoid unphysical instabilities based on the linear response (LR) formalism in infi-nite nuclear medium. This solution is particularly simple and very efficient especially for some particular terms of thefunctional that are odd under time reversal symmetry and give very little contribution to masses of odd-systems [8].However, avoiding unphysical instabilities is not the only requirement to have an effectient functional : one also hasto check how it performs to describe nuclear observables. On this point, the UNEDF collaboration [17] has recentlystudied much in detail the properties of Skyrme functionals against a large set of nuclear observables [18–20]. Themain conclusion in their last article [20] is that the standard Skyrme functional [2] has reached its limits. If we wantto improve the description of experimental data (as masses, radii, fission barriers,...) we need to follow two paths:explore different functional forms or develop functionals at multi-reference level [21].Following the idea of Carlsson and collaborators [3, 22], we have decided to explore the first path and to study theimpact of additional gradient terms into the Skyrme pseudo-potential [23]. The gradient terms have been introducedin a systematic way by considering all possible combinations allowed by the symmetries of the problem up to 6thpower. The resulting pseudo-potential has been called N ℓ LO which by definition incorporates gradients up to order2 ℓ . Within this language, the standard Skyrme interaction [24] is named N1LO. In Ref. [25], we have shown theexplicit connection between the Taylor momentum expansion of any finite range interaction and the actual form ofthe N ℓ LO pseudo-potential [3]. In that article, we have also proven that such an expansion works fairly well in infinitenuclear medium and that the main properties of the Equation of State (EoS) of a finite-range interaction can be fairlyreproduced by truncating the momentum expansion to fourth order (N2LO). The result is coherent with previousfindings based on Density Matrix Expansion (DME) [26]: the role of fourth order terms is important and it leads to ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: navarro@ific.uv.es § Electronic address: [email protected] a remarkable improvement of the DME results when compared to finite-range interactions. Higher order terms canthus be neglected as a first step since their contribution becomes systematically less important.At present, the only existing parametrizations of the extended Skyrme N2LO/N3LO pseudo-potentials have beenobtained by considering only properties of infinite nuclear medium [27, 28], that is without taking into accountproperties of finite nuclei. In order to remedy this aspect, we present here a new Skyrme Hartree-Fock-Bogoliubov(HFB) code that incorporates higher order derivatives terms appearing in N2LO. It is worth remainding at thispoint that an alternative code named HOSPHE [29] has already been published. This code, based on Harmonic-Oscillator (HO) basis also considers the most general functional form of the N3LO functional [22] using sphericalbasis representation. However, following our previous findings of Ref. [23], we have decided to express the N ℓ LOpseudo-potential in Cartesian coordinates and to develop for this specific case a numerical code to work in coordinatespace: the r-space representation is in fact more convenient to be used in a fitting procedure since we do not need touse a very large number of basis states to achieve convergence. See Ref. [30] for more details.The article is organized as follows: in Sec. II we present the general functional formalism for the N2LO pseudo-potential and in Sec. III we specialize the formalism for the spherically symmetric case. In Sec. IV we present in detailthe generalization of the Hartree-Fock-Bogoliubov equations to include the N2LO pseudo-potential. In Sec. V wepresent the fitting protocol to determine the parameters of the new N2LO functionals. Finally we give our conclusionsin Sec. VI.
II. N2LO SKYRME FUNCTIONAL
The N2LO Skyrme pseudo-potential as described in Refs. [3, 22] is a generalization of the standard Skyrme inter-action, corresponding to the expansion of the momentum space matrix elements of a generic interaction in powers ofthe relative momenta k , k ′ up to the fourth order. Following [31], the form considered in this article respects bothGalilean and local gauge invariance [32]. It is written as the sum of three terms V N2LO = V C N2LO + V LS N1LO + V DD N1LO . (1)The central term reads V C N2LO = t (1 + x P σ )+ 12 t (1 + x P σ )( k + k ′ )+ t (1 + x P σ )( k · k ′ )+ 14 t (4)1 (1 + x (4)1 P σ ) h ( k + k ′ ) + 4( k ′ · k ) i + t (4)2 (1 + x (4)2 P σ )( k ′ · k )( k + k ′ ) . (2)In these expressions, a Dirac function δ ( r − r ) is to be understood, but has been omitted for the sake of clarity.See Ref. [1] for details on the adopted notations. The spin-orbit term V LS N1LO is not affected by the inclusion of higherorder gradient terms: in Ref. [25], we have shown that other possible spin-orbit terms are suppressed once the localgauge invariance [3, 33] is imposed. In Ref. [25], we have discussed in details the problem of local gauge invariancefor spin-orbit term and in particular the possible violation of such a symmetry for finite-range spin-orbit terms. Thedensity-dependent term V DD N1LO has also exactly the same structure as in the standard Skyrme interaction [24], sinceits nature is to mimic the effect of a three-body term [5, 34]. Tensor terms should be also included into Eq. (2).In Ref. [27], we have discussed them based on the partial-wave decomposition of the total EOS. In finite nuclei it isactually very difficult to constrain them in NEDF [35] because of their strong competition with the spin-orbit termin modifying the underlying single-particle structure [36]. For this preliminary exploration, we have thus decided toneglect them. Finally, it is worth mentioning that in the present article we will always use the complete interaction inthe sense that we will not discard the so-called J tensor terms [36] as often done in the literature. For the Coulombinteraction between protons, we adopt the same procedure as described in Ref. [24] i.e. using the standard Slaterapproximation for the exchange term [37].Starting from Eq. (2), it is possible to derive the explicit form of the Skyrme functional in Cartesian coordinates.We write it as E = X t E (1) , even t + E (1) , odd t + E (2) , even t + E (2) , odd t , (3)where t = 0 , ℓ LO terms E ( ℓ =1 , . The standard terms E (1) t read [36] E (1) , even t = C ρt [ ρ ] ρ t + C ∆ ρt ρ t ∆ ρ t + C τt ρ t τ t − C Tt z X µ,ν = x J t,µν J t,µν + C ∇ Jt ρ t ∇ · J t , (4) E (1) , odd t = C st [ ρ ] s t − C τt j t + C ∆ st s t · ∆ s t + C Tt s t · T t + C ∇ Jt s t · ∇ × j t , (5)while the new terms can be written as E (2),even t = C (∆ ρ ) t (∆ ρ t ) + C Mρt M Mρ, even t + C Mst M Ms, even t , (6) E (2),odd t = C (∆ s ) t (∆ s t ) + C Mρt M Mρ, odd t + C Mst M Ms, odd t , (7)where M Mρ, even = (cid:8) ρ Q + τ (cid:9) + 2 [Re( τ µν )Re( τ µν ) − Re( τ µν ) ∇ µ ∇ ν ρ ] , (8) M Ms, even = n ( ∇ µ J µν ) + 4 J µν V µν − Im( K µνκ )Im( K µνκ ) o , (9) M Mρ, odd = n ( ∇ · j ) + 4 j · Π o , (10) M Ms, odd = (cid:8) s · S + T (cid:9) + 2 [Re( K µνκ )Re( K µνκ ) − Im( τ µν )Im( τ µν ) − Re( K µνκ ) ∇ µ ∇ ν s κ ] . (11)These terms contain six new densities: τ µν , V µν , Π , K µνκ , Q and S . Their explicit definition is given in Appendix B. III. N2LO FUNCTIONAL IN SPHERICAL SYMMETRY
In the present section, we limit ourselves to the case of spherical symmetry. In this case, the single-particle wavefunction can be written as follows ψ nℓjmq ( r ) = 1 r R nℓjq ( r ) Ω ℓjm (ˆ r ) , (12)where n is the principal quantum number, Ω ℓjm (ˆ r ) is a solid spherical harmonic [38] and ℓjm refer respectively to theorbital angular momentum, the total angular momentum and its relative projection along the z -axis. Here q ≡ n, p stands for proton (p) or neutron (n). In our formalism the two nuclear species are not mixed explicitly [2, 39]. Byconsidering only even-even systems, we can further simplify the expressions given in Eq. (3) E (1) = C ρ ρ + C ρ ρ + C ∆ ρ ρ ∆ ρ + C ∆ ρ ρ ∆ ρ (13)+ C τ ρ τ + C τ ρ τ − C T J − C T J + C ∇ J ρ ∇ · J + C ∇ J ρ ∇ · J ; E (2) = C (4)∆ ρ (∆ ρ ) + C (4)∆ ρ (∆ ρ ) + C (4) Mρ n (cid:2) ρ Q + τ (cid:3) o + C (4) Mρ n [Re( τ ,µν )Re( τ ,µν ) − Re( τ ,µν ) ∇ µ ∇ ν ρ ] o + C (4) Mρ n (cid:2) ρ Q + τ (cid:3) o + C (4) Mρ n [Re( τ ,µν )Re( τ ,µν ) − Re( τ ,µν ) ∇ µ ∇ ν ρ ] o − C (4) Ms h ( ∇ µ J ,µν ) + 4 J ,µν V ,µν − Im( K ,µνκ )Im( K ,µνκ ) i − C (4) Ms h ( ∇ µ J ,µν ) + 4 J ,µν V ,µν − Im( K ,µνκ )Im( K ,µνκ ) i . (14) A. Local densities
Let us introduce the short-hand notation α = { nℓjq } and C α = j ( j + 1) − ℓ ( ℓ + 1) − . The explicit expressionsof the densities in spherical symmetry (we limit ourselves to systems that are even under time-reversal) up to secondorder take the form [40] ρ ( r ) = X α (2 j + 1)4 π R α ( r ) r , (15) τ ( r ) = X α (2 j + 1)4 πr "(cid:18) R ′ α ( r ) − R α ( r ) r (cid:19) + ℓ ( ℓ + 1) r R α ( r ) , (16) J ( r ) = X α (2 j + 1)4 π C α R α ( r ) r . (17) τ ( r ) can be conveniently decomposed in a radial and centrifugal part as τ = τ R, + τ C, where τ R, ( r ) = X α (2 j + 1)4 πr (cid:20) R ′ α ( r ) − R α ( r ) r (cid:21) , (18) τ C, ( r ) = X α (2 j + 1)4 π ℓ ( ℓ + 1) r R α ( r ) r . (19)Eq. (17) corresponds to the radial part of the J µν, ( r ) spin-orbit vector density defined as J µν, = 12 ǫ µνκ J κ = 12 ǫ µνκ X κ r J ( r ) , (20)where X µ represents the Cartesian coordinates. If we now come to fourth order, the explicit expressions of the newdensities in spherical symmetry take the form τ µν, ( r ) = 12 τ C, ( r ) δ µν + X µ X ν r (cid:20) τ R, ( r ) − τ C, ( r ) (cid:21) , (21) V ( r ) = X α (2 j + 1)4 πr C α (cid:20) R α r [ ℓ ( ℓ + 1) + 2] + R ′ α ( r ) r − R ′ α ( r ) R α ( r ) r (cid:21) , (22) Q ( r ) = X α (2 j + 1)4 πr (cid:20) R ′′ α ( r ) − ℓ ( ℓ + 1) R α ( r ) r (cid:21) , (23) K µνκ, ( r ) = i K1 ( r ) ǫ µνκ + i K2 (cid:20) ǫ µκM X M X ν r + ǫ µνM X M X κ r + ǫ κνM X M X µ r (cid:21) . (24)We have defined K and K asK1 ( r ) = X α (2 j + 1)16 πr C α R ′ α ( r ) R α ( r ) , (25)K2 ( r ) = X α (2 j + 1)16 πr C α (cid:20) r R α ( r ) − R ′ α ( r ) R α ( r ) (cid:21) . (26) τ µν, (r) is the kinetic density tensor. The usual N1LO τ (r) density is given by its trace X µ τ µµ, ( r ) = τ . (27)The even part of the N2LO functional only receives a non-vanishing contribution from the real part of this density(Eq. 14). Given that the imaginary part is zero under spherical symmetry, we will write τ µν, (r) instead of Re( τ µν, (r))in the following. Similarly to J ( r ), V ( r ) is the radial part of the vector density V µν, V µν, = 12 ǫ µνκ X κ r V ( r ) , (28)and it can be decomposed in a radial and centrifugal part as V = V R, + V C, where V R, ( r ) = X α (2 j + 1)4 πr C α (cid:20) R ′ α ( r ) − r R ′ α ( r ) R α ( r ) + 2 r R α ( r ) (cid:21) . (29) V C, ( r ) = X α (2 j + 1)4 πr C α (cid:20) ℓ ( ℓ + 1) r R α ( r ) (cid:21) . (30)Since the K µνκ, ( r ) density is imaginary in spherical symmetry, the N2LO functional (Eq. 14) only receives acontribution of this density multiplied by itself. As for the τ µν, (r), we will use K µνκ, ( r ) without mentioninganymore that it actually stands for the imaginary part of this density.Some additional expressions which represent the new contributions to the functional are also written below forcompleteness. τ µν, ( r ) τ µν, ( r ) = τ R, ( r ) − τ C, ( r ) (31) τ µν, ∇ µ ∇ ν ρ = ρ (2) τ R + ρ (1) r τ C (32) J µν, V µν, = 12 J ( r ) V ( r ) (33) K µνκ, K µνκ, = 6K1 ( r ) + 6K2 ( r ) − ( r )K2 ( r ) . (34)In order to have a qualitative and quantitative idea of all these densities, we represent in Fig. 1, the isoscalar densitiesin Pb. These densities have been determined using a single particle basis obtained from a fully-converged Hartree-Fock (HF) solution based on the SLy5 functional [24]. We observe that all the densities used here are well-behavedat the origin of the coordinate system. D en s i t i e s r [fm] ρ [fm -3 ] τ [fm -5 ]J[fm -4 ]Q[fm -7 ]V[fm -6 ]K1 [fm -5 ]K2 [fm -5 ] FIG. 1: (Colors online) Isoscalar densities in
Pb calculated using single particle wave functions obtained by a SLy5 mean-fieldsolution. See text for details.
IV. HARTREE-FOCK-BOGOLIUBOV EQUATIONS IN SPHERICAL SYMMETRY
In this section we describe the method used to solve the complete Hartree-Fock-Bogoliubov (HFB) equations andthe numerical tests we have performed.
A. Hartree-Fock
We start considering closed-shell nuclei for which the HFB equations can be safely reduced to the standard Hartree-Fock (HF) equations. They read [5, 41] h q ( r ) R nljq ( r ) = ε qnlj R nljq ( r ) , (35)where R nljq ( r ) is the radial part of the single-particle wave-function given in Eq. (12). The corresponding Hamiltonianis derived as a functional derivative as h q ( r ) = A q d dr + A q d dr + A q R d dr + A q R ddr + A q R + ℓ ( ℓ + 1) r (cid:20) A q C d dr + A q C ddr + ℓ ( ℓ + 1) r A q CC + A q C (cid:21) + W q R d dr + W q R ddr + W q R + ℓ ( ℓ + 1) r W q C . (36)We observe that the inclusion of 4th order term in the interaction translates into a fourth order differential equation.Although this is quite unusual in nuclear physics, a 4th order differential equation is routinely solved in other physicalsystems, as for example to describe the behaviour of a bending solid beam [42].The coefficients in Eq. (36) are defined as A q = C Mρ − ρ + 2 C Mρ ρ q , (37) A q = 2 C Mρ − ρ (1)0 + 4 C Mρ ρ (1) q , (38) A q R = − ~ m − C τ − ρ − C τ ρ q + C Mρ − h ρ (2)0 − τ R, − τ C, i + 2 C Mρ h ρ (2) q − τ R,q − τ C,q i , (39) A q C = − C Mρ − ρ − C Mρ ρ q , (40) A q R = − C τ − ρ (1)0 − C τ ρ (1) q + 2 C Mρ − h ρ (3)0 − τ (1) R, − τ (1) C, i + 4 C Mρ h ρ (3) q − τ (1) R,q − τ (1) C,q i , (41) A q C = 2 C Mρ − (cid:16) − ρ (1)0 + 2 ρ r (cid:17) + 4 C Mρ (cid:16) − ρ (1) q + 2 ρ q r (cid:17) , (42) A q R = U q ( r ) + C τ − ρ (1)0 r + 2 C τ ρ (1) q r , + 2 C Mρ − " τ (1) R, r + τ (1) C, r − ρ (3)0 r + 4 C Mρ " τ (1) R,q r + τ (1) C,q r − ρ (3) q r , (43) A q C = ~ m + C τ − ρ + 2 C τ ρ q + C Mρ − " τ R, + 4 τ C, + 2 ρ (1)0 r − ρ (2)0 − ρ r + 2 C Mρ " τ R,q + 4 τ C,q + 2 ρ (1) q r − ρ (2) q − ρ q r , (44) A q CC = C Mρ − ρ + 2 C Mρ ρ q . (45)Here we used the shorthand notation C x − = C x − C x with x = ρ, ∆ ρ, . . . . The exponent ( i = 1 , , ,
4) in the densitiesstands for the derivative order. Finally, the central field appearing in the previous equation reads U q ( r ) = 2 C ρ − ρ + 4 C ρ ρ q + 2 C ∆ ρ − ∆ ρ + 4 C ∆ ρ ∆ ρ q + C τ − τ + 2 C τ τ q + 2 C (∆ ρ ) − ∆∆ ρ + 4 C (∆ ρ ) ∆∆ ρ q C Mρ − [ Q − ∇ µ ∇ ν τ µν, ] + 2 C Mρ [ Q q − ∇ µ ∇ ν τ µν,q ]+ C ∇ J − ∇ · J + 2 C ∇ J ∇ · J q . (46)This field is obtained through the variational principle varying the matter density ρ , and it receives contributionsfrom both N1LO and N2LO terms. In Fig. 2 we show the coefficients A qR and the central field U q obtained with afully converged HF calculation (cf Tab. V) in Pb using a N2LO pseudo-potential. We refer the reader to Sec. Vfor more details on this parametrisation. On the same figure we also report the corresponding values obtained withSLy5. As it should be, SLy5 induces non-zero contributions only for the terms originating from the N1LO part of the -0.35-0.3-0.25-0.2-0.15-0.1-0.05 0 0 2 4 6 8 10 12 A q4 (r) [ M e V f m ] r [fm] SN2LO1 protonsneutrons -0.1-0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0 2 4 6 8 10 12 A (r) [ M e V f m ] r [fm]SN2LO1 protonsneutrons -31-30-29-28-27-26-25-24-23-22-21-20 0 2 4 6 8 10 12 A R q (r) [ M e V f m ] r [fm] SN2LO1 protonsneutronsSLy5 protonsneutrons -2-1 0 1 2 3 4 5 0 2 4 6 8 10 12 A R q (r) [ M e V f m ] r [fm] SN2LO1 protonsneutronsSLy5 protonsneutrons -2-1.5-1-0.5 0 0.5 1 1.5 2 2.5 0 2 4 6 8 10 12 A R q (r) [ M e V ] r [fm] SN2LO1 protonsneutronsSLy5 protonsneutrons -90-80-70-60-50-40-30-20-10 0 0 2 4 6 8 10 12 U q (r) [ M e V ] r [fm] SN2LO1 protonsneutronsSLy5 protonsneutrons
FIG. 2: (Colors online) Radial dependence of the coefficients defined in Eq. (36) for
Pb obtained using the SN2LO1 andSLy5 interactions. See text for details. functional. In Fig. 3 we show the other set of fields appearing in Eq. (36) and corresponding to the centrifugal parts.These fields are active only for non-zero orbital momentum states. All the fields behave normally around r = 0 apartfrom the A q c , A q c that present a divergency. Such a behaviour, which already exists at N1LO level for the centrifugalfield, is actually not a problem as we will see in Sec. IV B when we examine the asymptotic properties of our 4thorder differential equation. We will then demonstrate that there exists a particular solution of Eq. (36) that exhibitsno divergency. Although we have only one explicit spin-orbit term in the effective interaction, we obtain four distinct A C q (r) [ M e V f m ] r [fm] SN2LO1 protonsneutrons -50-45-40-35-30-25-20-15-10-5 0 0 2 4 6 8 10 12 A C q (r) [ M e V f m ] r [fm] SN2LO1 protonsneutrons A C q (r) [ M e V f m ] r [fm] SN2LO1 protonsneutronsSLy5 protonsneutrons A CC q (r) [ M e V f m ] r [fm] SN2LO1 protonsneutrons
FIG. 3: (Colors online) Same as Fig. 2, but for centrifugal fields given in Eq. (36). contributions to the mean-field equation W q R ( r ) = − C α " C T − J r + 2 C T J q r + C ∇ J − ρ (1)0 r + 2 C ∇ J ρ (1) q r (47)+ C α " C (4) Ms − J r − J (1)0 r − V ( r ) r + 2 K ( r ) r ! + 4 C (4) Ms J q r − J (1) q r − V q ( r ) r + 2 K q ( r ) r ! ,W q C ( r ) = C α (cid:20) − C Ms − J ( r ) r − C Ms J q ( r ) r (cid:21) , (48) W q R ( r ) = C α " C Ms − J (1)0 ( r ) r − J ( r ) r ! + 4 C Ms J (1) q ( r ) r − J q ( r ) r ! , (49) W q R ( r ) = C α (cid:20) C Ms − J ( r ) r + 4 C Ms J q ( r ) r (cid:21) . (50)This is a very interesting feature of our functional which appears to have more flexibility than N1LO. This newdependence could be of particular interest in different situations, by instance in adjusting centroids of single particlestates without the need of using an explicit tensor term. Moreover, these terms are associated with the first twoderivatives in the differential equation, contrary to the standard Skyrme interaction, and one of them is a centrifugalterm. Such a term could thus allow to act on the single-particle levels with a new dependency in l . It is worthmentioning that several Skyrme functionals use different coupling constants in the spin-orbit sector to enrich thefreedom of the corresponding field [43]. In such a case, the link with the underlying interaction is then broken. Thenew N2LO functional presented here has the advantage of keeping such a link and also gaining a more complexspin-orbit structure, thus making it a suitable candidate for multi-reference calculations. In Fig. 4, we show thedifferent spin-orbit contributions. The current parametrisation SN2LO1 leads to relative small values, but we shouldnot exclude a priori the possibility of finding significative corrections with a different set of parameters. -0.035-0.03-0.025-0.02-0.015-0.01-0.005 0 0.005 0.01 0 2 4 6 8 10 12 W R q (r) [ M e V f m ] r [fm] SN2LO1 protonsneutrons -0.03-0.025-0.02-0.015-0.01-0.005 0 0.005 0.01 0.015 0.02 0 2 4 6 8 10 12 W R q (r) [ M e V f m ] r [fm] SN2LO1 protonsneutrons -4-3-2-1 0 1 2 3 0 2 4 6 8 10 12 W R q (r) [ M e V ] r [fm] SN2LO1 protonsneutronsSLy5 protonsneutrons -0.01-0.005 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0 2 4 6 8 10 12 W C q (r) [ M e V f m ] r [fm] SN2LO1 protonsneutrons
FIG. 4: (Colors online) Same as Fig. 2 but for the spin-orbit fields given in Eq. (36).
B. Asymptotic properties
Before entering the numerical details of the solution of Eq. (35), we want to prove that a solution with a well-behaved asymptotic behaviour (origin and infinity) exists. It has been well established for the standard Skyrmesecond-order differential equation [5] that the radial part of the wave-function Eq. (12) behaves as R α ∝ r l +1 at theorigin so that it compensates the behavior of the centrifugal term which diverges as 1 /r . In the case of the presentfourth-order differential equation, this result is a priori no longer true. We thus assume that R α ( r ) ∝ r β around r = 0and determine the possible physical value for β . We insert it in HF equations given in Eq. (35) and we obtain ǫ α r = β ( β − β − β − A + β ( β − β − A r + β ( β − A R r + βA R r + A R r + ℓ ( ℓ + 1) (cid:2) β ( β − A C + βA C r + A C r + ℓ ( ℓ + 1) A CC (cid:3) + (cid:18) j ( j + 1) − l ( l + 1) − (cid:19) (cid:2) W R r + ℓ ( ℓ + 1) W C r + βW R r + β ( β − W R r (cid:3) . (51)All non relevant single-particle quantum numbers are omitted in this discussion to make the notation lighter. Byinspecting the formal expressions of the coefficients A i in Eqs. (37-45), we observe that some fields diverge aroundorigin A C −−−→ r → r , (52) A C −−−→ r → r . (53)0The term A R does not diverge since the derivative of the density is zero at the origin. This is typically the case ofnuclear densities, even in the case of strong shell effects [44]. The spin-orbit fields have no divergence, so we can dropthem. To have a well-behaved wave-function at r = 0 we thus need to check that only the following terms give zero β ( β − β − β − A + β ( β − β − A r + β ( β − A R r + βA R r + A R r + ℓ ( ℓ + 1) (cid:2) β ( β − A C + βA C r + A C r + ℓ ( ℓ + 1) A CC (cid:3) ≈ . (54)First we notice that A , A R and A R do not diverge at the origin. When multiplied by powers of r , they thus go tozero at the origin. By inspecting Eqs. (37-45), we can then notice that to leading order the following relations hold A C = − A A C = 4 A A C = − A A CC = A , (55)so that we can simplify β ( β − β − β − A + ℓ ( ℓ + 1) [ β ( β − A C + βA C + A C + ℓ ( ℓ + 1) A CC ] ≃ . (56)We finally obtain β − β + β (cid:0) − ℓ − ℓ + 11 (cid:1) + 6 β (cid:0) ℓ + ℓ − (cid:1) + ℓ ( ℓ + 1) (cid:0) ℓ + ℓ − (cid:1) ≃ . (57)This equation has 4 solutions β = 2 − ℓ, β = − ℓ, β = ℓ + 1 , β = ℓ + 3 . (58)The first two solutions diverge for some specific values of ℓ and can not represent the physical behaviour of the radialwave function. The last two solutions are physically well-behaved but since the nuclear density needs to be non-zeroat the center of the nucleus, only the solution β = ℓ + 1 can be accepted. The radial part has therefore the samebehaviour for N1LO and N2LO. At infinity, all the fields vanish as one can easily see from Figs.2-4, thus we canrecover the typical asymptotic behaviour of the solutions of the N1LO functional. C. Numerical methods to solve 4th order equations
The solution of HF equations with 4th order derivative terms represent a major numerical challenge. The standardtechnique for N1LO is usually to project the HF equations on an Harmonic Oscillator basis, since one can use particularproperties of orthonormal polynomials to avoid the explicit numerical derivation [29]. However, the main inconvenientis the slow convergence as a function of the number of basis states, as compared to the solution of the HF equationsvia direct integration [30]. We have thus decided to develop a new numerical solver named
WHISKY [45]: the code hasbeen built in a modular way so it can accept the central part of the N ℓ LO Skyrme pseudo-potential with ℓ = 1 , , ℓ, j, q )-block. As a consequence the number of basis functionsgrows quite quickly, especially when we include pairing correlations (see Sec. IV D) so that we introduced an auxiliaryWood-Saxon (WS) basis and an additional energy cut-off. Since the WS wave-functions are reasonably close to thefinal single-particle solutions, the number of basis states to ensure convergence is quite reduced. An alternative tothe WS basis would be the use of the self-consistent HF basis. However, we did not explore this possibility: since weare not currently working with very neutron rich nuclei, a WS approximation is expected to give a result close to thefinal solution. We plan to add this option to explore the properties of the extended N2LO functional close to stabilityin the next version of the code.In Fig. 5, we compare the accuracy of our HF code against the HF code named LENTEUR [48, 49] as a function of theintermediate WS basis size. The calculations are done in both cases using SLy5 interaction [24] with Coulomb includedand a mesh of h = 0 .
05 fm within a box of 20 fm. It is worth reminding that the code
LENTEUR works with a similartwo-basis method: HF and r-space representation with direct integration of HF equation in coordinate space [40]. Thetotal energy difference for different nuclei obtained with the two codes is defined as ∆ E = | E WHISKY − E LENTEUR | .We observe that the accuracy of our code is very good in a reasonably small basis size. By considering states upto 300 MeV we obtain an accuracy of ≈ Pb using
LENTEUR and
WHISKY witha cut-off of 300 MeV in WS basis. We see that the agreement is very good (8 keV at worst). We conclude that the1 -8-6-4-2 0 2 0 200 400 600 800 1000 l og ( ∆ E ) Energy cutoff on the Wood-Saxon basis [MeV] Ca Pb FIG. 5: (Colors online) Precision obtained with WHISKY against LENTEUR as a function of the cutoff energy in the Wood-Saxon basis for Ca (+) and
Pb ( × ). See text for details. Pb[MeV]
WHISKY LENTEUR
Total Energy -1636.10 -1636.10 Kinetic energy 3874.7 Field energy -6209.6 -6209.6 Spin-orbit -99.081 -99.081Direct Coulomb 829.143 829.143Exchange Coulomb -31.314 -31.314TABLE I: Energies obtained by WHISKY and LENTEUR with self-consistent HF calculations using the SLy5 interaction. Thedifferences appear on the last digits and are written in bold. basis size we have chosen is clearly an excellent compromise of efficiency and accuracy since all energy contributionsare described by the two codes at the keV level of accuracy. This cutoff is consequently used in the fit.The code
LENTEUR accepts only N1LO Skyrme-like functionals. Therefore, in order to test the energy contributionof the terms originated from higher order derivatives, we benchmarked our code against the latest version of
MOCCA [50,51].
MOCCA is a 3D solver working in a cubic box and using imaginary-time algorithm to solve the HF equations [52].For the current comparison, we used a mesh of dx = 0 . D. Pairing correlations
Once we move away from closed-shell nuclei, we need to consider extra pairing correlations [53]. To this purpose,we have generalized the
WHISKY code to solve the complete HFB equations. Since we use a two-basis method, we firstsolve the HF equations in coordinate space and then we transform back to the WS basis. The HFB equations in this2
Pb[MeV]
WHISKY MOCCA
Total Energy -1539.2 -1539.2 Total energy N2LO 89.
E[(∆ ρ ) ] 4.39 E[ ρQ ] 37.4 E[ τ ] 27.2 E[ τ µν τ µν − τ µν ∇ µ ∇ ν ρ ] 19.8 E[ K µνκ K µνκ ] 0.0546 E[ J µν V µν ] 0.3385 TABLE II: Comparison of the results for WHISKY and MOCCA: different N2LO functional contributions to the total energyafter a self-consistent calculation with a toy N2LO interaction. The discrepancies are presented in bold. basis read [54] X α ′ ( h lj,qα ′ α − µ qF ) U nlj,qα ′ + X α ′ ∆ lj,qαα ′ V nlj,qα ′ = E nlj,q U nlj,qα , (59) X α ′ ∆ lj,qαα ′ U nlj,qα ′ − X α ′ ( h lj,qα ′ α − µ qF ) V nlj,qα ′ = E nlj,q V nlj,qα , (60)where µ qF is the chemical potential and U nlj,qα and V nlj,qα are the Bogoliubov amplitudes for the quasiparticle of energy E nlj,q , α is the index of the WS basis and n is the index of the quasi-particle state. The field h lj,qα ′ α is derived fromEq. (36) via a unitary transformation.For the pairing channel we used a simple pairing interaction of the form [55, 56] v ( r , r ) = V q (cid:20) − η (cid:18) ρ ( R ) ρ sat (cid:19)(cid:21) δ ( r ) , (61)where R = ( r + r ) / r = r − r is their mutual distance.In the present article we use the so-called volume shape [57] with parameter V n = V p = −
200 MeV.fm , η = 0,and ρ sat = 0 .
16 fm − . Since this interaction has an ultraviolet divergency [58], we use a simple cut-off procedure inquasi-particle space E cut = 60 MeV. For more details on this topic we refer to Ref. [59]. The choice of the pairinginteraction is crucial to determine properties of nuclei far from stability [60, 61]. At present we followed the Saclay-Lyon fitting protocol, so we decoupled the problem in two steps. After the complete fit of the N2LO functional, the V parameters can be fixed to pairing effects. In this article we have used prefixed values of the V parameters, butwe plan to extend our fitting procedure to take into account also pairing effects more precisely [62]. The pairinginteraction for protons and neutrons is not necessary the same, since Coulomb effects should also taken into accountin the calculation of proton Cooper pairs [63]. We plan to include such effects in the next version of the code.In Tab.III, we compare WHISKY against
LENTEUR for
Sn and SLy5 interaction plus volume pairing Eq.61. Weobserve that the accuracy is remarkably high. The small discrepancy of 4 keV originates from a different definition ofcut-off in single particle states:
LENTEUR operates with a cut-off on the total angular momentum j of the quasi-particlestates entering the calculation, while WHISKY operates with a cut-off on the orbital angular momentum.In Fig.6, we compare the total density ρ ( r ) for Sn obtained with the two codes and also the pairing density ˜ ρ .Following Refs. [40, 64] we define it as˜ ρ q ( r ) = − X nlj (2 j + 1)4 π V nlj,q ( r ) U nlj,q ( r ) r , (62)where V nlj,q ( r ) , U nlj,q ( r ) are the quasi-particle amplitudes expressed in r-space. The agreement is excellent, thusdemonstrating the very high accuracy of our new NEDF solver.3 Sn[MeV]
WHISKY LENTEUR
Total energy -1018.81 -1018.81 Kinetic energy 2188.1 Field energy -3485.1 -3485.1 Spin-orbit energy -55.00 -55.00 Coulomb (direct) 367.336 367.336Coulomb (exchange) -19.147 -19.147Neutron pairing energy -15.01 -15.01 TABLE III: Comparaison between the energies obtained by WHISKY and lenteur with self-consistent HFB calculations usingthe SLy5 interaction. See text for details [f m - ] r [fm] SN2LO1 ρ ρ ~SLy5 ρ ρ ~ FIG. 6: (Colors online) Isoscalar particle density and pairing density for
Sn obtained with a self-consistent mean-fieldcalculation with the SLy5 interaction.
V. FIT OF N2LO INTERACTION
To fit the N2LO pseudo-potential we adopted a modified version of the Saclay-Lyon fitting protocol [24, 65]:the protocol includes here both properties of some selected double-magic nuclei and some basic properties of theinfinite nuclear medium as saturation density, incompressibility and the equation of state of pure neutron matter(PNM) derived from realistic nucleon-nucleon interactions [66]. We consider all terms of the interaction, and we treatspurious center of mass motion with the usual one-body approximation [1, 24]. We also assume equal neutron andproton masses and we use the value ~ m = 20 . [24]. A. Fitting protocol
To obtain the parameters of the pseudo-potential we need to minimize the following penalty function [14] χ = M X i =1 ( O i − f i ( p )) ∆ O i , (63)where the sum runs over all the M (pseudo)-observables O i we want to constraint in our fit, f i is the value obtainedwith our solver for a given array of parameters p = { t , t , t , . . . } , while ∆ O i is the weight we give to each pointin the fit. Let’s mention that ∆ O i does not correspond necessarily to the experimental uncertainty. In Tab. IV,we give the actual constraints we used to build the χ function in Eq. (63). On top of this constraints, we paidparticular attention in tuning the spin-orbit parameter W to some specific range of acceptable values. Finally, itis worth noticing that during the χ minimisation the parameters p cannot vary freely: in order to avoid finite-size4 Fit Constraints O i ∆ O i Units Reference
Infinite nuclear matter ρ sat − [68, 69]E/A ( ρ sat ) -16.0000 0.2 MeV [68, 69] m ∗ /m K ∞ J EoS PNM [66]E/N ( ρ =0.1) 11.88 2.0 MeVE/N ( ρ =0.3) 35.94 7.0 MeVE/N ( ρ =0.35) 44.14 9.0 MeV Stability [10]INM(S,M,T) ρ crit ≥ .
24 asymmetric fm − constraint Finite nuclei
Binding energies [72] Ca -342.02300 1.5 MeV Ca -415.98300 1.0 MeV Ni -483.95300 1.5 MeV
Sn -825.13000 1.5 MeV
Sn -1102.67300 1.0 MeV
Pb -1635.86100 1.0 MeV
Proton radii [73] Ca 3.38282 0.03 fm Ca 3.39070 0.02 fm Ni 3.66189 0.03 fm
Sn 4.64745 0.02 fm
Pb 5.45007 0.02 fm
Parameter W TABLE IV: Constraints O i used in the fitting procedure and the associated error ∆ O i . See text for details. instabilities [10], the critical densities in all channels are computed at each iteration, and an asymmetric constraint isimposed in terms of a penalty function χ fs = X α exp − β ( O α − ρ crit ) , (64)where O α =( S,M,T ) is the lowest density at which an instability appears in symmetric nuclear matter (SNM). ρ crit isan empirical value defined in Refs. [10, 11] to avoid unphysical instabilities. β is an arbitrary parameter ( β = 10 here)fixed in such a way that the penalty function grows very fast when we approach the critical density from below, butgives no contribution when above it. This constraint is applied in all channels for which we calculate the responsefunction of the system (see Sec. V B). Finite-size instabilities may also have important impact at high density onastrophysical applications such as the neutrino mean free path [67]. However, in this work, we concentrate ourselveson finite-size instabilities only in densities ranges that are relevant for finite nuclei. In other words, we allow in thispreliminary work the appearance of instabilities at densities above ρ crit which is slightly above saturation density.At the end of the minimisation procedure, we have obtained the parameters p = { t , t , t , . . . } given in Tab. V.Notice that the exponent α of the density dependent term has been fixed from the beginning (see Sec. V C). Fromthe table, it is difficult to judge the quantitative relative importance of the different parameters. A way to bypass theproblem is to use the concept of naturalness . Following Ref. [74] we multiply each N2LO coupling constant by S = f l − π Λ n + l − , (65)where f π = 93 MeV is the pion decay constant, Λ = 687 MeV, l is the power of the density of the correspondingterm and n is the order. Special treatment is required for the density dependent coupling constant. See Ref. [74]5 SN2LO1 n i t ( n ) i [MeVfm n ] x ( n ) i t = 13707 . α ) ] x = 0 . α = 1/6 W = 117.904418 [MeVfm ]TABLE V: Numerical values for N2LO parameters.SN2LO1Natural units C ρ -1.06 C ρ C ρ [ ρ α ] 13.0 C ρ [ ρ α ] -12.1 C τ C τ C ∆ ρ -1.06 C ∆ ρ C ∇ J -1.22 C ∇ J -0.406 C T -0.0882 C T -0.816 C (∆ ρ ) -0.115 C (∆ ρ ) C Mρ -0.288 C Mρ C Ms C Ms -0.0162TABLE VI: Values of the parameters of the N2LO pseudo-potential expressed in natural units. for details. It is important to keep in mind that the value of Λ is somehow arbitrary since it has been derived inRef. [74] by observing the behaviour of several N1LO functionals. The results are presented in Tab. VI. Owing tothe arbitrariness of the value of Λ, one should not look too close to the actual numbers, but only to the order ofmagnitude. By inspecting the table, we clearly observe that there is a natural hierarchy in the coupling constants: theN2LO coupling constants are one order of magnitude smaller than the N1LO ones. This is a very important aspectsince the entire idea behind the N ℓ LO expansion is to have a fast convergence: from these results, we can expect thatwithin this scheme, the N3LO coupling constants would be another order of magnitude smaller.
B. Finite-size instabilities
As discussed in the introduction, several effective interactions are biased by spurious instabilities [12, 13, 75]. Toavoid such a problem, we have developed in Ref. [16] a new fitting protocol based on the LR formalism [76]. Fromprevious analysis of Refs [10, 11], we have noticed that when a pole in the response function appears at densities lowerthan ≈ . F parameters, all the Landau inequalities [75] are respected up to two times the saturation density. The onlyinstability appears in the G ′ parameter at ρ ≈ .
35 fm − . This does not represent a major issue for this study sincewe do consider only finite-nuclei and and not astrophysical applications [83]. In Fig. 8, we show the position of the -2.0-1.00.01.02.03.0 F l y G l l=0l=1l=2 -1.0-0.50.00.51.01.52.0 0 0.2 0.4 0.6 0.8 F l’ y ρ [fm -3 ] 0 0.2 0.4 0.6 0.8 G l’ ρ [fm -3 ] FIG. 7: (Colors online) Landau parameters in SNM for the SN2LO1 pseudo-potential as a function of the density of the system.See text for details. critical densities obtained in SNM as a function of the transferred momentum q . The LR is calculated for each spin(S) spin projection (M) and isospin (I) channel (S,M,I). See Ref [12] for more details on the adopted notation. Weobserve no finite-size instabilities, apart from the physical spinodal one [77], around saturation density. This meansthat our interaction is well stable in all spin-isospin channels [10, 11]. This results confirm our preliminary findingsin Ref. [16]: the LR formalism can be considered as a very simple tool to be added in a fitting procedure to avoidexploring regions of parameters that induce unphysical instabilities. ρ ρ crit ρ [f m - ] q [fm -1 ] Canal S=0 M=0 I=0Canal S=0 M=0 I=1Canal S=1 M=0 I=0Canal S=1 M=0 I=1Canal S=1 M=1 I=0Canal S=1 M=1 I=1
SN2LO1
FIG. 8: (Colors online) Critical densities in SNM as a function of transferred momentum q . The horizontal dashed linesrepresent the saturation density ρ and the critical density ρ crit . -20-10010203040 Symmetric matter E / A [ M e V ] Neutron matter E / N [ M e V ] SLy5SN2LO1BHF
Spin-polarized matter E / A [ M e V ] ρ [fm -3 ] 0 0.1 0.2 0.3 0.4 0 20 40 60 80 100 120 Neutron spin-polarized matter E / N [ M e V ] ρ [fm -3 ] FIG. 9: (Colors online) Equation of state for SNM and PNM obtained with the N2LO Skyrme interaction. The squaresrepresent the values obtained from BHF calculations.
C. Infinite nuclear matter
In our fitting protocol, we include information of the infinite nuclear medium. Following Ref. [24], we have used asa constrain three points of the EoS in PNM dervied in Ref [66]. We can now benchmark our results against other wellknown EoS as the one derived via Brueckner-Hartree-Fock (BHF) [84]. In Fig. 9, we compare the EoS for symmetricmatter and neutron matter obtained with BHF and the SN2LO1 interaction. For completeness the results with SLy5are also given. The SN2LO1 follows quite closely the BHF results, and in particular the EoS of PNM up to 3 timessaturation density. Beyond this point the EoS becomes slightly softer. We remind the reader that SLy5 and SN2LO1follow each other quite closely in PNM at low-density since they have been constrained on the same points in thisdensity region.On the same figure, we also give the results for spin-polarised symmetric matter and spin-polarised pure neutronmatter and compare SLy5 and SN2LO1 results. Although these two quantities have not been fitted explicitly, weobserve a qualitative similar behaviour in the two functionals. For completeness, in Tab. VII, we give the main featuresof the EoS of SN2LO1, i.e. saturation density ρ , incompressibility K ∞ , symmetry energy J and slope of symmetryenergy L (not fitted). The values we obtained are in agreement with the existing constraints [85].As already discussed in Ref. [24], there is a strong model correlation for N1LO between the nuclear incompressibilityand the effective mass. In our case, the correlation between K ∞ and m/m ∗ is of course different since the newparameters give us more freedom in adjusting these two values. It can be calculated analytically in infinite matterwith the result K ∞ = − α + 1) EA ( ρ ) + 35 ~ m k F (cid:16) α − − α − mm ∗ (cid:17) + 3140 C (4)0 ρk F (3 α + 10) . (66)In Fig. 10, we observe that to obtain a reasonable value of the nuclear incompressibility, the allowed range for α is α ∈ [1 / , / D. Finite nuclei
In this section, we analyse the properties of finite nuclei obtained with the extended Skyrme pseudo-potential. InFig. 11, we show the energy difference ∆ E between the experimental values and the ones calculated using either8 SN2LO1 SLy5 ρ [fm − ] 0.162 0.1603E/A( ρ ) [MeV] -15.948 -15.98 K ∞ [MeV] 221.9 229.92 J [MeV] 31.95 32.03 L [MeV] 48.9 48.15 m ∗ /m m * / m K ∞ [MeV] α = 1/6 α = 1/3 α = 1/2 α = 1SN2LO1 FIG. 10: (Colors online) Correlation between the effective mass and the nuclear incompressibility in infinite nuclear matter fordifferent values of the power of the density dependent term.
SLy5 or SN2LO1 for the few selected double-magic nuclei used in the fit. The results obtained with SN2LO1 areof the same quality as SLy5. Moreover they are all very close to the tolerance ∆ O i we used for the fit given inTab. IV. In Fig. 12, we compare the differences of proton radii ∆ r p obtained with SLy5 and our new pseudo-potentialSN2LO1. In this case we see that SN2LO1 behaves marginally better than SLy5 giving a result typically closer tothe experimental values. It is worth noticing that compared to SLy5, we have few additional constraints concerningfinite-size instabilities that were not present in the original fitting protocol of SLy5. The closest functional to SN2LO1,in terms of fitting protocol, is represented by SLy5 ∗ [16]. We do not report here the direct comparison, but we havechecked that the results are qualitatively the same.In Figs.13, we compare the differences between the binding energies calculated for isotopic (isotonic) chains with ∆E = E th - E exp [MeV] Ca Ca Ni Sn Sn Pb -5 -4 -3 -2 -1 0 1 2 3 4 5 • SN2LO1 • SLy5
FIG. 11: Difference of binding energies obtained with SN2LO1 and SLy5 and experimental values extracted from Ref. [72]. ∆r p = r th - r exp [10 − fm] Ca Ca Ni Sn Pb • •••• • •••• • SN2LO1 • SLy5 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7FIG. 12: Proton radii difference of two interactions (SN2LO1/SLy5) calculated with WHISKY with experimental radii obtainedin Ref. [73].
Z(N)=20, 28, 50, 82 for our extended Skyrme interaction. The experimental measurements are taken from Ref. [72].On the same figure, we also report the values obtained with SLy5. Notice that we did not optimised the value ofthe pairing strength to improve the reproduction of experimental data. Moreover, since the effective masses arenumerically quite similar for SLy5 and SN2LO1, we used exactly the same pairing interaction. The main feature we -6-5-4-3-2-101 10 15 20 25 30 35 40 ∆ E t o t Neutron number Ca
50 55 60 65 70 75 80 85 90-4-3-2-1 0 1 2 3 4 ∆ E t o t Neutron numberSN2LO1SLy5 Sn -7-6-5-4-3-2-1012 20 25 30 35 40 45 50 ∆ E t o t Neutron number Ni
95 100 105 110 115 120 125 130 135-2 0 2 4 6 8 ∆ E t o t Neutron number Pb -1.5-1.0-0.50.00.51.01.5 20 22 24 26 28 30 32 ∆ E t o t Proton number
N=28
25 30 35 40 45 50-5-4-3-2-1 0 1 2 ∆ E t o t Proton numberSN2LO1SLy5
N=50 -3.0-2.0-1.00.01.02.03.04.0 45 50 55 60 65 70 ∆ E t o t Proton number
N=82
80 82 84 86 88 90 92 0 5 10 15 ∆ E t o t Proton number
N=126
FIG. 13: (Colors online) Systematic comparison of binding energies, expressed in MeV, for isotopic (isotonic) chains calculatedwith our extended Skyrme interaction SN2LO1 and experimental ones. On the same figure we also compare with the SLy5parametrisation. See text for details. observe is the strong arch-like structures. This is the main drawback of a fitting protocol that fixes a very limitednumber of nuclei. A better fitting protocol has been designed for example for UNEDF functionals [18–20] and weplan to use it for a systematic exploration of the parameter space of higher order terms. In Fig. 14, we compare theproton radii. The data are taken from Ref. [73]. The new interaction is fairly closer to experimental data than theoriginal SLy5 and the main trends are reproduced. One of the biggest discrepancy we observe in the data is relatedto the anomalous isotopic dependence of proton radii of calcium isotopes. With the current parametrisation we havenot been able to reproduce both Ca and Ca. A recent article [88] suggests that a different form of the pairingfunctionals based on Fayans form [89] may be the key to solve this anomaly, while the specific form of the functionalused for the calculation of the central potential is not relevant. Since we did not fix any particular pairing functionalin our fit, we plan to test the results of Ref. [88] with our new functional.Finally, we have explored the behaviour of single particle spectra. In Fig. 15, we compare the Hartree-Fock neutronsingle particle states for Ca obtained using SLy5 and SN2LO1. The values are compared with the experimentalvalues extracted from Ref. [90]. The HF states obtained with the two functionals are very close to each other. SN2LO1shows a slight compression of the spectrum, but this is simply related to a slightly larger effective mass (see Tab. VII).Similar behaviour is also observed in Fig. 16 for neutron single-particle states in
Pb.As discussed in Sec. IV, the higher order gradient terms induce three extra spin-orbit fields Eqs.47-50. In principle0this should provide us with some extra flexibility compared to a standard Skyrme interaction. However, the majorproblem encountered in this first analysis is to find the right observable that may let us explore a new region ofparameter space that may increase their importance. We recall that we neglected completely tensor terms at N2LOlevel, which means two extra tensor parameters [25, 27]. This could also give an extra freedom to correct some knownanomaly in the shell evolution of some particular states [91]. The exploration of this particular aspect is currentlyunder investigation. r p Neutron number Ca
50 55 60 65 70 75 80 85 90 4.4 4.45 4.5 4.55 4.6 4.65 4.7 r p Neutron number
SN2LO1SLy5Exp Sn r p Neutron number Ni
95 100 105 110 115 120 125 130 135 5.25 5.3 5.35 5.4 5.45 5.5 5.55 r p Neutron number Pb FIG. 14: (Colors online) Systematic comparison of proton radii. Experimental data are taken from Ref. [73].
VI. CONCLUSIONS
In the present article, we have discussed the formalism to include fourth-order gradient terms of the N2LO Skyrmeinteraction. We have derived the functional, the complete expression of the densities in the case of spherical symmetryand the corresponding HF equation. The resulting 4-th order differential equation has been solved with a newnumerical code named
WHISKY . This code has been tested against two different HFB solvers to check numerical [ M e V ] -24-20-16-12-8-40 Exp SLy5 SN2LO1 d / s / d / f / p / p / f / FIG. 15: Neutron single-particle energies around the Fermi energy in the Ca for SLy5 and SN2LO1 parametrisations. Theexperimental values are taken from Ref. [90]. See text for details. [ M e V ] -12-10-8-6-4-20 Exp SLy5 SN2LO1 h / f / i / p / f / p / g / i / FIG. 16: Same as Fig. 15, but for
Pb. accuracy of the new solver. Thanks to this new code, we have been able to perform for the very first time a completefit of a stable N2LO Skyrme interaction including finite-nuclei. This achievement has been made possible by the useof the Linear Response formalism as a tool to prevent unphysical instabilities.For the very first time, we thus have been able to prove that it is possible to go beyond the standard Skyrmeinteraction by including physically motivated terms. Thanks to the work on the foundations of various non-relativisticeffective interactions [25], we have been able to clarify the inner nature of the higher order gradient terms in theextended N ℓ LO Skyrme pseudo-potential. The LR formalism we have been able to solve also the long-standingproblem of finite-size instabilities in effective functionals. Finite-size instabilities seem to appear in various functionalsnot only the Skyrme-like ones [13]. The LR formalism thus represents a simple tool that should be included in allmodern fitting protocol to avoid the appearance of non physical results.Combining all the previous results, we have been able to derive the complete set of parameters of the N2LOpseudo-potential, named SN2LO1 in this paper. We have compared its performances on both infinite nuclear matter(pseudo)-observables as well as ground state properties of some selected nuclei. The global performances are of thesame quality as the standard SLy5. However, it is very important to underline here that since SN2LO1 has fouradditional parameters compared to SLy5, we have imposed extra stability constraints to our functional: SLy5 has afinite-size instability in the spin-channel and thus can not be used to perform calculations where the time-odd channelis open. To the best of our knowledge, SN2LO1 is free from pathologies and it can be safely used in various numericalcodes.Finally we insist on the fact that the higher order terms introduce several new features as for example three newspin-orbit fields that have not been completely investigated in this article and may give rise to new properties of thefunctional: N2LO clearly offers some new degrees of freedom and goes beyond N1LO.
Acknowledgments
We are grateful to W. Ryssens for providing us with the
MOCCA results, as well to K. Bennaceur for providing uswith the
LENTEUR code as well for fruitful discussion. We also acknowledge interesting discussions with M. Bender.The work of J.N. has been supported by grant FIS2014-51948-C2-1-P, Mineco (Spain).2
Appendix A: Coupling constants
In this section we give the explicit expressions of the new coupling constants of N2LO functional in terms of Skyrmeparameters. The expression of the coupling constants for the standard Skyrme functional can be found in Ref. [36]. C (∆ ρ ) = h t (4)1 − t (4)2 (cid:16) x (4)2 (cid:17)i (A1) C (∆ ρ ) = − h t (4)1 (cid:16) x (4)1 (cid:17) + t (4)2 (cid:16) x (4)2 (cid:17)i (A2) C Mρ = h t (4)1 + t (4)2 (cid:16) x (4)2 (cid:17)i (A3) C Mρ = − h t (4)1 (cid:16) x (4)1 (cid:17) − t (4)2 (cid:16) x (4)2 (cid:17)i (A4) C Ms = − h t (4)1 (cid:16) − x (4)1 (cid:17) − t (4)2 (cid:16) x (4)2 (cid:17)i (A5) C Ms = − h t (4)1 − t (4)2 i (A6) Appendix B: Densities in Cartesian representation
We define the density matrix in coordinate space as in [41] ρ q ( r σ, r ′ σ ′ ) = 12 ρ q ( r , r ′ ) δ σσ ′ + 12 s q ( r , r ′ ) h σ ′ | ˆ σ | σ i , (B1)where ρ q ( r , r ′ ) = X σ ρ q ( r σ, r ′ σ ′ ) (B2) s q ( r , r ′ ) = X σσ ′ ρ q ( r σ, r ′ σ ′ ) h σ ′ | ˆ σ | σ i . (B3)The Skyrme energy density functional up to 2nd order is composed by seven local densities whose explict expressioncan be found, for example, in Ref. [36]. The extension to fourth order requires the definition of six additional localdensities τ µν,q ( r ) = ∇ µ ∇ ′ ν ρ q ( r , r ′ ) | r = r ′ (B4) K µνκ,q ( r ) = ∇ µ ∇ ′ ν s κq ( r , r ′ ) | r = r ′ (B5)Π µ,q ( r ) = ∇ · ∇ ′ j µ,q ( r , r ′ ) | r = r ′ (B6) V µν,q ( r ) = ∇ · ∇ ′ J µν,q ( r , r ′ ) | r = r ′ (B7) Q q ( r ) = ∆∆ ′ ρ q ( r , r ′ ) | r = r ′ (B8) S µ,q ( r ) = ∆∆ ′ s µ,q ( r , r ′ ) | r = r ′ (B9)Similarly to the spin-current pseudo-tensor J µν,q ( r ), the density τ µν,q ( r ) can be decomposed into a pseudo-scalar,vector and traceless pseudo-tensor term. For more details we refer to Ref. [76]. [1] M. Bender, P.-H. Heenen, and P.-G. Reinhard, Rev. Mod. Phys. , 121 (2003).[2] E. Perli´nska, S. G. Rohozi´nski, J. Dobaczewski, and W. Nazarewicz, Phys. Rev. C , 014316 (2004).[3] F. Raimondi, B. G. Carlsson, and J. Dobaczewski, Phys. Rev. C , 054311 (2011).[4] T. H. R. Skyrme, Nucl. Phys. , 615 (1959).[5] D. Vautherin and D. M. Brink, Phys. Rev. C , 626 (1972).[6] S. Goriely, N. Chamel, and J. M. Pearson, Phys. Rev. Lett. , 152503 (2009).[7] T. Lesinski, K. Bennaceur, T. Duguet, and J. Meyer, Phys. Rev. C , 044315 (2006).[8] N. Schunck, J. Dobaczewski, J. McDonnell, J. Mor´e, W. Nazarewicz, J. Sarich, and M. Stoitsov, Phys. Rev. C , 024316(2010).[9] S. Fracasso, E. B. Suckling, and P. Stevenson, Physical Review C , 044303 (2012). [10] V. Hellemans, A. Pastore, T. Duguet, K. Bennaceur, D. Davesne, J. Meyer, M. Bender, and P.-H. Heenen, Phys. Rev. C , 064323 (2013).[11] A. Pastore, D. Tarpanov, D. Davesne, and J. Navarro, Phys. Rev. C , 024305 (2015).[12] A. Pastore, D. Davesne, and J. Navarro, Phys. Reports , 1 (2015).[13] A. De Pace and M. Martini, Phys. Rev. C , 024342 (2016).[14] J. Dobaczewski, W. Nazarewicz, and P. Reinhard, J. Phys. G: Nucl. Part. Phys. , 074001 (2014).[15] T. Nikˇsi´c and D. Vretenar, Phys. Rev. C , 024333 (2016).[16] A. Pastore, D. Davesne, K. Bennaceur, J. Meyer, and V. Hellemans, Phys. Scripta T154 , 014014 (2013).[17] G. Bertsch, D. Dean, and W. Nazarewicz, SciDAC Review , 42 (2007).[18] M. Kortelainen, T. Lesinski, J. Mor´e, W. Nazarewicz, J. Sarich, N. Schunck, M. V. Stoitsov, and S. Wild, Phys. Rev. C , 024313 (2010).[19] M. Kortelainen, J. McDonnell, W. Nazarewicz, P.-G. Reinhard, J. Sarich, N. Schunck, M. Stoitsov, and S. Wild, Phys.Rev. C , 024304 (2012).[20] M. Kortelainen, J. McDonnell, W. Nazarewicz, E. Olsen, P.-G. Reinhard, J. Sarich, N. Schunck, S. M. Wild, D. Davesne,J. Erler, et al., Phys. Rev. C , 054314 (2014).[21] T. Duguet, M. Bender, J.-P. Ebran, T. Lesinski, and V. Som`a, Eur. Phys. J. A , 1 (2015).[22] B. G. Carlsson, J. Dobaczewski, and M. Kortelainen, Phys. Rev. C , 044326 (2008).[23] D. Davesne, A. Pastore, and J. Navarro, J. Phys. G: Nucl. Part. Phys. , 095104 (2013).[24] E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer, Nucl. Phys. A , 710 (1997).[25] D. Davesne, P. Becker, A. Pastore, and J. Navarro, Ann. Phys. (NY) , 288 (2016).[26] B. G. Carlsson and J. Dobaczewski, Phys. Rev. Lett. , 122501 (2010).[27] D. Davesne, J. Navarro, P. Becker, R. Jodon, J. Meyer, and A. Pastore, Phys. Rev. C , 064303 (2015).[28] D. Davesne, A. Pastore, and J. Navarro, Astron. Astrophys. , A83 (2015).[29] B. Carlsson, J. Dobaczewski, J. Toivanen, and P. Vesel`y, Comp. Phys. Comm. , 1641 (2010).[30] N. Schunck, J. D. McDonnell, J. Sarich, S. M. Wild, and D. Higdon, J. Phys. G: Nucl. Part. Phys. , 034024 (2015).[31] D. Davesne, A. Pastore, and J. Navarro, J. Phys. G: Nucl. Part. Phys. , 065104 (2014).[32] J. Dobaczewski and J. Dudek, Phys. Rev. C , 1827 (1995).[33] F. Raimondi, B. G. Carlsson, J. Dobaczewski, and J. Toivanen, Phys. Rev. C , 064303 (2011).[34] J. Sadoudi, T. Duguet, J. Meyer, and M. Bender, Phys. Rev. C , 064326 (2013).[35] H. Sagawa and G. Col`o, Progr. Part. Nucl. Phys. , 76 (2014).[36] T. Lesinski, M. Bender, K. Bennaceur, T. Duguet, and J. Meyer, Phys. Rev. C , 014312 (2007).[37] J. Skalski, Phys. Rev. C , 024312 (2001).[38] D. Varshalovich, A. N. Moskalev, and V. K. Khersonskii, Quantum theory of angular momentum (World Scientific, Singa-pore, 1988).[39] K. Sato, J. Dobaczewski, T. Nakatsukasa, and W. Satu la, Physical Review C , 061301 (2013).[40] K. Bennaceur and J. Dobaczewski, Comp. Phys. Comm. , 96 (2005).[41] P.Ring and P.Schuck, The Nuclear Many-Body Problem (Springer-Verlag Berlin Heidelberg, 1980).[42] J. Banerjee, J. Sound Vib. , 97 (2001).[43] P.-G. Reinhard, D. J. Dean, W. Nazarewicz, J. Dobaczewski, J. A. Maruhn, and M. R. Strayer, Phys. Rev. C , 014316(1999), URL https://link.aps.org/doi/10.1103/PhysRevC.60.014316 .[44] J.-F. Berger, L. Bitaud, J. Decharg´e, M. Girod, and K. Dietrich, Nucl. Phys. A , 1 (2001).[45] P. Becker, Ph.D. thesis, Lyon (2017).[46] N. Schunck, J. Dobaczewski, J. McDonnell, W. Satu la, J. Sheikh, A. Staszczak, M. Stoitsov, and P. Toivanen, Comp. Phys.Com. , 166 (2012).[47] R. Hooverman, Nucl. Phys. A , 155 (1972).[48] V. Rotival and T. Duguet, Phys. Rev. C , 054308 (2009).[49] V. Rotival, K. Bennaceur, and T. Duguet, Phys. Rev. C , 054309 (2009).[50] W. Ryssens, V. Hellemans, M. Bender, and P.-H. Heenen, Comp. Phys. Comm. , 175 (2015).[51] W. Ryssens, Ph.D. thesis, Universit´e libre de Bruxelles (2016).[52] P. Bonche, H. Flocard, and P.-H. Heenen, Comp. Phys. Comm. , 49 (2005).[53] D. M. Brink and R. A. Broglia, Nuclear superfluidity: pairing in finite systems , vol. 24 (Cambridge University Press, 2005).[54] A. Pastore, F. Barranco, R. Broglia, and E. Vigezzi, Phys. Rev. C , 024315 (2008).[55] G. Bertsch and H. Esbensen, Ann. Phys. (NY) , 327 (1991).[56] E. Garrido, P. Sarriguren, E. Moya de Guerra, and P. Schuck, Phys. Rev. C , 064312 (1999).[57] N. Sandulescu, P. Schuck, and X. Vi˜nas, Phys. Rev. C , 054303 (2005).[58] A. Bulgac and Y. Yu, Phys. Rev. Lett. , 042504 (2002).[59] P. Borycki, J. Dobaczewski, W. Nazarewicz, and M. Stoitsov, Phys. Rev. C , 044319 (2006).[60] J. Dobaczewski, I. Hamamoto, W. Nazarewicz, and J. Sheikh, Phys. Rev. Lett. , 981 (1994).[61] A. Pastore, J. Margueron, P. Schuck, and X. Vi˜nas, Phys. Rev. C , 034314 (2013).[62] J. G´omez, C. Prieto, and J. Navarro, Nucl. Phys. A , 125 (1992).[63] H. Nakada and M. Yamagami, Physical Review C , 031302 (2011).[64] J. Dobaczewski, H. Flocard, and J. Treiner, Nuclear Physics A , 103 (1984).[65] K. Washiyama, K. Bennaceur, B. Avez, M. Bender, P.-H. Heenen, and V. Hellemans, Phys. Rev. C , 054309 (2012).[66] R. B. Wiringa, V. Fiks, and A. Fabrocini, Phys. Rev. C , 1010 (1988). [67] A. Pastore, M. Martini, D. Davesne, J. Navarro, S. Goriely, and N. Chamel, Phys. Rev. C , 025804 (2014).[68] C. S. Wang, K. C. Chung, and A. J. Santiago, Phys. Rev. C , 034310 (1999).[69] R. Nayak, V. S. Uma Maheswari, and L. Satpathy, Phys. Rev. C , 711 (1995).[70] J. Blaizot, Phys. Reports , 171 (1980).[71] P.-G. Reinhard, Nucl. Phys. A , 305 (1999).[72] M. Wang, G. Audi, A. Wapstra, F. Kondev, M. MacCormick, X. Xu, and B. Pfeiffer, Chin. Phys. C , 1603 (2012).[73] I. Angeli, At. Data Nucl. Data Tables , 185 (2004).[74] M. Kortelainen, R. J. Furnstahl, W. Nazarewicz, and M. V. Stoitsov, Phys. Rev. C , 011304 (2010).[75] J. Margueron, J. Navarro, and N. Van Giai, Phys. Rev. C , 014303 (2002).[76] P. Becker, D. Davesne, J. Meyer, A. Pastore, and J. Navarro, J. Phys. G: Nucl. Part. Phys , 034001 (2014).[77] C. Ducoin, P. Chomaz, and F. Gulminelli, Nucl. Phys. A , 403 (2007).[78] D. Davesne, A. Pastore, and J. Navarro, Phys. Rev. C , 044302 (2014).[79] L. Landau, Sov. Phys. JETP , 70 (1959).[80] D. Davesne, J. Holt, A. Pastore, and J. Navarro, Phys. Rev. C , 014323 (2016).[81] S.-O. B¨ackman, A. Jackson, and J. Speth, Phys. Lett. B , 209 (1975).[82] W. Zuo, C. Shen, and U. Lombardo, Phys. Rev. C , 037301 (2003).[83] S. Goriely, N. Chamel, and J. Pearson, Phys. Rev. C , 035804 (2010).[84] M. Baldo, I. Bombaci, and G. Burgio, Astron. Astrophys. , 274 (1997).[85] M. Dutra, O. Louren¸co, J. S. S´a Martins, A. Delfino, J. R. Stone, and P. D. Stevenson, Phys. Rev. C , 035201 (2012).[86] D. Lacroix, T. Duguet, and M. Bender, Phys. Rev. C , 044318 (2009).[87] M. Bender, T. Duguet, and D. Lacroix, Phys. Rev. C , 044319 (2009).[88] P.-G. Reinhard and W. Nazarewicz, arXiv: 1704.07430 (2017).[89] S. Fayans, S. Tolokonnikov, E. Trykov, and D. Zawischa, Nucl. Phys. A , 49 (2000).[90] N. Schwierz, I. Wiedenhover, and A. Volya, arXiv: 0709.3525 (2007).[91] G. Col`o, H. Sagawa, S. Fracasso, and P. Bortignon, Phys. Lett. B646