Local nuclear energy density functional at next-to-next-to-next-to-leading order
aa r X i v : . [ nu c l - t h ] O c t Local nuclear energy density functional at next-to-next-to-next-to-leading order
B.G. Carlsson, J. Dobaczewski,
1, 2 and M. Kortelainen Department of Physics, P.O. Box 35 (YFL), FI-40014 University of Jyv¨askyl¨a, Finland Institute of Theoretical Physics, University of Warsaw, ul. Ho˙za 69, 00-681 Warsaw, Poland. (Dated: September 21, 2018)We construct nuclear energy density functionals in terms of derivatives of densities up to sixth,next-to-next-to-next-to-leading order (N LO). A phenomenological functional built in this way con-forms to the ideas of the density matrix expansion and is rooted in the expansions characteristic toeffective theories. It builds on the standard functionals related to the contact and Skyrme forces,which constitute the zero-order (LO) and second-order (NLO) expansions, respectively. At N LO,the full functional with density-independent coupling constants, and with the isospin degree of free-dom taken into account, contains 376 terms, while the functionals restricted by the Galilean andgauge symmetries contain 100 and 42 terms, respectively. For functionals additionally restricted bythe spherical, space-inversion, and time-reversal symmetries, the corresponding numbers of termsare equal to 100, 60, and 22, respectively.
PACS numbers: 21.60.Jz, 21.30.Fe, 71.15.Mb
I. INTRODUCTION
The Skyrme force was introduced into nuclear physicsmore than half a century ago [1, 2] and it is still a con-cept that is widely used in methods striving to determineproperties of nuclei irrespective of their mass number andisospin. However, at present we understand this conceptin a significantly different way than it was originally pro-posed. Indeed, instead of using this force as an effectiveinteraction within the Hartree-Fock (HF) approximation,we rather focus on the underlying Skyrme energy densityfunctional (EDF), without direct references to the effec-tive interaction or HF approximation.In electronic systems, the use of functionals of den-sity is motivated by formal results originating fromthe Hohenberg-Kohn [3] and Kohn-Sham theorems [4],whereby exact ground-state energies of many-fermionsystems can be obtained by minimizing a certain exactfunctional of the one-body density. This led to numer-ous extensions and applications, now collectively knownunder the name of density functional theory (DFT) [5–7].The fact that properties of electronic systems are gov-erned by the well-known Coulomb interaction allows forderivations of functionals from first principles, by whichtoken this approach can proudly be called a theory. Fornuclear systems, the luxury of knowing the exact interac-tion is not there, so the analogous approaches developedin this domain of physics carry the name of the EDFmethods.In this article we construct a phenomenological nuclearEDF based on strategies that are proper to effective the-ories [8]. There, guiding principles [9] are based on: (i)appropriate choice of effective fields, (ii) building effec-tive Lagrangian or Hamiltonian densities restricted onlyby symmetry principles, (iii) employing ideas of powercounting. In the low-energy nuclear structure, correctfields can probably be associated with nonlocal one-bodynuclear densities. Then, functionals of densities acquirethe meaning of effective Hamiltonian densities. Although a formal construction of power-counting schemes is notyet available, ideas based on the density matrix expan-sion (DME) [10–16] (see also a recent example of an ap-plication to electronic systems in Ref. [17]) can be usedto propose expansions in terms of moments of effectivenuclear interactions, or equivalently, in orders of deriva-tives acting on the one-body densities. This is preciselythe strategy we are going to follow in the present study.Effective field theories (EFT) were recently extensivelyapplied in analyzing properties of nuclear systems. Herewe are not able to give an even shortest possible review ofthis rapidly developing area of physics, but let us mentiontwo specific examples.Firstly, the nucleon-nucleon (NN) scattering propertieswere very successfully described by employing the effec-tive nucleon-pion Lagrangian at next-to-next-to-next-to-leading order (N LO), see Ref. [18] and references citedtherein. This showed that the EFT expansion is capa-ble of grasping the main features of nuclear interactionsat low energies, without explicitly invoking microscopicfoundations in terms of, e.g., heavy meson exchanges.Secondly, methods using the harmonic-oscillator effec-tive operators have been developed up to N LO, to beemployed within the shell-model approaches [19]. There,the N LO expansion was explicitly expressed in the formof pseudopotentials that contain derivatives up to sixthorder. Evidently, such pseudopotentials are exact equiv-alents of higher-order Skyrme-like forces. When averagedwithin the HF approximation, they would lead to EDFsdepending on derivatives of densities up to sixth order.This allows us to label our approach with the traditionalname of the N LO expansion too.There is also a recent significant effort in deriving nu-clear EDFs directly from low-energy QCD within chi-ral perturbation theory, see, e.g., Refs. [20–22]. Thismay have a potential of providing new important insightinto the precise structure of terms in the EDF, whileat present we are bound to proceed phenomenologically,with only the symmetry constraints available, as is donein the present study.In all rigorous EFT expansions, one strives to achieveconvergence in describing physical observables by goingto higher and higher orders of expansion. This is bestillustrated by the so-called Lepage plots [19, 23], wherefor theories cut at different orders, relative errors of ob-servables are plotted as functions of energy. In nuclearEDF methods, this kind of convergence tests were neverperformed—simply because the functionals beyond thesecond order (NLO) of the standard Skyrme type hadnever been constructed or studied. The present studyconstitutes the first step towards this goal.Our paper is organized as follows. In Sec. II we de-fine basic building blocks for our construction, and thenwe construct local densities up to N LO. In Sec. III weconstruct terms in the EDF up to N LO and evaluateconstraints imposed by the Galilean and gauge symme-tries. In Sec. IV we derive results for the case of con-served spherical, space-inversion, and time-reversal sym-metries. After formulating conclusions of the presentstudy in Sec. V, in Appendix A we discuss general sym-metry properties of the energy density, in Appendix B wepresent details of the adopted choice of the phase conven-tion, and in Appendix C we list results pertaining to theGalilean and gauge symmetries.
II. CONSTRUCTION OF LOCAL DENSITIESA. Building blocks
Let ρ ( r σ, r ′ σ ′ ) denote the one-body density matrix inspace-spin coordinates. In what follows, in order to sim-plify the notation, we omit the isospin degree of freedom,because in the particle-hole channel all densities appearin the isoscalar and isovector forms [24], and general-ization to proton-neutron systems does not present anyproblem. Within this assumption, the EDF we considerhas the form: E = Z d r H E ( r ) , (1) where the energy density H E ( r ) can be represented as asum of the kinetic and potential energies, H E ( r ) = ~ m τ + H ( r ) . (2)In the present study, we focus on the potential energydensity H ( r ) only.First, using the Pauli matrices σ a , where index a = { x, y, z } enumerates the Cartesian components of a vec-tor, the density matrix is separated into the standardscalar and vector parts [25], ρ ( r σ, r ′ σ ′ ) = ρ ( r , r ′ ) δ σσ ′ + X a h σ | σ a | σ ′ i s a ( r , r ′ ) , (3)where ρ ( r , r ′ ) = X σ ρ ( r σ, r ′ σ ) , (4) s ( r , r ′ ) = X σσ ′ ρ ( r σ, r ′ σ ′ ) h σ ′ | σ | σ i . (5)These two nonlocal densities will be used as buildingblocks of the functional together with the derivative op-erator ∇ and the relative momentum operator k , k = 12 i ( ∇ − ∇ ′ ) . (6)To most easily satisfy the constraints imposed by therotational invariance, in our method, the building blocksare represented as spherical tensor operators [26], i.e., ρ λµ ( r , r ′ ) for λ = 0 and s λµ ( r , r ′ ), ∇ λµ , and k λµ for λ =1. In this notation, λ is the rank of the tensor, and µ = − λ, . . . , + λ is its tensor component. In the present studywe use the following definitions of the building blocks inthe spherical representation: ρ ( r , r ′ ) = ρ ( r , r ′ ) , (7) s ,µ = {− , , } ( r , r ′ ) = − i n √ ( s x ( r , r ′ ) − is y ( r , r ′ )) , s z ( r , r ′ ) , − √ ( s x ( r , r ′ ) + is y ( r , r ′ )) o , (8) ∇ ,µ = {− , , } = − i n √ ( ∇ x − i ∇ y ) , ∇ z , − √ ( ∇ x + i ∇ y ) o , (9) k ,µ = {− , , } = − i n √ ( k x − ik y ) , k z , − √ ( k x + ik y ) o . (10)In what follows, we most often omit indices and argu-ments of these spherical tensors and we simply write ρ , s , ∇ , and k to lighten the notation. In principle, arbitrary phase factors could be used infront of the spherical tensors. In Appendix B, we discusspossible choices of such phase conventions, and determinethe particular ones selected in Eqs. (7)–(10). These phaseconventions, which are not the standard ones, are usedthroughout the paper and define the phase properties ofall other objects that we construct by using the buildingblocks above. B. Higher-order derivative operators
We begin by constructing all possible higher-order andhigher-rank tensor operators from powers of the deriva-tive ∇ µ , where µ = − , , +1 are the spin-projectioncomponents of the vector (rank-1) operator ∇ . It is ob-vious that all possible n th-order powers of the derivativecan be written as sums of terms ∇ µ . . . ∇ µ n . There-fore, any ( n + 1)th-order power is simply obtained bymultiplying some n th-order power by a sum of ∇ µ op-erators. Then, powers of a given rank can be obtainediteratively by vector coupling.In the second order, the two nabla operators can becoupled to angular momenta 0 and 2. The coupling to an-gular momentum 0, [ ∇∇ ] = ∆ / √
3, corresponds to theLaplacian operator. Furthermore, the coupling to angu-lar momentum 2, [ ∇∇ ] , gives the second-order, rank-2derivative operator. The rank-1 coupling, [ ∇∇ ] = 0,vanishes because the derivatives commute. Similarly, ineach one higher order, a rank- L symmetric operator canbe coupled with ∇ only to L − L + 1. Hence, all the n th-order powers have the form of ∆ ( n − L ) / multipliedby the L th-order rank- L (stretched) coupled operatorsfor L = n, n − , . . . , (1)0. Then, up to N LO, one ob-tains 16 different operators D nL listed in Table I. Anyarbitrary tensor formed by coupled operators ∇ can al-ways be rewritten as a sum of operators D nL through therepeated use of the 6 j symbols.Exactly in the same way, we define 16 different oper-ators K nL , which are spherical tensors built of the rel-ative momentum operators k coupled up to N LO, i.e.,for n ≤ L ≤
6. In the remainder of this section, weonly discuss operators D nL , while all the results mutatismutandis also pertain to operators K nL .The stretched coupled operators D LL for n = L , D LL = [ ∇ . . . [ ∇ [ ∇∇ ] ] . . . ] L , (11)play a central role in our derivations below. They corre-spond to irreducible symmetric traceless Cartesian ten-sors built of the derivative ∇ . They have 2 L + 1 tensorcomponents D LLM numbered by the quantum number M = − L, . . . , L that we most often do not show belowexplicitly. Moreover, since terms in the EDF up to N LOdepend only on operators D LL and K LL up to fourthorder, L ≤
4, see Sec. III A, below we do not discussstretched coupled operators of fifth or sixth orders.Equivalently, derivative operators D LL can be writtenin the Cartesian representation, in which their compo-nents are numbered by L Cartesian indices, D LL,a ...a L , a i = x, y, z . The order of these indices does not matter TABLE I: Derivative operators D nL up to N LO as expressedthrough spherical tensor representation of the operator ∇ .No. tensor D nL order n rank L ∇ ∇∇ ] ∇∇ ] ∇∇ ] ∇ ∇ [ ∇∇ ] ] ∇∇ ] ∇∇ ] [ ∇∇ ] ∇ [ ∇ [ ∇∇ ] ] ] ∇∇ ] ∇ ∇∇ ] [ ∇ [ ∇∇ ] ] ∇ [ ∇ [ ∇ [ ∇∇ ] ] ] ] ∇∇ ] ∇∇ ] [ ∇∇ ] ∇∇ ] [ ∇ [ ∇ [ ∇∇ ] ] ] ∇ [ ∇ [ ∇ [ ∇ [ ∇∇ ] ] ] ] ] (totally symmetric tensors) and all traces vanish, X a D LL,aaa ...a L = 0 . (12)The Cartesian components D LL,a ...a L can be calcu-lated by using the detracer operator defined in Sec. 5of Ref. [27]. Up to fourth order they read: D = 1 , (13) D ,a = ∇ a , (14) D ,a a = ∇ a ∇ a − δ a a ∆ , (15) D ,a a a = ∇ a ∇ a ∇ a − ∆ (cid:16) ∇ a δ a a + ∇ a δ a a + ∇ a δ a a (cid:17) , (16) D ,a a a a = ∇ a ∇ a ∇ a ∇ a − ∆ (cid:16) ∇ a ∇ a δ a a + ∇ a ∇ a δ a a + ∇ a ∇ a δ a a + ∇ a ∇ a δ a a + ∇ a ∇ a δ a a + ∇ a ∇ a δ a a (cid:17) + · ∆ (cid:16) δ a a δ a a + δ a a δ a a + δ a a δ a a (cid:17) . (17) TABLE II: Spherical components of the derivative operators D nLM expressed through the Cartesian derivatives. Expres-sions for negative components can be obtained as D nL, − M =( − L − M D ∗ nLM , see Eqs. (B20) and (B22). D nLM Cartesian derivatives D = − i∂ z D = i √ ( ∂ x + i∂ y ) D = √ ∆ D = √ ` ∂ x + ∂ y − ∂ z ´ D = ( ∂ x + i∂ y ) ∂ z D = − ( ∂ x + i∂ y ) D = i √ ∂ z ` − ∂ x − ∂ y + 2 ∂ z ´ D = i q ( ∂ x + i∂ y ) ` ∂ x + ∂ y − ∂ z ´ D = i √ ∂ x + i∂ y ) ∂ z D = − i √ ( ∂ x + i∂ y ) D = √ ` ∂ x + 6 ` ∂ y − ∂ z ´ ∂ x + 3 ∂ y + 8 ∂ z − ∂ y ∂ z ´ D = √ ( ∂ x + i∂ y ) ∂ z ` ∂ x + 3 ∂ y − ∂ z ´ D = − √ ( ∂ x + i∂ y ) ` ∂ x + ∂ y − ∂ z ´ D = − √ ( ∂ x + i∂ y ) ∂ z D = ( ∂ x + i∂ y ) We note here in passing that we could have equallywell used the Cartesian derivative operators with tracesnot subtracted out, i.e., D = 1 , (18) D ,a = ∇ a , (19) D ,a a = ∇ a ∇ a , (20) D ,a a a = ∇ a ∇ a ∇ a , (21) D ,a a a a = ∇ a ∇ a ∇ a ∇ a . (22)Representations (13)–(17) and (18)–(22) are equivalentin the sense that each operator D LL,a ...a L is evidentlya linear combination of operators ∆ ( L − L ′ ) / D L ′ L ′ ,a ...a L ′ for L ′ = L , L − , . . . , (1)0.In principle, below one could replace the spherical rep-resentations of derivative operators shown in Table I bytheir Cartesian counterparts (13)–(17) or (18)–(22), andwork entirely in the Cartesian representation. However,in our opinion, the use of the spherical representation issuperior and more economical. Moreover, whenever cal-culation of the Cartesian derivatives is more suitable, wemay express spherical components of the derivative op-erators through the Cartesian derivatives, as shown inTable II. An example of using the Cartesian representa-tion (18)–(22) is given in Sec. IV. TABLE III: Local primary densities (23) up to N LO builtfrom the scalar nonlocal density ρ ( r , r ′ ) ( v = 0). To simplifythe notation the limit of r ′ = r is not shown explicitly. Stars( ⋆ ) mark densities that enter the EDF up to N LO. Bullets ( • )mark densities that enter the EDF up to N LO for conservedspherical, space-inversion, and time-reversal symmetries, seeSec. IV. The last two columns show the T and P paritiesdefined in Eqs. (25) and (26), respectively. In addition, thetime-even densities are marked by using the bold-face font.No. ρ nLvJ = density n L v J T P ⋆ • ρ = [ ρ ] ⋆ ρ = [ kρ ] − − ⋆ • ρ = [[ kk ] ρ ] ⋆ • ρ = [[ kk ] ρ ] ⋆ ρ = [[ kk ] kρ ] − − ⋆ ρ = [[ k [ kk ] ] ρ ] − − ⋆ • ρ = [[ kk ] ρ ] ⋆ • ρ = [[ kk ] [ kk ] ρ ] ρ = [[ k [ k [ kk ] ] ] ρ ] ⋆ ρ = [[ kk ] kρ ] − − ρ = [[ kk ] [ k [ kk ] ] ρ ] − − ρ = [[ k [ k [ k [ kk ] ] ] ] ρ ] − − ⋆ • ρ = [[ kk ] ρ ] ρ = [[ kk ] [ kk ] ρ ] ρ = [[ kk ] [ k [ k [ kk ] ] ] ρ ] ρ = [[ k [ k [ k [ k [ kk ] ] ] ] ] ρ ] C. Local densities
Local densities are formed by acting several times onthe scalar and vector nonlocal densities with the relativemomentum operator k and taking the limit of r ′ = r .Using the spherical representation, the possible coupled k -tensors (10) (up to sixth order in derivatives) K nL arethose given in Table I (replacing ∇ with k ).Acting with K nL on the scalar nonlocal density ρ ( r , r ′ )gives 16 different local densities up to N LO (one forevery term in Table I). They are listed in Table III.When acting with K nL on the vector nonlocal densities s ( r , r ′ ), one has to construct all possible ways of couplingthe k -tensors with the vector density. Obviously, each ofthe 4 scalar ( L = 0) derivative operators gives one localdensity, while each of the 12 non-scalar ( L >
0) derivativeoperators gives three local densities. Altogether, from thevector density one obtains 40 local densities up to N LO.They are listed in Table IV.Finally, all local densities can be denoted by four inte-gers nLvJ as ρ nLvJ ( r ) = { [ K nL ρ v ( r , r ′ )] J } r ′ = r , (23)where the n th-order and rank- L relative derivative oper-ator K nL acts on the scalar ( v = 0) or vector ( v = 1)nonlocal density, and ranks L and v are then vector cou-pled to J . We call these local densities primary densities. TABLE IV: Same as in Table III but for densities built fromthe vector nonlocal density s ( r , r ′ ) ( v = 1).No. ρ nLvJ = density n L v J T P ⋆ ρ = [ s ] − ⋆ ρ = [ ks ] − ⋆ • ρ = [ ks ] − ⋆ ρ = [ ks ] − ⋆ ρ = [[ kk ] s ] − ⋆ ρ = [[ kk ] s ] − ⋆ ρ = [[ kk ] s ] − ⋆ ρ = [[ kk ] s ] − ⋆ ρ = [[ kk ] ks ] − ⋆ • ρ = [[ kk ] ks ] − ⋆ ρ = [[ kk ] ks ] − ⋆ ρ = [[ k [ kk ] ] s ] − ⋆ • ρ = [[ k [ kk ] ] s ] − ⋆ ρ = [[ k [ kk ] ] s ] − ⋆ ρ = [[ kk ] s ] − ⋆ ρ = [[ kk ] [ kk ] s ] − ⋆ ρ = [[ kk ] [ kk ] s ] − ⋆ ρ = [[ kk ] [ kk ] s ] − ⋆ ρ = [[ k [ k [ kk ] ] ] s ] − ρ = [[ k [ k [ kk ] ] ] s ] − ρ = [[ k [ k [ kk ] ] ] s ] − ρ nLvJ = density n L v J T P ⋆ ρ = [[ kk ] ks ] − ⋆ • ρ = [[ kk ] ks ] − ⋆ ρ = [[ kk ] ks ] − ⋆ ρ = [[ kk ] [ k [ kk ] ] s ] − ρ = [[ kk ] [ k [ kk ] ] s ] − ρ = [[ kk ] [ k [ kk ] ] s ] − ρ = [[ k [ k [ k [ kk ] ] ] ] s ] − ρ = [[ k [ k [ k [ kk ] ] ] ] s ] − ρ = [[ k [ k [ k [ kk ] ] ] ] s ] − ⋆ ρ = [[ kk ] s ] − ⋆ ρ = [[ kk ] [ kk ] s ] − ρ = [[ kk ] [ kk ] s ] − ρ = [[ kk ] [ kk ] s ] − ρ = [[ kk ] [ k [ k [ kk ] ] ] s ] − ρ = [[ kk ] [ k [ k [ kk ] ] ] s ] − ρ = [[ kk ] [ k [ k [ kk ] ] ] s ] − ρ = [[ k [ k [ k [ k [ kk ] ] ] ] ] s ] − ρ = [[ k [ k [ k [ k [ kk ] ] ] ] ] s ] − ρ = [[ k [ k [ k [ k [ kk ] ] ] ] ] s ] − The tensor components corresponding to the total rank J are not explicitly shown.One can also act on each of the local densities withderivative operators D mI of Table I, and then coupleranks I and J to the total rank Q , i.e., ρ mI,nLvJ,Q ( r ) = [ D mI ρ nLvJ ( r )] Q . (24)For m >
0, we call these local densities secondary den-sities. We do not explicitly list them, because they canbe obtained in a straightforward way from the primarydensities corresponding to m = 0, ρ nLvJ = ρ ,InLvJ,J ,which are listed in Tables III and IV.In Tables III and IV, for completeness we also show thetime-reversal ( T ) and space-inversion ( P ) parities definedas, T = ( − n + v , (25) P = ( − n . (26)These definitions are based on the analysis of symmetryproperties, which we present in Appendix A. To bettervisualize the time-reversal properties of the local den-sities, in Tables III and IV the time-even densities areshown in bold face.Local densities constructed above are complex. Tak-ing the complex conjugations gives relations derived inAppendix B: ρ ∗ mI,nLvJ,Q,M = ( − Q − M ρ mI,nLvJ,Q, − M , (27)where the tensor components, denoted M , are shownexplicitly. These relations allow for expressing positivetensor components through negative ones or vice versa .Therefore, complete information is contained in non-negative or non-positive tensor components only. The M = 0 components are either real (for even Q ) or imagi-nary (for odd Q ), and hence 2 Q + 1 real functions alwayssuffice to describe a given local density of rank Q . More-over, all scalar densities are real, which was the basis ofchoosing this particular phase convention, as describedin Appendix B. III. CONSTRUCTION OF THE ENERGYDENSITYA. Terms in the energy density
Terms in the EDF we construct here are required tobe quadratic in densities, invariant with respect to timereversal (Sec. A 1), and covariant with respect to spaceinversion and rotations (Sec. A 2). All terms up to theN LO order in derivatives fulfilling these restrictions areconstructed below.Using notation of Eq. (24), a general term in the energydensity can be written in the following form, T m ′ I ′ ,n ′ L ′ v ′ J ′ mI,nLvJ,Q ( r ) = [ ρ m ′ I ′ ,n ′ L ′ v ′ J ′ ,Q ( r ) ρ mI,nLvJ,Q ( r )] , (28) TABLE V: Numbers of terms defined in Eq. (28) of differentorders in the EDF up to N LO. Numbers of terms dependingon the time-even and time-odd densities are given separately.The last two columns give numbers of terms when the Galileanor gauge invariance is assumed, respectively, see Sec. III B. Totake into account both isospin channels, the numbers of termsshould be multiplied by a factor of two.order T-even T-odd Total Galilean Gauge0 1 1 2 2 22 8 10 18 12 124 53 61 114 45 296 250 274 524 129 54N LO 312 346 658 188 97
Order in derivatives N r . o f t e r m s o f d i ff e r e n t o r d e r Eq. (30)Eq. (30), Galilean inv.Eq. (30), Gauge inv.Eq. (28)Eq. (28), Galilean inv.Eq. (28), Gauge inv.
FIG. 1: (Color online) Numbers of terms (28) and (30) shownin Tables V and VI, respectively, plotted in logarithmic scaleas a function of the order in derivatives. where both densities must have the same rank Q to becoupled to a scalar. Moreover, their time-reversal andspace-inversion parities ( T and P ) must be the same.Again, at N LO only terms with m ′ + n ′ + m + n ≤ H ( r ) = X m ′ I ′ ,n ′ L ′ v ′ J ′ mI,nLvJ,Q C m ′ I ′ ,n ′ L ′ v ′ J ′ mI,nLvJ,Q T m ′ I ′ ,n ′ L ′ v ′ J ′ mI,nLvJ,Q ( r ) , (29)where C m ′ I ′ ,n ′ L ′ v ′ J ′ mI,nLvJ,Q are coupling constants and the sum-mation runs over all allowed indices.Had we considered the case of coupling constants de-pending on density, all terms in Eq. (28) would have beenindependent of one another (up to a possible exchange ofthe two densities). Table V lists numbers of such inde-pendent terms, and they are also plotted in Fig. 1.In the present study, we concentrate on the case ofdensity-independent coupling constants, in which caseone can perform integrations by parts, so that the deriva-tive operators D m ′ I ′ are transferred from one density to TABLE VI: Same as in Table V, but for numbers of termsdefined in Eq. (30).order T-even T-odd Total Galilean Gauge0 1 1 2 2 22 6 6 12 7 74 22 23 45 15 66 64 65 129 26 6N LO 93 95 188 50 21 the other. That this can always be done is obvious bythe fact that the coupled derivative operators D m ′ I ′ canalways be expressed as sums of products of uncoupledderivatives ∇ µ or ∇ a . As a result of the integration byparts, the integral (28) can now be written as a sum ofterms, where each term has the form: T n ′ L ′ v ′ J ′ mI,nLvJ ( r ) = [ ρ n ′ L ′ v ′ J ′ ( r )[ D mI ρ nLvJ ( r )] J ′ ] , (30)where ranks I and J are coupled to J ′ . Here, at N LO(i) only terms with n ′ + m + n ≤ T , and(iii) their space-inversion parities must differ by factors( − I . Finally, in order to avoid double-counting onetakes only terms with n ′ < n , and for n ′ = n only thosewith L ′ < L , and for L ′ = L only those with v ′ < v , andfor v ′ = v only those with J ′ ≤ J . Then, the total energydensity reads H ( r ) = X n ′ L ′ v ′ J ′ mI,nLvJ,J ′ C n ′ L ′ v ′ J ′ mI,nLvJ T n ′ L ′ v ′ J ′ mI,nLvJ ( r ) , (31)where C n ′ L ′ v ′ J ′ mI,nLvJ are coupling constants and the summa-tion again runs over all allowed indices. As we did for thelocal densities above, to better visualize the time-reversalcharacteristics of terms in the EDF, the coupling con-stants C n ′ L ′ v ′ J ′ mI,nLvJ corresponding to terms that depend ontime-even densities are shown in bold face.Based on the results obtained in Secs. A 1 and A 2, andon Eqs. (25) and (26), we see that time-reversal invari-ance and space-inversion covariance require that( − n ′ + v ′ + n + v = 1 , (32)( − n ′ + m + n = 1 , (33)respectively. This means that integers v ′ + v , n ′ + n ,and m must be simultaneously either even or odd. Thenumbers of all such allowed terms are given in Table VIand plotted in Fig. 1. The space-inversion covariance (33)requires that for all terms, the total orders in derivativesare even numbers, which defines our classification of theEDF up to LO (0), NLO (2), NNLO (4), and N LO (6).In Appendix B, we presented terms in the EDF up toNLO, i.e., for zero and second orders, see Table XXII.The EDF at NLO is exactly equivalent to the standard
TABLE VII: Terms in the EDF (30) that are of fourth or-der, depend on time-even densities, and are built from thescalar nonlocal density ρ ( r , r ′ ). Coupling constants corre-sponding to terms that depend on time-even densities aremarked by using the bold-face font. Bullets ( • ) mark cou-pling constants corresponding to terms that do not vanish forconserved spherical, space-inversion, and time-reversal sym-metries, see Sec. IV.No. C n ′ L ′ v ′ J ′ mI,nLvJ ρ n ′ L ′ v ′ J ′ D mI ρ nLvJ • C , [ ρ ] [ ∇∇ ] [ ρ ] • C , [ ρ ] [ ∇∇ ] [[ kk ] ρ ] • C , [ ρ ] [ ∇∇ ] [[ kk ] ρ ] • C , [ ρ ] kk ] ρ ] • C , [[ kk ] ρ ] kk ] ρ ] • C , [[ kk ] ρ ] kk ] ρ ] TABLE VIII: Same as in Table VII but for terms that arebuilt from the vector nonlocal density s ( r , r ′ ).No. C n ′ L ′ v ′ J ′ mI,nLvJ ρ n ′ L ′ v ′ J ′ D mI ρ nLvJ C , [ ks ] [ ∇∇ ] [ ks ] C , [ ks ] [ ∇∇ ] [ ks ] C , [ ks ] kk ] ks ] • C , [ ks ] [ ∇∇ ] [ ks ] • C , [ ks ] [ ∇∇ ] [ ks ] C , [ ks ] [ ∇∇ ] [ ks ] • C , [ ks ] kk ] ks ] C , [ ks ] [ ∇∇ ] [ ks ] C , [ ks ] [ ∇∇ ] [ ks ] C , [ ks ] kk ] ks ] C , [ ks ] k [ kk ] ] s ] Skyrme functional [28, 29], generalized to include alltime-odd terms [24, 25, 30]. In both representations thefunctional depends, in general, on 14 coupling constants,and both sets are related by simple expressions given inEqs. (B6)–(B19).In Tables VII–XVIII, we list all 45 and 129 terms inthe EDF that are of fourth and sixth order, respectively.Together with 14 terms at NLO, listed in Table XXII,this constitutes the full list of 188 terms in the EDF atN LO.After the complete list of terms in the EDF at N LOis constructed, one can check that not all of the localdensities listed in Tables III and IV appear in the finalEDF at N LO. This is so, because it is not possible tocouple all these densities to scalars, and simultaneouslyfulfill conditions (32) and (33), without obtaining morethan total sixth order in derivatives. It turns out thatout of the 56 local densities at N LO, which are listedin Tables III and IV, only 35 occur in the final EDF
TABLE IX: Same as in Table VII but for terms that are builtfrom the scalar nonlocal density ρ ( r , r ′ ) and vector nonlocaldensity s ( r , r ′ ).No. C n ′ L ′ v ′ J ′ mI,nLvJ ρ n ′ L ′ v ′ J ′ D mI ρ nLvJ • C , [ ρ ] [ ∇∇ ] ∇ [ ks ] • C , [ ρ ] ∇ [[ kk ] ks ] • C , [[ kk ] ρ ] ∇ [ ks ] • C , [[ kk ] ρ ] ∇ [ ks ] C , [[ kk ] ρ ] ∇ [ ks ] TABLE X: Terms in the EDF (30) that are of fourth order,depend on time-odd densities, and are built from the scalarnonlocal density ρ ( r , r ′ ).No. C n ′ L ′ v ′ J ′ mI,nLvJ ρ n ′ L ′ v ′ J ′ D mI ρ nLvJ C , [ kρ ] [ ∇∇ ] [ kρ ] C , [ kρ ] [ ∇∇ ] [ kρ ] C , [ kρ ] kk ] kρ ] at N LO. In Tables III and IV such densities are markedwith stars ( ⋆ ). Table XIX gives their numbers determinedseparately at each order. TABLE XI: Same as in Table X but for terms that are builtfrom the vector nonlocal density s ( r , r ′ ).No. C n ′ L ′ v ′ J ′ mI,nLvJ ρ n ′ L ′ v ′ J ′ D mI ρ nLvJ C , [ s ] [ ∇∇ ] [ s ] C , [ s ] [ ∇∇ ] [ ∇∇ ] [ s ] C , [ s ] [ ∇∇ ] [[ kk ] s ] C , [ s ] [ ∇∇ ] [[ kk ] s ] C , [ s ] [ ∇∇ ] [[ kk ] s ] C , [ s ] [ ∇∇ ] [[ kk ] s ] C , [ s ] [ ∇∇ ] [[ kk ] s ] C , [ s ] [ ∇∇ ] [[ kk ] s ] C , [ s ] kk ] s ] C , [ s ] kk ] [ kk ] s ] C , [[ kk ] s ] kk ] s ] C , [[ kk ] s ] kk ] s ] C , [[ kk ] s ] kk ] s ] C , [[ kk ] s ] kk ] s ] C , [[ kk ] s ] kk ] s ] TABLE XII: Same as in Table X but for terms that are builtfrom the scalar nonlocal density ρ ( r , r ′ ) and vector nonlocaldensity s ( r , r ′ ).No. C n ′ L ′ v ′ J ′ mI,nLvJ ρ n ′ L ′ v ′ J ′ D mI ρ nLvJ C , [ kρ ] [ ∇∇ ] ∇ [ s ] C , [ kρ ] ∇ [[ kk ] s ] C , [ kρ ] ∇ [[ kk ] s ] C , [ kρ ] ∇ [[ kk ] s ] C , [[ kk ] kρ ] ∇ [ s ] TABLE XIII: Terms in the EDF (30) that are of sixth or-der, depend on time-even densities, and are built from thescalar nonlocal density ρ ( r , r ′ ). Coupling constants corre-sponding to terms that depend on time-even densities aremarked by using the bold-face font. Bullets ( • ) mark cou-pling constants corresponding to terms that do not vanish forconserved spherical, space-inversion, and time-reversal sym-metries, see Sec. IV.No. C n ′ L ′ v ′ J ′ mI,nLvJ ρ n ′ L ′ v ′ J ′ D mI ρ nLvJ • C , [ ρ ] [ ∇∇ ] [ ρ ] • C , [ ρ ] [ ∇∇ ] [[ kk ] ρ ] • C , [ ρ ] [ ∇∇ ] [ ∇∇ ] [[ kk ] ρ ] • C , [ ρ ] [ ∇∇ ] [[ kk ] ρ ] • C , [ ρ ] [ ∇∇ ] [[ kk ] [ kk ] ρ ] • C , [ ρ ] kk ] ρ ] • C , [[ kk ] ρ ] [ ∇∇ ] [[ kk ] ρ ] • C , [[ kk ] ρ ] [ ∇∇ ] [[ kk ] ρ ] • C , [[ kk ] ρ ] kk ] ρ ] • C , [[ kk ] ρ ] [ ∇∇ ] [[ kk ] ρ ] • C , [[ kk ] ρ ] [ ∇∇ ] [[ kk ] ρ ] • C , [[ kk ] ρ ] kk ] [ kk ] ρ ] B. Galilean and gauge invariance
In the previous Sec. III A, the functional has been re-quired to be consistent with time reversal invariance,invariance under space reflections and rotational invari-ance. These constraints arise from symmetries of theNN interaction, see e.g. Refs. [31, 32]. Derivations abovewere much easier to perform in a general form, withoutimposing any other additional symmetry conditions. Inthis section, we treat such additional constraints comingfrom imposing the Galilean and gauge invariance.Assumption of the Galilean instead of Lorentz invari-ance goes hand in hand with using the Schr¨odinger equa-tion as a starting point and relies on the assumptionthat relativistic effects are negligible. This symmetry en-sures that the collective translational mass, calculatedwithin the time-dependent HF or random-phase approx-imations, is correctly equal to the total mass, M = Am . TABLE XIV: Same as in Table XIII but for terms that arebuilt from the vector nonlocal density s ( r , r ′ ).No. C n ′ L ′ v ′ J ′ mI,nLvJ ρ n ′ L ′ v ′ J ′ D mI ρ nLvJ C , [ ks ] [ ∇∇ ] [ ks ] C , [ ks ] [ ∇∇ ] [ ∇∇ ] [ ks ] C , [ ks ] [ ∇∇ ] [[ kk ] ks ] C , [ ks ] [ ∇∇ ] [[ kk ] ks ] C , [ ks ] [ ∇∇ ] [[ k [ kk ] ] s ] C , [ ks ] kk ] ks ] • C , [ ks ] [ ∇∇ ] [ ks ] • C , [ ks ] [ ∇∇ ] [ ∇∇ ] [ ks ] C , [ ks ] [ ∇∇ ] [ ∇∇ ] [ ks ] • C , [ ks ] [ ∇∇ ] [[ kk ] ks ] • C , [ ks ] [ ∇∇ ] [[ kk ] ks ] C , [ ks ] [ ∇∇ ] [[ kk ] ks ] C , [ ks ] [ ∇∇ ] [[ k [ kk ] ] s ] • C , [ ks ] [ ∇∇ ] [[ k [ kk ] ] s ] • C , [ ks ] kk ] ks ] Therefore, in principle, the Galilean invariance should al-ways be imposed. However, in many phenomenologicalapproaches, like the non-interacting or interacting shellmodel, the Galilean symmetry is not considered, becausethe translational motion is not within the scope of suchmodels. The question of whether the Galilean symmetrymust be imposed in phenomenological models is not yetresolved, and in the present study we keep this questionopen.For a local interaction, v ( r ′ , r ′ , r , r ) = δ ( r ′ − r ) δ ( r ′ − r ) v ( r , r ), the HF interaction energy is in-variant with respect to the local gauge [36]. Therefore,for the total energy (1) obtained within the EDF method,one may also consider constraints resulting from assum-ing local gauge invariance. As mentioned, this symmetryis only fulfilled when the forces involved are local. Anexample of a local approximation is the well known lo-cal one pion-exchange (OPE) potential, which is only anapproximate representation of the correct nonlocal OPEFeynman amplitude [33]. This approximation is goodas long as the relative momenta of interacting particlesare about the same in the initial and final states (seeFig. 10 in Ref. [33]). Some of the fitted NN potentialslike the Argonne V , Nijm-II, and Reid93 use this localapproximation while others (CD-Bonn) use the full non-local OPE amplitude [33]. The most important nonlocalterm is the two-body spin orbit interaction [32], whichviolates the assumption of gauge invariance. However,in this case a gauge invariant spin-orbit term (used inthe Skyrme and Gogny forces) can be obtained in theshort-range limit [32, 34, 36].For the EDF derived in this work it is, however, thesymmetries of effective forces rather than bare forces thatshould be considered. One of the methods to obtain an TABLE XIV: continued.No. C n ′ L ′ v ′ J ′ mI,nLvJ ρ n ′ L ′ v ′ J ′ D mI ρ nLvJ C , [ ks ] [ ∇∇ ] [ ks ] C , [ ks ] [ ∇∇ ] [ ∇∇ ] [ ks ] C , [ ks ] [ ∇ [ ∇ [ ∇∇ ] ] ] [ ks ] C , [ ks ] [ ∇∇ ] [[ kk ] ks ] C , [ ks ] [ ∇∇ ] [[ kk ] ks ] C , [ ks ] [ ∇∇ ] [[ kk ] ks ] C , [ ks ] [ ∇∇ ] [[ kk ] ks ] C , [ ks ] [ ∇∇ ] [[ k [ kk ] ] s ] C , [ ks ] [ ∇∇ ] [[ k [ kk ] ] s ] C , [ ks ] [ ∇∇ ] [[ k [ kk ] ] s ] C , [ ks ] [ ∇∇ ] [[ k [ kk ] ] s ] C , [ ks ] kk ] ks ] C , [ ks ] kk ] [ k [ kk ] ] s ] C , [[ kk ] ks ] kk ] ks ] • C , [[ kk ] ks ] kk ] ks ] C , [[ kk ] ks ] kk ] ks ] C , [[ kk ] ks ] k [ kk ] ] s ] C , [[ k [ kk ] ] s ] k [ kk ] ] s ] • C , [[ k [ kk ] ] s ] k [ kk ] ] s ] C , [[ k [ kk ] ] s ] k [ kk ] ] s ] TABLE XV: Same as in Table XIII but for terms that are builtfrom the scalar nonlocal density ρ ( r , r ′ ) and vector nonlocaldensity s ( r , r ′ ).No. C n ′ L ′ v ′ J ′ mI,nLvJ ρ n ′ L ′ v ′ J ′ D mI ρ nLvJ • C , [ ρ ] [ ∇∇ ] ∇ [ ks ] • C , [ ρ ] [ ∇∇ ] ∇ [[ kk ] ks ] • C , [ ρ ] [ ∇ [ ∇∇ ] ] [[ k [ kk ] ] s ] • C , [ ρ ] ∇ [[ kk ] ks ] • C , [[ kk ] ρ ] [ ∇∇ ] ∇ [ ks ] • C , [[ kk ] ρ ] ∇ [[ kk ] ks ] • C , [[ kk ] ρ ] [ ∇∇ ] ∇ [ ks ] • C , [[ kk ] ρ ] [ ∇ [ ∇∇ ] ] [ ks ] C , [[ kk ] ρ ] [ ∇∇ ] ∇ [ ks ] C , [[ kk ] ρ ] [ ∇ [ ∇∇ ] ] [ ks ] • C , [[ kk ] ρ ] ∇ [[ kk ] ks ] C , [[ kk ] ρ ] ∇ [[ kk ] ks ] C , [[ kk ] ρ ] ∇ [[ k [ kk ] ] s ] • C , [[ kk ] ρ ] ∇ [[ k [ kk ] ] s ] • C , [[ kk ] ρ ] ∇ [ ks ] • C , [[ kk ] [ kk ] ρ ] ∇ [ ks ] C , [[ kk ] [ kk ] ρ ] ∇ [ ks ] TABLE XVI: Terms in the EDF (30) that are of sixth order,depend on time-odd densities, and are built from the scalarnonlocal density ρ ( r , r ′ ).No. C n ′ L ′ v ′ J ′ mI,nLvJ ρ n ′ L ′ v ′ J ′ D mI ρ nLvJ C , [ kρ ] [ ∇∇ ] [ kρ ] C , [ kρ ] [ ∇∇ ] [ ∇∇ ] [ kρ ] C , [ kρ ] [ ∇∇ ] [[ kk ] kρ ] C , [ kρ ] [ ∇∇ ] [[ kk ] kρ ] C , [ kρ ] [ ∇∇ ] [[ k [ kk ] ] ρ ] C , [ kρ ] kk ] kρ ] C , [[ kk ] kρ ] kk ] kρ ] C , [[ k [ kk ] ] ρ ] k [ kk ] ] ρ ] TABLE XVII: Same as in Table XVI but for terms that arebuilt from the vector nonlocal density s ( r , r ′ ).No. C n ′ L ′ v ′ J ′ mI,nLvJ ρ n ′ L ′ v ′ J ′ D mI ρ nLvJ C , [ s ] [ ∇∇ ] [ s ] C , [ s ] [ ∇∇ ] [ ∇∇ ] [ s ] C , [ s ] [ ∇∇ ] [[ kk ] s ] C , [ s ] [ ∇∇ ] [ ∇∇ ] [[ kk ] s ] C , [ s ] [ ∇∇ ] [[ kk ] s ] C , [ s ] [ ∇∇ ] [ ∇∇ ] [[ kk ] s ] C , [ s ] [ ∇∇ ] [ ∇∇ ] [[ kk ] s ] C , [ s ] [ ∇∇ ] [ ∇∇ ] [[ kk ] s ] C , [ s ] [ ∇ [ ∇ [ ∇∇ ] ] ] [[ kk ] s ] C , [ s ] [ ∇∇ ] [[ kk ] s ] C , [ s ] [ ∇∇ ] [[ kk ] s ] C , [ s ] [ ∇∇ ] [[ kk ] [ kk ] s ] C , [ s ] [ ∇∇ ] [[ kk ] [ kk ] s ] C , [ s ] [ ∇∇ ] [[ kk ] [ kk ] s ] C , [ s ] [ ∇∇ ] [[ kk ] [ kk ] s ] C , [ s ] [ ∇∇ ] [[ k [ k [ kk ] ] ] s ] C , [ s ] kk ] s ] C , [ s ] kk ] [ kk ] s ] effective NN force from the bare NN force is the unitarycorrelation operator method (UCOM) [35]. The use ofthe UCOM scheme, however, leads to a nonlocal effectiveinteraction even if the bare interaction would have beena local one.But rather than discussing to which extent gauge sym-metry is conserved or broken in nuclei we aim to providea theoretical framework where different choices can beaccommodated. Because several successful phenomeno-logical forces (e.g. Skyrme and Gogny [36]) are invariantunder the gauge transformation, this symmetry consti-tutes a natural starting point in the search of improvedEDFs. To which extent gauge symmetry is violated for0 TABLE XVII: continued.No. C n ′ L ′ v ′ J ′ mI,nLvJ ρ n ′ L ′ v ′ J ′ D mI ρ nLvJ C , [[ kk ] s ] [ ∇∇ ] [[ kk ] s ] C , [[ kk ] s ] [ ∇∇ ] [[ kk ] s ] C , [[ kk ] s ] [ ∇∇ ] [[ kk ] s ] C , [[ kk ] s ] [ ∇∇ ] [[ kk ] s ] C , [[ kk ] s ] [ ∇∇ ] [[ kk ] s ] C , [[ kk ] s ] [ ∇∇ ] [[ kk ] s ] C , [[ kk ] s ] kk ] s ] C , [[ kk ] s ] kk ] [ kk ] s ] C , [[ kk ] s ] [ ∇∇ ] [[ kk ] s ] C , [[ kk ] s ] [ ∇∇ ] [[ kk ] s ] C , [[ kk ] s ] [ ∇∇ ] [[ kk ] s ] C , [[ kk ] s ] [ ∇∇ ] [[ kk ] s ] C , [[ kk ] s ] kk ] s ] C , [[ kk ] s ] kk ] [ kk ] s ] C , [[ kk ] s ] [ ∇∇ ] [[ kk ] s ] C , [[ kk ] s ] [ ∇∇ ] [[ kk ] s ] C , [[ kk ] s ] [ ∇∇ ] [[ kk ] s ] C , [[ kk ] s ] kk ] [ kk ] s ] C , [[ kk ] s ] [ ∇∇ ] [[ kk ] s ] C , [[ kk ] s ] [ ∇∇ ] [[ kk ] s ] C , [[ kk ] s ] kk ] [ kk ] s ] C , [[ kk ] s ] k [ k [ kk ] ] ] s ] TABLE XVIII: Same as in Table XVI but for terms thatare built from the scalar nonlocal density ρ ( r , r ′ ) and vectornonlocal density s ( r , r ′ ).No. C n ′ L ′ v ′ J ′ mI,nLvJ ρ n ′ L ′ v ′ J ′ D mI ρ nLvJ C , [ kρ ] [ ∇∇ ] ∇ [ s ] C , [ kρ ] [ ∇∇ ] ∇ [[ kk ] s ] C , [ kρ ] [ ∇∇ ] ∇ [[ kk ] s ] C , [ kρ ] [ ∇∇ ] ∇ [[ kk ] s ] C , [ kρ ] [ ∇ [ ∇∇ ] ] [[ kk ] s ] C , [ kρ ] [ ∇ [ ∇∇ ] ] [[ kk ] s ] C , [ kρ ] ∇ [[ kk ] s ] C , [ kρ ] ∇ [[ kk ] [ kk ] s ] C , [ kρ ] ∇ [[ kk ] [ kk ] s ] C , [[ kk ] kρ ] [ ∇∇ ] ∇ [ s ] C , [[ kk ] kρ ] ∇ [[ kk ] s ] C , [[ kk ] kρ ] ∇ [[ kk ] s ] C , [[ kk ] kρ ] ∇ [[ kk ] s ] C , [[ k [ kk ] ] ρ ] [ ∇ [ ∇∇ ] ] [ s ] C , [[ k [ kk ] ] ρ ] ∇ [[ kk ] s ] C , [[ k [ kk ] ] ρ ] ∇ [[ kk ] s ] C , [[ kk ] kρ ] ∇ [ s ] TABLE XIX: Numbers of local densities ρ nLvJ , Eq. (23), ofdifferent orders, which enter into the EDF up to N LO. Num-bers of local densities constructed from the scalar ρ ( r , r ′ ) orvector s ( r , r ′ ) nonlocal densities, and numbers of time-evenand time-odd local densities, are given separately. In Ta-bles III and IV these densities are marked with stars ( ⋆ ).order from ρ from s T-even T-odd total0 1 1 1 1 21 1 3 3 1 42 2 4 2 4 63 2 6 6 2 84 2 5 2 5 75 1 4 4 1 56 1 2 1 2 3total 10 25 19 16 35 effective renormalized interactions is a question that canbe investigated by comparing models using preserved andbroken symmetries (see e.g. Ref. [37]).
1. Local gauge transformations of the nonlocal densities
The gauge-transformed nonlocal densities read [24, 25,36] ρ ′ ( r , r ′ ) = e i [ φ ( r ) − φ ( r ′ ) ] ρ ( r , r ′ ) , (34) s ′ ( r , r ′ ) = e i [ φ ( r ) − φ ( r ′ ) ] s ( r , r ′ ) . (35)Since the local gauge transformations form a U(1) group,invariance with respect to transformations that are of thefirst-order in gauge angles, [1 + iG ] ρ ( r , r ′ ), where G ( r , r ′ ) = φ ( r ) − φ ( r ′ ) , (36)is enough to ensure full gauge invariance. By Taylor ex-panding the exponential functions in eq. 34 and 35 afterthey are inserted in the functional one may, of course,also prove this fact explicitly.One specific type of gauge transformation is theGalilean transformation, for which the gauge angles de-pend linearly on positions, i.e., G ( r , r ′ ) = p · ( r − r ′ ) / ~ ,and which corresponds to a transformation to a referenceframe that moves with velocity p /m . For this transfor-mation, only first-order derivatives of G survive, whichmakes Galilean invariance less restrictive than the fullgauge invariance.
2. Local gauge transformations of the local densities
Let indices β, γ, . . . = 1 , . . . ,
35 label primary local den-sities ρ nLvJ (23) listed with stars ( ⋆ ) in Tables III andIV, which enter the EDF at N LO, as shown in Eqs. (30)1and (31). Using this notation, the linearized gauge trans-formation of one of the local densities can be written as ρ ′ β ( r ) = { [ K nL (1 + iG ( r , r ′ )) ρ v ( r , r ′ )] J } r = r ′ = ρ β ( r ) + { [ K nL iG ( r , r ′ ) ρ v ( r , r ′ )] J } r = r ′ = ρ β ( r ) + ρ Gβ ( r ) , (37)where the first term is the untransformed local densityand the second term is the part affected by the gaugetransformation.As an illustration, let us begin by considering the sim-pler case of Galilean transformation, and look at the termwith n = 2, where only two relative momentum operators k appear. Operator k can be written as k ρ + k G , with thefirst term acting only on ρ v ( r , r ′ ) and the second termacting only on G ( r , r ′ ). Then we have. K L = [ kk ] L = [ k ρ k ρ ] L + 2 [ k ρ k G ] L + [ k G k G ] L . (38)When this is inserted into the expression for ρ Gβ , the lastterm can be dropped since only the first-order derivativesof G ( r , r ′ ) survive for the Galilean transformation, andthe first term disappears when one takes the limit of r = r ′ since G ( r , r ) = 0. Thus in this case we obtain ρ Gβ ( r ) = i (cid:8) (cid:2) [ k ρ k G ] l G ( r , r ′ ) ρ v ( r , r ′ ) (cid:3) J , (cid:9) r = r ′ = i X J ′ c J ′ (cid:2)(cid:8) [ k ρ ρ v ] J ′ (cid:9) r = r ′ ( ∇ φ )( r ) (cid:3) J , (39)where the second equation results from the vector re-coupling (note that G ( r , r ′ ) is a scalar and ( ∇ φ )( r ) is avector) and c J ′ are the ensuing numerical coefficients.This example illustrates the main features of thederivation, namely, (i) for the Galilean transformationsonly terms with first-order derivatives of G ( r , r ′ ) occurin the final expression for ρ Gβ , (ii) local densities appear-ing in the sum are of one order less in derivatives thanthe density being transformed, and (iii) the tensor or-der is preserved so that a local density is transformedinto a sum of densities which can couple with a vector to the same tensor order. This leads to the ansatz for theGalilean transformation, ρ Gβ ( r ) = i X γ c ( β, γ ) [ ρ γ ( r )( ∇ φ )( r )] J , (40)where c ( β, γ ) are numerical coefficients. Similarly, for thefull gauge transformation the corresponding ansatz reads ρ Gβ ( r ) = i X γmI c mI ( β, γ ) [ ρ γ ( r ) [ D mI φ ] I ( r )] J . (41)In both cases, the numerical coefficients can be foundby using the method outlined above, combined with arepeated use of the 6 j -symbols.However, instead of using this method, it turned outto be more efficient to proceed in another way. First,by using symbolic programming [38], we constructed thetransformed densities ρ Gβ ( r ) explicitly in terms of deriva-tives of the density matrices and the gauge angle. Then,from the resulting expression the ansatz (40) or (41) wassubtracted, which gave equations for the numerical coeffi-cients by requesting that these differences must be identi-cally equal to zero. Because these equations must hold forall density matrices and gauge angles, we could randomlyassume arbitrary values for these quantities and theirderivatives. In this way, all linearly independent equa-tions for the coefficients could be obtained and solved an-alytically, again by using symbolic programming. The so-lutions were then double-checked by using the full formsof the densities.
3. Galilean or gauge invariant EDF
A Galilean or gauge-invariant EDF is the one whichdoes not change upon inserting Galilean or gauge-transformed densities (37) into the energy density ofEq. (31). Since terms quadratic in G ( r , r ′ ) can bedropped, the condition for the Galilean or gauge invari-ance reads Z d r X C βmI,γ (cid:16)(cid:2) ρ Gβ ( r ) [ D mI ρ γ ( r )] J ′ (cid:3) + h ρ β ( r ) (cid:2) D mI ρ Gγ ( r ) (cid:3) J ′ i (cid:17) = 0 , (42)where C βmI,γ is a short-hand notation for the couplingconstants C n ′ L ′ v ′ J ′ mI,nLvJ , and the sum runs over all the termsin the energy density.The task now is to group together all proportionalterms in Eq. (42). In doing so, we do not aim at obtain-ing an invariant energy density but an invariant EDF andtotal energy. Therefore, after densities (40) or (41) are inserted into Eq. (42), all terms must be integrated byparts to obtain some standard form, where terms equalthrough integration by parts become identical.Finally, Eq. (42) can be transformed into a sum ofindependent terms using recoupling. In this expressioneach term is multiplied by a specific linear combinationof coupling constants C n ′ L ′ v ′ J ′ mI,nLvJ . The Galilean or gaugeinvariance of the EDF then means that these linear com-2binations must all vanish. This gives a set of linear equa-tions that must be fulfilled for an invariant EDF. On theone hand, if a given coupling constant appears in noneof these linear equations, the corresponding term of theEDF is invariant on its own, and the corresponding cou-pling constant is not restricted by the Galilean or gaugesymmetry. On the other hand, for some coupling con-stants the only solution can be the value of zero, andthen the corresponding term cannot appear in the in-variant energy density.Among all the remaining coupling constants, we mayalways select a subset of those that we will call the de-pendent ones, and express them as linear combinationsof the other ones, which we will call the independentones. This procedure is highly non-unique and can berealized in very many different ways, However, when thedependent coupling constants in function of independentones are inserted back into the energy density (31), lin-ear combinations of terms appearing at each independentcoupling constant will all be invariant with respect to theGalilean or gauge transformations.Then, the energy density of Eq. (31) takes the form H ( r ) = X n ′ L ′ v ′ J ′ mI,nLvJ,J ′ C n ′ L ′ v ′ J ′ mI,nLvJ G n ′ L ′ v ′ J ′ mI,nLvJ ( r ) , (43)where the sum runs over indices that correspond to un-restricted and independent coupling constants, which wejointly call free coupling constants. For a term in Eq. (43)that corresponds to an unrestricted coupling constant,we have G n ′ L ′ v ′ J ′ mI,nLvJ ( r ) = T n ′ L ′ v ′ J ′ mI,nLvJ ( r ), i.e., one term inthe energy density of Eq. (31) is Galilean or gauge in-variant. For a term in Eq. (43) that corresponds to anindependent coupling constant, G n ′ L ′ v ′ J ′ mI,nLvJ ( r ) is equal toa specific linear combination of terms T from the originalenergy density (31). These linear combinations are listedin Appendix C.We performed the analysis along these lines for energydensities of orders 0, 2 , 4, and 6, and the obtained resultsare listed in Appendix C. Derivations were performed byusing symbolic programming [38] and employed the tech-nique of forming linear equations by randomly assign-ing values to local densities and their derivatives, whichwe also used above. Numbers of linearly independentGalilean or gauge invariant terms are listed in Tables Vand VI, and plotted in Fig. 1.It turns out that only at orders 0 and 2, i.e., for thestandard Skyrme functional, all Galilean invariant com-binations of terms are also gauge invariant. At orders 4and 6, there are only 6 gauge invariant terms available,while the numbers of Galilean invariant terms equal to15 and 26, respectively. This is much less than the totalnumbers of terms at these orders, which are equal to 45and 129, respectively. Altogether, at N LO we obtain theEDF parametrized in general by 188 coupling constants,and by 50 or 21 coupling constants if Galilean or gaugeinvariance is assumed. If isoscalar and isovector channelsare included, all these numbers must be multiplied by a factor of two.
IV. ENERGY DENSITY AT N LO WITHCONSERVED SPHERICAL, SPACE-INVERSION,AND TIME-REVERSAL SYMMETRIES
In this section, we apply the results obtained above tothe simplest case of spherical even-even nuclei [28], whereone can assume that the spherical symmetry, along withthe space inversion and time reversal, are simultaneouslyconserved symmetries. In this case, all primary densities ρ nLvJ (23), which we listed in Tables III and IV, musthave the form [39]: ρ nLvJ ( r ) = R JJ ( r ) ρ nLvJ ( | r | ) , (44)where R JJ ( r ) = [ r [ r . . . , [ rr ] , . . . , ] J − ] J (45)is the J th-order, rank- J stretched coupled tensor builtfrom the position vector r in exactly the same way asthe derivative operators D nL of Table I are built fromthe derivative ∇ in the spherical representation (9), and ρ nLvJ ( | r | ) is a scalar function depending only on thelength | r | of the position vector r .Indeed, due to the generalized Cayley-Hamilton(GCH) theorem, a rank- J tensor function of a rank- k tensor must be a linear combination of all independentrank- J tensors built from that rank- k tensor, with scalarcoefficients. In the GCH theorem, tensors that differ byscalar factors are not independent. In our case, onlyone independent rank- J function R JJ ( r ) can be builtfrom the rank-1 tensor (position vector r ), which givesEq. (44). The spherical symmetry assumed here is essen-tial for this argument to work, because many more inde-pendent rank- J tensors can be built when other “mate-rial” tensors (like, e.g., the quadrupole deformation ten-sor) are available.The spherical form of ρ nLvJ ( r ) (44) requires that thefollowing selection rule is obeyed: P = ( − J , (46)where P = ( − n is the space-inversion parity defined inEq. (26). For the time-even densities ( T = 1), selectionrule (46) does not impose any new restriction on localdensities built from ρ ( r , r ′ ) ( v = 0), see Table III. On theother hand, for local densities built from s ( r , r ′ ) ( v = 1),see Table IV, only the densities with L = J are allowed.In Tables III and IV, all densities allowed by the con-served spherical, space-inversion, and time-reversal sym-metries are marked with bullets ( • ). One can see thatthey correspond to quantum numbers LvJ being equalto 000 or 202 [for densities built from ρ ( r , r ′ )] and 111or 313 [for densities built from s ( r , r ′ )]. Then, it is easyto select all allowed terms in the energy density—in Ta-bles VII–XVIII and XXII these are also marked with bul-lets ( • ). Numbers of such terms are listed in Table XX3 TABLE XX: Numbers of terms defined in Eq. (30) of differentorders in the EDF up to N LO, evaluated for the conservedspherical, space-inversion, and time-reversal symmetries. Thelast two columns give numbers of terms when the Galilean orgauge invariance is assumed, respectively, see Sec. III B. Totake into account both isospin channels, the numbers of termsshould be multiplied by a factor of two.order Total Galilean Gauge0 1 1 12 4 4 44 13 9 36 32 16 3N LO 50 30 11 together with those obtained by imposing, in addition,the Galilean or gauge invariance.All results for the EDF restricted by the spherical,space-inversion, and time-reversal symmetries can nowbe extracted from the general results presented in Secs. IIand III and Appendices B and C. However, in the re-mainder of this section we give an example of how theseresults can be translated into those based on the Carte-sian representations of derivative operators (18)–(22). In-deed, in this representation, all non-zero densities can bedefined as: R = ρ, (47) R = k ρ, (48) R ↔ ab = k a k b ρ, (49) R = k ρ, (50) R ↔ ab = k k a k b ρ, (51) R = k ρ, (52)and J a = ( k × s ) a , (53) J a = k ( k × s ) a , (54) J ↔ abc = k a k b ( k × s ) c + k b k c ( k × s ) a + k c k a ( k × s ) b , (55) J a = k ( k × s ) a , (56)where k = X a k a k a , (57)and the Cartesian indices are defined as a, b, c = x, y, z .To lighten the notation, in these definitions we have omit-ted the arguments of local densities, r , and limits of r ′ = r .The six local densities (47)–(52) are the Cartesian ana-logues of densities marked in Table III with bullets ( • ),and the four local densities (53)–(56) are analogues of those marked in Table IV. However, one should notethat rank-2 densities R ↔ ab and R ↔ ab are not proportionalto ρ and ρ , respectively, and the rank-3 density J ↔ abc is not proportional to ρ . This is so, becausethey are defined in terms of the derivative operators(18)–(22), where appropriate traces have not been sub-tracted out. Nevertheless, linear relations between den-sities (47)–(56) and their spherical-representation coun-terparts ρ nLvJ can easily be worked out and will not bepresented here.Note also that the scalar densities R and R can beexpressed as the corresponding sums of the rank-2 den-sities R ↔ ab and R ↔ ab , and the vector density J a as thatof J ↔ abc . However, based on the results obtained in thespherical representation, we know that they have to betreated separately to give separate terms in the energydensity.Again, based on the results obtained in the sphericalrepresentation, we can write the N LO energy densityas a sum of contributions from zero, second, fourth, andsixth orders: H = H + H + H + H , (58)where H = C R R , (59) H = C R ∆ R + C R R + C R ∇ · J , + C J , (60)Energy densities H and H correspond, of course, tothe standard Skyrme functional [24, 36] with C = C ρ , C = C ∆ ρ + C τ , C = C τ , C = C ∇ J , and C = C J . At fourth and sixth orders, these energy densitiesread H = C R ∆ R + C R ∆ R + C R R + C R R + D R X ab ∇ a ∇ b R ↔ ab + D X ab R ↔ ab R ↔ ab + C J · ∆ J + C J · J + D J · ∇ ( ∇ · J )+ C R ∆ ( ∇ · J ) + C R ( ∇ · J )+ C R ( ∇ · J ) + D X ab R ↔ ab ∇ a J b , (61)4 H = C R ∆ R + C R ∆ R + C R ∆ R + C R R + C R ∆ R + C R R + D R ∆ X ab ∇ a ∇ b R ↔ ab + D R X ab ∇ a ∇ b R ↔ ab + D R X ab ∇ a ∇ b R ↔ ab + E X ab R ↔ ab ∆ R ↔ ab + F X abc R ↔ ab ∇ a ∇ c R ↔ cb + E X ab R ↔ ab R ↔ ab + C J · ∆ J + C J · ∆ J + C J · J + C J · J + D J · ∆ ∇ ( ∇ · J ) + D J · ∇ ( ∇ · J )+ E X abc J a ∇ b ∇ c J ↔ abc + D X abc J ↔ abc J ↔ abc + C R ∆ ( ∇ · J ) + C R ∆ ( ∇ · J )+ C R ( ∇ · J ) + C R ∆ ( ∇ · J )+ C R ( ∇ · J ) + C R ( ∇ · J )+ D R X abc ∇ a ∇ b ∇ c J ↔ abc + D X abc R ↔ ab ∇ c J ↔ abc + D X ab R ↔ ab ∆ ∇ a J b + E X ab R ↔ ab ∇ a J b + D X ab R ↔ ab ∇ a J b + E X ab R ↔ ab ∇ a ∇ b ( ∇ · J ) . (62)The energy densities above are given in terms of thecoupling constants C n ′ mn , D n ′ mn , E n ′ mn , and F n ′ mn . Theindices correspond to orders of derivatives indicated inthe same way as for the spherical-representation couplingconstants C n ′ L ′ v ′ J ′ mI,nLvJ . Linear relations between both setsof coupling constants can easily be derived and are notreported here. V. CONCLUSIONS
In the present study, we constructed nuclear energydensity functionals in terms of derivatives of densitiesup to sixth order. This constitutes the next-to-next-to-next-to-leading order (N LO) expansion of the func-tional, whereby, in this scheme, the contact and standardSkyrme forces provide the zero-order (LO) and second-order (NLO) expansions, respectively. The higher-orderterms were built to provide tools for testing convergenceproperties of methods based on energy density function-als, within the spirit of effective field theories.At N LO, depending on several options of using theenergy density functionals, the numbers of free couplingconstants are as follows. If one would like to includethe density dependence of all the coupling constants (theoption, which is not advocated here), one would haveto use 658 different terms in the functional. Full func- tional with density-independent coupling constants con-tains 188 terms, while functionals restricted by Galileanand gauge symmetries contain 50 and 21 terms, respec-tively. If both isoscalar and isovector channels are in-cluded, all these numbers must be multiplied by a factorof twoAt the present stage of searching for precise,spectroscopic-quality nuclear functionals, extensions be-yond the standard Skyrme NLO form are mandatory, seethe analysis in Ref. [40]. These may include richer densitydependencies [41, 42], higher-order derivative terms. asconstructed in the present study, terms of higher powersin densities, richer forms of functional dependence be-yond simple power expansions, and possibly many othermodifications.Further studies of higher-order energy density func-tionals requires constructing appropriate codes to solveself-consistent equations. Although this is a complicatedproblem, various techniques have already been developedthat can be used here. First of all, expressions for meanfields must be derived by using the standard methodspresented, e.g., in Refs. [24, 25, 43]. Obviously, suchmean fields will involve derivative operators up to sixthorder, so the connection with the one-body Schr¨odingerequation, discussed, e.g., in Ref. [28], will be lost. Nev-ertheless, all basis-expansions methods can still be usedand their implementation will not be basically differentthan it was done up to now at NLO. Work along theselines is now in progress.This work was supported in part by the Academy ofFinland and University of Jyv¨askyl¨a within the FIDIPROprogramme, by the Polish Ministry of Science underContract No. N N202 328234, and by the UNEDF Sci-DAC Collaboration under the U.S. Department of En-ergy grant No. DE-FC02-07ER41457.
APPENDIX A: SYMMETRIES
Within the HF approximation, symmetry propertiesof the interaction carry over to symmetry properties ofthe total HF energy. It means that whenever the inter-action is invariant with respect to a symmetry group,the total HF energy is invariant with respect to trans-forming densities (or density matrices, in general) by thesame symmetry group. However, the well-known phe-nomenon of the spontaneous symmetry breaking [32] mayrender densities and energy density themselves not invari-ant with respect to the symmetry group in question. Wemust, therefore, consider energy densities for symmetry-breaking densities.This is best illustrated by the EDF derived for theSkyrme interaction [25, 28]. The standard derivationtreats the time-reversal and isospin symmetries in a dif-ferent way than space symmetries (space inversion andother point symmetries or space rotation) [24, 25, 36].Indeed, for the time reversal, the nonlocal densities are5first split into the time-even and time-odd parts as ρ ( r , r ′ ) = ρ + ( r , r ′ ) + ρ − ( r , r ′ ) , (A1) s ( r , r ′ ) = s + ( r , r ′ ) + s − ( r , r ′ ) , (A2)where ρ ± ( r , r ′ ) = (cid:2) ρ ( r , r ′ ) ± ρ T ( r , r ′ ) (cid:3) , (A3) s ± ( r , r ′ ) = (cid:2) s ( r , r ′ ) ± s T ( r , r ′ ) (cid:3) , (A4)such that ρ T ± ( r , r ′ ) = ± ρ ± ( r , r ′ ) , (A5) s T ± ( r , r ′ ) = ± s ± ( r , r ′ ) . (A6)The superscript T here means that the nonlocal densitiesare calculated for the time-reversed many-body states.Then, in the derivation of the HF energy density, onlysquares of the time-even and time-odd densities appear,because, due to the time-reversal symmetry of the in-teraction, the cross terms do not contribute. As a con-sequence, the energy density itself is time-even. Simi-larly, for the isospin symmetry, the densities are first splitinto the isoscalar and isovector parts, for which no crossterms contribute, and the obtained energy density is anisoscalar. In what follows, we call this kind of deriva-tion ’derivation after separation of symmetries’, whichimplies that the symmetry-breaking terms are absent inthe energy density.The derivation after separation of symmetries can be il-lustrated by considering the simplest term of the Skyrmeinteraction – just the contact force, V δ ( r , r ) = t δ ( r − r ) , (A7)for which the energy density reads (we neglected theisospin degree of freedom, so only one type of particles isconsidered), H δ ( r ) = t ρ ( r ) − t s ( r ) . (A8)This energy density is invariant with respect to the timereversal of local densities, ρ T ( r ) = ρ ( r ) and s T ( r ) = − s ( r ).Here, the coupling constant multiplying the time-evendensity, C ρ = t , is not independent of the couplingconstant multiplying the time-odd density, C s = − t .This fact is not related to the time-reversal symmetrybut results from the vanishing range of the contact force.Proper treatment of the finite range corrections renderthese two coupling constants independent of one another.[15, 16]. Irrespective of zero or finite range, the isovectorand isoscalar coupling constants are also independent ofone another [24, 36].For space symmetries, the standard derivation pro-ceeds in another way, namely, the energy density is deter-mined directly for the broken-symmetry HF state. Forthe space-inversion symmetry, for example, this meansthat both parity-even and parity-odd densities, ρ P = ± ( r ) = (cid:2) ρ ( r ) ± ρ P ( r ) (cid:3) = [ ρ ( r ) ± ρ ( − r )] , (A9) s P = ± ( r ) = (cid:2) s ( r ) ± s P ( r ) (cid:3) = [ s ( r ) ± s ( − r )] , (A10) appear in ρ ( r ) and s ( r ) in the energy density of Eq. (A8), ρ ( r ) = ρ P =+1 ( r ) + ρ P = − ( r ) , (A11) s ( r ) = s P =+1 ( r ) + s P = − ( r ) . (A12)The superscript P here means that the nonlocal and lo-cal densities are calculated for the space-inversed many-body states. We call this kind of derivation ’derivationbefore separation of symmetries’, which implies that thesymmetry-breaking terms are then explicitly present inthe energy density. If the symmetry is broken, which inthe case of space inversion corresponds to ρ P = − ( r ) = 0or s P = − ( r ) = 0, then the energy density is not invariantwith respect to space inversion.The total energy, i.e., the integral of the energy density(1) is, of course, invariant with respect to space inversion,because the integration then picks up only the space-inversion-invariant parts of the integrand. Therefore, theenergy densities derived before and after separation ofsymmetries are not equal, but they are equivalent. In thecase of the space-inversion symmetry, the energy density(A8) derived after separation of symmetries reads H ′ δ ( r ) = t ρ P =+1 ( r ) + t ρ P = − ( r ) − t s P =+1 ( r ) − t s P = − ( r ) . (A13)This energy density is invariant with respect to space in-version and the coupling constants multiplying densitiesof opposite parities are not independent of one another(see also discussion in Ref. [30]).The same principle applies to other broken spatial sym-metries. For example, for the broken rotational sym-metry, the density can be split into the sum of termsbelonging to different irreducible representations of therotational group, which in this case corresponds to thestandard multipole series [44], ρ ( r ) = X λµ ρ λµ ( r ) , (A14)for ρ λµ ( r ) = ρ λµ ( r ) Y ∗ λµ ( θ, φ ) , (A15)where the multipole densities ρ λµ ( r ), ρ λµ ( r ) = Z d θ d φρ ( r ) Y λµ ( θ, φ ) , (A16)depend only on the radial coordinate r . Then, the firstterm in the energy density (A8), which is derived beforeseparation of rotational symmetries, and which is not ro-tationally invariant, is equivalent to the following energydensity derived after separation of rotational symmetries, H ρ ( r ) = t X λ √ λ + 1[ ρ λ ( r ) ρ λ ( r )] , (A17)(see Ref. [45] for an example application of this series).This energy density is rotationally invariant and the cou-pling constants multiplying different multipole densitiesare again not independent of one another.6We have presented a detailed analysis of the problemto arm ourselves with proper tools for discussing con-struction of EDFs in situations where there is no un-derlying interaction. Then, the only consideration is therequirement of invariance of the total energy with respectto all symmetries usually conserved by nuclear interac-tions. We proceed with such a construction in two dif-ferent ways as described below.
1. Symmetry-invariant energy density
Based on the derivation after separation of symme-tries, which we introduced above, it is clear that we canproceed by separating densities into irreducible represen-tations of all required symmetries and then building theenergy density by taking scalar products separately ineach of the representations. Such a construction gives asymmetry-invariant energy density, H S ( r ) = H ( r ) , (A18)where H S ( r ) denotes the energy density calculated for amany-body state transformed by the symmetry operator S . This guarantees the invariance of the EDF and to-tal energy (1) with respect to all considered symmetries.Such a strategy would also allow for using arbitrary, un-related to one another coupling constants in each of theirreducible representations.However, in practical applications, such a strategy wasnever up to now fully implemented—neither at N LO, forwhich the present study is the first attempt in the liter-ature, nor at NLO, which corresponds to the standardSkyrme functionals (see Sec. III A below). Only the timereversal and isospin symmetries were up to now treatedin this way, and below we are going to follow the samepath.For the time reversal, all local densities discussed inSec. II C are either time-even or time-odd. Indeed, thissimply follows from the facts [25] that ρ T ( r , r ′ ) = ρ ∗ ( r , r ′ ) = ρ ( r ′ , r ) , s T ( r , r ′ ) = − s ∗ ( r , r ′ ) = − s ( r ′ , r ) , (A19)which give the time-even and time-odd parts in Eqs. (A3)and (A4) as ρ + ( r , r ′ ) = ρ ∗ + ( r , r ′ ) = ρ + ( r ′ , r ) ,ρ − ( r , r ′ ) = − ρ ∗− ( r , r ′ ) = − ρ − ( r ′ , r ) , s + ( r , r ′ ) = − s ∗ + ( r , r ′ ) = − s + ( r ′ , r ) , s − ( r , r ′ ) = s ∗− ( r , r ′ ) = s − ( r ′ , r ) , (A20)i.e., ρ + ( r , r ′ ) and s − ( r , r ′ ) are real symmetric functionsand ρ − ( r , r ′ ) and s + ( r , r ′ ) are imaginary antisymmetricfunctions. Moreover, the relative momentum operator k (6), which defines derivative operators K nL , is imaginaryand antisymmetric with respect to exchanging variables r and r ′ . Altogether, it is easy to see that T -parities of primary densities ρ nLvJ ( r ) (23) are equal to ( − n + v , seeEq. (25) and columns denoted by T in Tables III and IV.Similarly, T -parities of secondary densities ρ mI,nLvJ,Q ( r )(24) are also equal to ( − n + v . Construction of the T -invariant energy density (A18) can now be realized bymultiplying densities that have identical T -parities.
2. Symmetry-covariant energy density
Treatment of space symmetries in the construction ofEDFs is another kettle of fish. Here, we base our consid-erations on the derivation before separation of symme-tries, which we introduced above, and on the fact thatinvariance of the energy density itself is not a prerequi-site for the invariance of the EDF. In fact, the EDF andtotal energy (1) are invariant with respect to symmetry S also when the energy density is covariant with S , i.e., H S ( r ) = H ( S + r S ) , (A21)where S + r S denotes the space point transformed by sym-metry S (see also discussion in Ref. [32]). Indeed, due tothe fact that space integrals are invariant, we have Z d r H ( S + r S ) = Z d r H ( r ) , (A22)which guarantees invariance of the EDF and total energy.For the space-inversion symmetry, we have H P ( r ) ≡ H (cid:2) ρ P ( r σ, r ′ σ ′ ) (cid:3) = H [ ρ ( − r σ, − r ′ σ ′ )] (A23) H ( P + r P ) ≡ H ( − r ) , (A24)and the covariance condition (A21) reads H [ ρ ( − r σ, − r ′ σ ′ )] = H ( − r ) . (A25)It is now essential to realize that the arguments of thedensity matrix, ρ ( − r σ, − r ′ σ ′ ) on which the energy den-sities in Eq. (A25) depends, are the same on both sidesof Eq. (A25). The covariance condition then tests onlythe parity of all other operators that may appear in thedefinition of local densities. Therefore, to each primarydensity ρ nLvJ ( r ) (23) we may attribute P -parity corre-sponding to the P -parity of the operator K nL only, whichis equal to ( − n , see Eq. (26) and columns denoted by P in Tables III and IV. This attribution is performed re-gardless of space-inversion properties of the nonlocal den-sities, i.e., regardless of whether the parity of the many-body state is conserved or broken. Similarly, P -paritiesof secondary densities ρ mI,nLvJ,Q ( r ) (24) are equal to( − n + m . Construction of the P -covariant energy den-sity (A21) can now be realized by multiplying densitiesthat have identical P -parities.Construction of a rotationally covariant energy densitycan be performed in an entirely analogous way. We mustonly ensure, that all tensor operators used in constructing7all terms of the energy density are always coupled to totalangular momentum (rank) zero. This coupling proceeds regardless of any transformation properties of nonlocaldensities with respect to rotation, because again, theirrotated space arguments appear on both sides of the co-variance condition (A21).It is obvious that this is the correct procedure to followwhen the rotational symmetry is not broken, and nonlo-cal densities ρ ( r , r ′ ) and s ( r , r ′ ) are scalar and vectorfunctions of their arguments, respectively. In fact, thisis how we refer to these densities throughout the entirepaper, seemingly forgetting that the rotational symme-try can be broken, and that these functions can thenhave no good tensor properties with respect to rotation.Nevertheless, in view of the covariance condition (A21), these rotational properties of broken-symmetry nonlocaldensities are irrelevant for the construction of the energydensity. APPENDIX B: PHASE CONVENTIONS
In the present study, we use four elementary buildingblocks to construct the EDF, namely, the scalar and vec-tor nonlocal densities, ρ ( r , r ′ ) and s ( r , r ′ ), along withthe total derivative ∇ and relative momentum k (6).Spherical representations of the building blocks can bedefined by using standard convention of spherical tensors[26] as ρ ( r , r ′ ) = p ρ ρ ( r , r ′ ) , (B1) s ,µ = {− , , } ( r , r ′ ) = p s n √ ( s x ( r , r ′ ) − is y ( r , r ′ )) , s z ( r , r ′ ) , − √ ( s x ( r , r ′ ) + is y ( r , r ′ )) o , (B2) ∇ ,µ = {− , , } = p ∇ n √ ( ∇ x − i ∇ y ) , ∇ z , − √ ( ∇ x + i ∇ y ) o , (B3) k ,µ = {− , , } = p k n √ ( k x − ik y ) , k z , − √ ( k x + ik y ) o , (B4)where p ρ , p s , p ∇ , and p k , are arbitrary phase factors, | p ρ | = | p s | = | p ∇ | = | p k | = 1. These phase factors definethe phase convention of the building blocks, and can beused to achieve specific phase properties of densities andterms in the EDF, as discussed in this Appendix.In order to motivate the best suitable choice of thephase convention, in Tables XXI and XXII we presentrelations between the spherical and Cartesian represen-tations of densities and terms in the EDF, respectively.All NLO densities in the Cartesian representation, whichare listed in Table XXI, are real. It is then clear that thephase convention, which would render all NLO densitiesin the spherical representation real does not exist. How-ever, for the phase factors p ρ , p s , p ∇ , and p k equal to ± ± i , in the spherical representation all NLO densitiesand terms in the EDF are either real or imaginary.Among many options of choosing the phase convention,in the present study we set p ρ = +1 , p s = − i, p ∇ = − i, and p k = − i. (B5)This choice is unique in the fact that all scalar densi-ties and all terms in the EDF are then characterized byphase factors +1 connecting the spherical and Cartesianrepresentations, see the last columns in Tables XXI andXXII. This allows for the closest possible relationshipsbetween both representations, which may facilitate theuse of the spherical representation as it is introduced inthe present study. In particular, relations between cou-pling constants up to NLO (Table XXII) and standard coupling constants in the Cartesian representation [24]then read, for terms depending on time-even densities: C , = C ρ , (B6) C , = √ (cid:0) C ∆ ρ + C τ (cid:1) , (B7) C , = √ C τ , (B8) C , = 3 C J , (B9) C , = √ C J , (B10) C , = √ C J , (B11) C , = √ C ∇ J , (B12)and for terms depending on time-odd densities: C , = √ C s , (B13) C , = √ C j , (B14) C , = 3 (cid:0) C ∆ s + C T (cid:1) + C F − C ∇ s , (B15) C , = √ (cid:0) C F − C ∇ s (cid:1) , (B16) C , = 3 C T + C F , (B17) C , = √ C F , (B18) C , = √ C ∇ j . (B19)At the same time, all vector densities in Table XXI andvector operators in Eqs. (B2)–(B4) are consistently char-acterized by phase factors − i connecting the sphericaland Cartesian representations.8 TABLE XXI: Spherical and Cartesian representations of local densities (24) up to NLO. Only scalar densities and the µ = 0components of vector densities are shown. Numbers in the first column refer to numbers of primary densities (23) shown inTables III and IV. The last column shows factors preceding densities in the Cartesian representation evaluated for the phaseconventions of Eq. (B5). Time-even densities are marked by using the bold-face font.No. ρ mI,nLvJ,Qµ Cartesian representation [24, 25, 30] Phase1 ρ , , = [ ρ ] = p ρ ρ +1 ρ , , = [ ∇ ρ ] = p ∇ p ρ ∇ z ρ − i ρ , , = [[ ∇∇ ] ρ ] = − p ∇ p ρ √ ∆ ρ +12 ρ , , = [ kρ ] = p k p ρ j z − iρ , , = [ ∇ [ kρ ] ] = − p ∇ p k p ρ √ ∇ · j +1 ρ , , = [ ∇ [ kρ ] ] = ip ∇ p k p ρ √ ( ∇ × j ) z − i ρ , , = [[ kk ] ρ ] = − p k p ρ √ ` τ − ∆ ρ ´ +117 ρ , , = [ s ] = p s s z − iρ , , = [ ∇ s ] = − p ∇ p s √ ∇ · s +1 ρ , , = [ ∇ s ] = ip ∇ p s √ ( ∇ × s ) z − iρ , , = [[ ∇∇ ] s ] = − p ∇ p s √ ∆ s z − iρ , , = [[ ∇∇ ] s ] = − p ∇ p s √ (3 ∇ z ∇ · s − ∆ s z ) − i ρ , , = [ ks ] = − p k p s √ J (0) +1 ρ , , = [ ∇ [ ks ] ] = − p ∇ p k p s √ ∇ z J (0) − i ρ , , = [ ks ] = ip k p s √ J z − i ρ , , = [ ∇ [ ks ] ] = − ip ∇ p k p s √ ∇ · J +1 ρ , , = [ ∇ [ ks ] ] = − p ∇ p k p s ( ∇ × J ) z − i ρ , , = [[ kk ] s ] = − p k p s √ ` T z − ∆ s z ´ − i ρ , , = [[ kk ] s ] = − p k p s √ ` F z − ∇ z ∇ · s − T z + ∆ s z ´ − i TABLE XXII: Spherical and Cartesian representations of terms in the EDF (30) up to NLO. The last column shows factorspreceding terms in the Cartesian representation evaluated using the phase conventions of Eq. (B5). Integration by partswas used to transform s · ∇∇ · s into − ( ∇ · s ) , which is the term used previously in Refs. [24, 30]. Coupling constantscorresponding to terms that depend on time-even densities are marked by using the bold-face font. Bullets ( • ) mark couplingconstants corresponding to terms that do not vanish for conserved spherical, space-inversion, and time-reversal symmetries, seeSec. IV. No. C n ′ L ′ v ′ J ′ mI,nLvJ [ ρ n ′ L ′ v ′ J ′ [ D mI ρ nLvJ ] J ′ ] Cartesian representation Phase1 • C , [ ρρ ] = p ρ ρ +12 C , [ ss ] = − p s √ s +13 • C , [ ρ [[ ∇∇ ] ρ ] ] = − p ∇ p ρ √ ρ ∆ ρ +14 • C , [ ρ [[ kk ] ρ ] ] = − p k p ρ √ ` ρτ − ρ ∆ ρ ´ +15 C , [[ ks ] [ ks ] ] = p k p s “ J (0) ” +16 • C , [[ ks ] [ ks ] ] = p k p s √ J +17 C , [[ ks ] [ ks ] ] = p k p s √ P ab J (2) ab J (2) ab +18 • C , [ ρ [ ∇ [ ks ] ] ] = − ip ∇ p k p s p ρ √ ρ ∇ · J +19 C , [[ kρ ] [ kρ ] ] = − p k p ρ √ j +110 C , [ s [[ ∇∇ ] s ] ] = p ∇ p s s ∆ s +111 C , [ s [[ ∇∇ ] s ] ] = p ∇ p s − √ ` ( ∇ · s ) + s ∆ s ´ +112 C , [ s [[ kk ] s ] ] = p k p s ` s · T − s ∆ s ´ +113 C , [ s [[ kk ] s ] ] = p k p s √ “ s · F + ( ∇ · s ) − s · T + s ∆ s ” +114 C , [[ kρ ] [ ∇ s ] ] = − ip ∇ p k p s p ρ √ j · ∇ × s +1 A ∗ λµ = P A ( − λ − µ A λ, − µ , (B20)where P A = ±
1. For nonlocal densities (B1)–(B2),Eq. (B20) holds separately for their time-even and time-odd parts, split as in Eqs. (A1) and (A2). UsingEqs. (A20) we then have P ρ + = + p ρ , P ρ − = − p ρ ,P s + = + p s , P s − = − p s ,P ∇ = − p ∇ , P k = p k , (B21)which for the phase convention of Eq. (B5) reads P ρ + = +1 , P ρ − = − ,P s + = − , P s − = +1 ,P ∇ = +1 , P k = − . (B22)Standard rule (B20) propagates through the angularmomentum coupling, i.e., if signs P A and P A ′ character-ize tensors A λ and A ′ λ ′ , respectively, then the coupledtensor, A ′′ λ ′′ µ ′′ = [ A λ A ′ λ ′ ] λ ′′ µ ′′ = X µµ ′ C λ ′′ µ ′′ λµλ ′ µ ′ A λµ A ′ λ ′ µ ′ , (B23)is characterized by the product of signs P A ′′ = P A P A ′ .Therefore, coupled higher-order densities (24) are char-acterized by signs, P ρ mI,nLvJ,Q = P m ∇ P nk P vT , (B24)where v = 0 or 1 denotes the scalar or vector density, ρ or s , respectively, and T = +1 or T = − k derivative determinethe time-reversal symmetry of each local density, so that T = ( − n + v . From Eqs. (B21) we then obtain P ρ mI,nLvJ,Q = ( − m + n + v p m ∇ p nk p v . (B25)which for the phase convention of Eq. (B5) reads P ρ mI,nLvJ,Q = +1 (B26)for all densities. Therefore, the phase convention ofEq. (B5) ensures that scalar densities and all terms inthe EDF are always real. APPENDIX C: RESULTS FOR THE GALILEANOR GAUGE INVARIANT ENERGY DENSITYFUNCTIONAL
As discussed in Sec. III B 3, when the Galilean or gaugeinvariance is imposed on the EDF, this induces specific
TABLE XXIII: Numbers of unrestricted, vanishing, indepen-dent, and dependent coupling constants in the EDF at zero,second, fourth, and sixth orders. Left and right columns cor-respond to the Galilean and gauge symmetries imposed, re-spectively. Galilean Gaugeorder 0 2 4 6 0 2 4 6unrestricted 2 3 3 3 2 3 3 3vanishing 0 0 0 0 0 0 27 100independent 0 4 12 23 0 4 3 3dependent 0 5 30 103 0 5 12 23 constraints on the coupling constants and terms of thefunctional. We pointed out that there can be three dis-connected classes of terms in the EDF with related prop-erties of the coupling constants:1. Terms that are invariant with respect to theGalilean or gauge transformation, and, therefore,the corresponding coupling constants are not re-stricted by the imposed symmetries.2. Terms that cannot appear in the energy densitywhen the Galilean or gauge symmetry is imposed,and, therefore, the corresponding coupling con-stants must be equal to zero.3. Terms that can appear in the energy density onlyin certain specific linear combinations with otherterms. This means that the coupling constants cor-responding to these terms must obey specific linearconditions. We then distinguish:(a) independent coupling constants, which mul-tiply invariant combinations of terms and,therefore, their values are not restricted by theimposed symmetries and(b) dependent coupling constants, which are equalto specific linear combinations of independentcoupling constants, and, therefore, their val-ues are in this way restricted by the imposedsymmetries.Division into the sets of independent and dependentcoupling constants is not unique, and below, in eachcase, we present only one specific choice thereof.In Table XXIII we show numbers of unrestricted, vanish-ing, independent, and dependent coupling constants thatappear at a given order when either Galilean or gaugesymmetry is imposed. In what follows, we use the nameof a free coupling constant to denote either the unre-stricted or independent one. Indeed, in the Galilean orgauge invariant energy density (43), these two groups ofcoupling constants become free parameters.0Below we simultaneously discuss the Galilean andgauge symmetries. In doing so, we use the fact that theGalilean symmetry is a special case of the gauge sym-metry, and, therefore, the latter may impose more re-strictions on the EDF than the former. At NLO, this isnot the case, and the Galilean and gauge symmetries im-pose, in fact, identical restrictions on the EDF [25, 36].However, at higher orders, restrictions imposed by theGalilean and gauge symmetries are very different.Both zero-order terms in the EDF, which correspond tothe contact interaction, are Galilean and gauge invariant,i.e., these symmetries do not restrict the form of the EDFat LO. In the three following sections we give results forsecond, fourth, and sixth orders, respectively.
1. Second order
At second order, we obtain the same restrictions ofthe EDF as those already identified for the Skyrme func-tional, see Ref. [24] for a complete list thereof. Then, 5dependent coupling constants are equal to specific linearcombinations of 4 independent ones: C , = − C , , (C1) C , = − C , − √ C , , (C2) C , = q C , − √ C , , (C3) C , = − √ C , − C , , (C4) C , = C , . (C5)These relations are obtained by imposing either Galileanor gauge invariance. In Eqs. (C1)–(C5), coupling con-stants corresponding to terms that depend on time-evendensities are marked by using the bold-face font. Thesame convention also applies below.At this order, the Galilean or gauge invariant energydensity of Eq. (43) is composed of three terms corre-sponding to unrestricted coupling constants: G , = T , , (C6) G , = T , , (C7) G , = T , , (C8)and of four terms corresponding to the independent cou-pling constants: G , = T , − T , , (C9) G , = T , + T , , (C10) G , = − T , − √ T , − √ T , + T , , (C11) G , = − √ T , + q T , − T , + T , . (C12) Again, terms that depend on time-even densities aremarked by using the bold-face font. Altogether, 7 freecoupling constants (3 unrestricted and 4 independent)define the Galilean or gauge invariant EDF at second or-der, cf. Table VI.
2. Fourth order
At fourth order, imposing either Galilean or gaugesymmetry forces 12 dependent coupling constants to bespecific linear combinations of 3 independent ones: C , = √ C , , (C13) C , = √ C , , (C14) C , = − √ C , , (C15) C , = − q C , − √ C , , (C16) C , = − √ C , , (C17) C , = − √ C , − C , , (C18) C , = − q C , , (C19) C , = q C , + √ C , , (C20) C , = √ C , + √ C , , (C21) C , = C , , (C22) C , = q C , + √ C , , (C23) C , = q C , + q C , . (C24)At this order, there are 3 stand-alone Galilean and gaugeinvariant terms: G , = T , , (C25) G , = T , , (C26) G , = T , , (C27)and 3 Galilean and gauge invariant linear combinationsof terms, corresponding to the 3 independent couplingconstants:1 G , = √ T , + T , − √ T , + √ T , , (C28) G , = √ T , − √ T , − T , − q T , + √ T , + T , + T , + √ T , + q T , , (C29) G , = √ T , + q T , − q T , − √ T , − √ T , + q T , + T , + q T , , (C30)Altogether, these 6 free coupling constants (3 unre-stricted and 3 independent) occur in both the Galileanand gauge invariant EDF at fourth order, cf. Table VI.Apart from these 6 free and 12 dependent couplingconstants, the gauge invariance requires that all the re-maining 27 coupling constants are equal to zero. These27 constants are allowed to be non-zero if the Galileansymmetry is imposed instead of the full gauge invariance.Then, there are 18 dependent coupling constants that areforced to be linear combinations of 9 independent ones: C , = − C , , (C31) C , = − C , , (C32) C , = − C , − √ C , , (C33) C , = − C , − √ C , − q C , , (C34) C , = q C , − √ C , , (C35) C , = − √ C , − q C , − q C , , (C36) C , = − √ C , − C , , (C37) C , = C , − √ C , + q C , , (C38) C , = √ C , + q C , + √ C , , (C39) C , = √ C , + 3 q C , +2 q C , , (C40) C , = C , , (C41) C , = − q C , , (C42) C , = − q C , , (C43) C , = √ C , , (C44) C , = C , , (C45) C , = − q C , , (C46) C , = √ C , , (C47) C , = − q C , . (C48)Finally, we list 9 combinations of terms that are in-variant with respect to the Galilean symmetry and cor-respond to the independent coupling constants:2 G , = − √ T , − T , + T , + T , + √ T , + √ T , , (C49) G , = − q T , − √ T , + T , − √ T , + q T , + 3 q T , , (C50) G , = T , − T , , (C51) G , = T , − T , , (C52) G , = T , + T , , (C53) G , = − q T , + √ T , + T , − q T , + √ T , − q T , + T , − q T , , (C54) G , = − T , − √ T , − √ T , + T , , (C55) G , = − √ T , + q T , − T , + T , , (C56) G , = − q T , − q T , + q T , + √ T , +2 q T , + T , . (C57)
3. Sixth order
An entirely analogous pattern of terms and cou-pling constants appears at sixth order. Imposing eitherGalilean or gauge symmetry forces 23 dependent cou-pling constants to be specific linear combinations of 3independent ones:. C , = − q C , , (C58) C , = − √ C , , (C59) C , = − q C , , (C60) C , = q C , , (C61) C , = √ C , , (C62) C , = − q C , − √ C , , (C63) C , = − √ C , , (C64) C , = − √ C , − C , , (C65) C , = − q C , , (C66) C , = − q C , − √ C , , (C67) C , = − √ C , , (C68) C , = − √ C , − C , , (C69) C , = − q C , , (C70) C , = − √ C , − C , , (C71) C , = − q C , , (C72) C , = − √ C , − √ C , , (C73) C , = q C , + √ C , , (C74) C , = q C , + √ C , , (C75) C , = 6 C , , (C76) C , = C , , (C77) C , = q C , + √ C , , (C78) C , = q C , + 18 q C , , (C79) C , = √ C , . (C80)At this order, there are 3 stand-alone Galilean and gaugeinvariant terms: G , = T , , (C81) G , = T , , (C82) G , = T , , (C83)and 3 Galilean and gauge invariant linear combinationsof terms, corresponding to the 3 independent couplingconstants:3 G , = √ T , + T , − √ T , − q T , + q T , − q T , , (C84) G , = − √ T , − T , − q T , − √ T , − T , − q T , − T , − √ T , + √ T , + 6 T , + √ T , + T , + T , + √ T , + 18 q T , + √ T , , (C85) G , = − q T , − √ T , − √ T , − q T , − √ T , − √ T , − √ T , − q T , − √ T , + q T , + q T , + q T , + T , + q T , . (C86)Altogether, 6 free coupling constants (3 unrestricted and3 independent) occur in both the Galilean and gauge in-variant EDF at sixth order, cf. Table VI.Apart from the 6 free and 23 dependent coupling con-stants, at sixth order the gauge invariance requires thatall the remaining 100 coupling constants are equal tozero. These 100 constants are allowed to be non-zeroif the Galilean symmetry is imposed instead of the fullgauge invariance. Then, there are 80 dependent couplingconstants that are forced to be linear combinations of 20independent ones: C , = − C , , (C87) C , = − C , , (C88) C , = √ C , , (C89) C , = − q C , , (C90) C , = √ C , , (C91) C , = − q C , , (C92) C , = − q C , , (C93) C , = − √ C , , (C94) C , = q C , , (C95) C , = √ C , − √ C , , (C96) C , = − C , − √ C , , (C97) C , = √ C , + q C , + q C , − √ C , , (C98) C , = q C , + q C , , (C99) C , = 2 √ C , + 2 √ C , , (C100) C , = q C , − √ C , , (C101) C , = √ C , − C , , (C102) C , = √ C , − √ C , , (C103) C , = √ C , − q C , , (C104) C , = − C , , (C105) C , = C , , (C106) C , = − C , , (C107) C , = − √ C , − C , , (C108) C , = q C , + √ C , − √ C , , (C109) C , = − q C , , (C110) C , = √ C , + q C , + q C , − q C , , (C111) C , = √ C , − √ C , , (C112) C , = √ C , + √ C , + √ C , − √ C , , (C113) C , = − √ C , − √ C , , (C114)4 C , = √ C , − C , + √ C , , (C115) C , = − √ C , − √ C , − q C , , (C116) C , = − √ C , − q C , − q C , , (C117) C , = C , , (C118) C , = − √ C , − C , + C , , (C119) C , = C , , (C120) C , = q C , − √ C , − √ C , , (C121) C , = C , , (C122) C , = − √ C , − q C , − q C , , (C123) C , = − q C , − √ C , − √ C , + q C , , (C124) C , = − √ C , − C , + C , , (C125) C , = √ C , − √ C , − √ C , , (C126) C , = − √ C , − q C , − q C , , (C127) C , = C , + √ C , , (C128) C , = − √ C , − C , − C , + q C , , (C129) C , = √ C , + √ C , − √ C , , (C130) C , = − √ C , − q C , + √ C , + q C , , (C131) C , = √ C , − √ C , , (C132) C , = √ C , − q C , − √ C , , (C133) C , = − q C , − q C , + C , , (C134) C , = q C , + √ C , , (C135) C , = − q C , − √ C , + q C , + q C , , (C136) C , = C , , (C137) C , = − q C , , (C138) C , = C , , (C139) C , = − √ C , , (C140) C , = − q C , , (C141) C , = − q C , , (C142) C , = √ C , , (C143) C , = √ C , , (C144) C , = C , , (C145) C , = −√ C , , (C146)5 C , = √ C , , (C147) C , = C , , (C148) C , = − q C , , (C149) C , = − q C , , (C150) C , = − q C , , (C151) C , = √ C , , (C152) C , = C , , (C153) C , = − q C , , (C154) C , = √ C , , (C155) C , = − √ C , , (C156) C , = − q C , , (C157) C , = √ C , , (C158) C , = C , , (C159) C , = − q C , , (C160) C , = − q C , , (C161) C , = √ C , , (C162) C , = C , , (C163) C , = − q C , , (C164) C , = − q C , , (C165) C , = − √ C , . (C166) Finally, we list 20 combinations of terms that are in-variant with respect to the Galilean symmetry and cor-respond to the independent coupling constants:6 G , = √ T , + T , − √ T , + √ T , , (C167) G , = T , + √ T , − T , − √ T , − q T , , (C168) G , = T , + T , + √ T , + q T , + q T , + √ T , − T , − √ T , − T , − √ T , − q T , − √ T , − T , − q T , − q T , − √ T , , (C169) G , = T , + 2 √ T , − q T , + √ T , − √ T , − √ T , , (C170) G , = √ T , − T , + T , − √ T , + q T , + q T , + T , + √ T , + T , − √ T , − q T , − √ T , − √ T , − T , + √ T , − q T , − q T , + √ T , + q T , − q T , , (C171) G , = − √ T , − √ T , − q T , + T , + √ T , − √ T , + q T , − √ T , − √ T , + q T , + T , + q T , + q T , , (C172) G , = √ T , + q T , + T , − √ T , , (C173) G , = T , − T , , (C174) G , = T , − T , , (C175) G , = − q T , − q T , + q T , + T , − q T , , (C176) G , = T , + T , , (C177) G , = − q T , + √ T , + T , − q T , + √ T , − q T , + T , − q T , , (C178) G , = T , + √ T , − √ T , − √ T , + T , + T , , (C179)7 G , = − q T , + √ T , + T , − q T , − q T , − √ T , − q T , + √ T , + T , − q T , + √ T , + T , − √ T , − q T , + √ T , + T , − q T , − q T , , (C180) G , = − √ T , + q T , − √ T , + T , , (C181) G , = 2 √ T , + √ T , − q T , + T , − q T , , (C182) G , = T , − T , , (C183) G , = T , + T , − T , − √ T , − √ T , + T , + √ T , + q T , , (C184) G , = T , − √ T , + q T , − T , − q T , + T , + √ T , − √ T , + √ T , , (C185) G , = − T , + q T , + √ T , + √ T , + √ T , − √ T , + √ T , − √ T , + q T , + T , − √ T , − q T , − √ T , + √ T , − √ T , + T , + √ T , − √ T , + √ T , − √ T , − q T , − q T , . (C186) [1] T.H.R. Skyrme, Phil. Mag. , 1043 (1956).[2] T.H.R. Skyrme, Nucl. Phys. , 615 (1959).[3] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864(1964).[4] W. Kohn and L.J. Sham, Phys. Rev. , A1133 (1965).[5] R.M. Dreizler and E.K.U. Gross, Density FunctionalTheory (Springer, Berlin, 1990).[6] I.Zh. Petkov and M.V. Stoitsov,
Nuclear Density Func-tional Theory , Oxford Studies in Physics Vol. 14, (Claren-don Press, Oxford, 1991).[7] G. Giuliani and G. Vignale,
Quantum theory of theelectron liquid (Cambridge University Press, Cambridge,2005).[8] E. Castellani, Stud. Hist. Phil. Mod. Phys. , 251(2002); arXiv:physics/0101039.[9] S. Weinberg, The Quantum Theory of Fields (CambridgeUniversity Press, Cambridge, 1996-2000) Vols. I-III.[10] J.W. Negele, Phys. Rev. C , 1260 (1970).[11] J.W. Negele and D. Vautherin, Phys. Rev. C , 1472(1972).[12] J.W. Negele, in Effective Interactions and Operators inNuclei , Lecture Notes in Physics 40 (Springer, Berlin,1975), p. 270.[13] X. Campi and A. Bouyssy, Phys. Lett. , 263 (1978).[14] F. Hofmann and H. Lenske, Phys. Rev. C , 2281 (1998).[15] J. Dobaczewski, in Trends in Field Theory Research , ed.O. Kovras (Nova Science Publishers, New York, 2005) p.157, nucl-th/0301069.[16] J. Dobaczewski, unpublished.[17] S. N. Maximoff and G. E. Scuseria, Jour. Chem. Phys. , 10591 (2001).[18] D.R. Entem and R. Machleidt, Phys. Rev. C , 041001(2003).[19] W. C. Haxton Phys. Rev. C , 034005 (2008).[20] N. Kaiser, S. Fritsch, and W. Weise, Nucl. Phys. A724 ,47 (2003).[21] S.J. Puglia, A. Bhattacharyya, and R.J. Furnstahl, Nucl.Phys. A , 145 (2003).[22] P. Finelli, N. Kaiser, D. Vretenar, and W. Weise, Nucl.Phys.
A770 , 1 (2006).[23] G.P. Lepage, nucl-th/9706029.[24] E. Perli´nska, S.G. Rohozi´nski, J. Dobaczewski, and W.Nazarewicz, Phys. Rev. C , 014316 (2004).[25] Y.M. Engel, D.M. Brink, K. Goeke, S.J. Krieger, and D.Vautherin, Nucl. Phys. A , 215 (1975).[26] D.A. Varshalovich, A.N. Moskalev, and V.K. Kherson-skii, Quantum Theory of Angular Momentum (WorldScientific, Singapore 1988).[27] J. Applequist, J. Phys. A , 4303 (1989). [28] D. Vautherin and D.M. Brink, Phys. Rev. C , 626(1972).[29] M. Beiner, H. Flocard, N. Van Giai, and P. Quentin,Nucl. Phys. A , 29 (1975).[30] J. Dobaczewski and J. Dudek, Acta Phys. Pol. B27 , 45(1996).[31] S.G. Nilsson and I. Ragnarsson,
Shapes and shells in nu-clear structure (Cambridge University Press, Cambridge,1995).[32] P. Ring and P. Schuck,
The Nuclear Many-Body Problem (Springer-Verlag, Berlin, 1980).[33] R. Machleidt and I. Slaus, J. Phys. G , R69 (2001).[34] J.S. Bell and T.H.R Skyrme. Phil. Mag. , 1055 (1956).[35] H. Feldmeier, T. Neff, R. Roth, and J. Schnack, Nucl.Phys. A , 61 (1998).[36] J. Dobaczewski and J. Dudek, Phys. Rev. C52 , 1827(1995);
C55 , 3177(E) (1997).[37] M. Bender, J. Dobaczewski, J. Engel, and W. Nazarewicz, Phys. Rev.
C65 , 054322 (2002).[38] S. Wolfram,
The Mathematica book / Stephen Wolfram ,4th ed. (Cambridge University Press, Cambridge, 1999).[39] S.G. Rohozi´nski, J. Dobaczewski and W. Nazarewicz, un-published.[40] M. Kortelainen, J. Dobaczewski, K. Mizuyama, and J.Toivanen, Phys. Rev. C , 064307 (2008).[41] B. Cochet, K. Bennaceur, P. Bonche, T. Duguet, and J.Meyer, Nucl. Phys. A , 34 (2004).[42] M. Kortelainen et al. , unpublished.[43] Ph. Chomaz, K.H.O. Hasnaoui, F. Gulminelli,arXiv:nucl-th/0610027.[44] J.D. Jackson, Classical Electrodynamics (Wiley, NewYork, 1975).[45] T.R. Werner, J. Dobaczewski, M.W. Guidry, W.Nazarewicz, and J.A. Sheikh, Nucl. Phys.