1 INTRODUCTION
###### Abstract

Random Phase Approximation (RPA) is the basic method for calculation of excited states of nuclei over the Hartree-Fock ground state, suitable also for energy density functionals (EDF or DFT). We developed a convenient formalism for expressing densities and currents in a form of reduced matrix elements, which allows fast calculation of spectra for spherical nuclei. All terms of Skyrme functional were taken into account, so it is possible to calculate electric, magnetic and vortical/toroidal/compression transitions and strength functions of any multipolarity. Time-odd (spin) terms in Skyrme functional become important for magnetic and isovector toroidal (i.e. second-order term in long-wave expansion of E1) transitions. It was also found that transition currents in pygmy region (low-lying part of E1 resonance) exhibit isoscalar toroidal flow, so the previously assumed picture of neutron-skin vibration is not the only mechanism present in pygmy transitions.

RPA calculations with heavy axially-symmetric nuclei now become feasible on ordinary PC. Detailed formulae for axial Skyrme RPA are given. Some numerical results are shown in comparison with the approximate approach of separable RPA, previously developed in our group for fast calculation of strength functions.

Skyrme RPA for spherical and axially symmetric nuclei

Anton Repko, Jan Kvasil, V.O. Nesterenko and P.-G. Reinhard

Institute of Particle and Nuclear Physics, Charles University, V Holešovičkách 2, 180 00 Prague 8, Czech Republic

Joint Institute for Nuclear Research, Dubna, Moscow region 141980, Russia

Institut für Theoretische Physik II, Universität Erlangen, D-91058 Erlangen, Germany

## 1 Introduction

Energy density functional in nuclear physics is a self-consistent microscopic phenomenological approach to calculate nuclear properties and structure over the whole periodic table. The method is analogous to Kohn-Sham density functional theory (DFT) used in electronic systems. Three types of functionals are frequently used nowadays: non-relativistic Skyrme functional [1] with zero-range two-body and density dependent interaction, finite-range Gogny force [2] and relativistic (covariant) mean-field [3]. Typical approach employs Hartree-Fock-Bogoliubov or HF+BCS calculation scheme to obtain ground state and single-(quasi)particle wavefunctions and energies. These results are then utilized to fit the parameters of the functionals, thus obtaining various parametrizations suitable for specific aims, such as: calculation of mass-table, charge radii, fission barriers, spin-orbit splitting and giant resonances.

Random Phase Approximation (RPA) is a textbook standard [4] to calculate one-phonon excitations of the nucleus. Increasing computing power has enabled to employ fully self-consistent residual interaction derived from the same density functional as the underlying ground state. While the spherical nuclei can be treated directly (by matrix diagonalization) [5, 6, 7], axially deformed nuclei still pose certain difficulties due to large matrix dimensions [8]. Our group developed a separable RPA (SRPA) approach [9, 10], which greatly reduces the computational cost by utilizing separable residual interaction, entirely derived from the underlying functional by means of multi-dimensional linear response theory.

Present article gives a convenient formalism for rotationally-invariant treatment of spherical full RPA, and also gives detailed expressions for matrix elements in axial symmetry. Both time-even and time-odd terms of Skyrme functional are employed, so the method is suitable for various electric and magnetic multipolarities. Finally, some results are shown by means of strength functions and transition currents.

## 2 Skyrme functional and density operators

Skyrme interaction is defined as

 ^VSk(→r1,→r2) =t0(1+x0^Pσ)δ(→r1−→r2)−18t1(1+x1^Pσ)[(←∇1−←∇2)2δ(→r1−→r2)+δ(→r1−→r2)(→∇1−→∇2)2] +14t2(1+x2^Pσ)(←∇1−←∇2)⋅δ(→r1−→r2)(→∇1−→∇2)+16t3(1+x3^Pσ)δ(→r1−→r2)ρα(→r1+→r22) +i4t4(→σ1+→σ2)⋅[(←∇1−←∇2)×δ(→r1−→r2)(→∇1−→∇2)] (1)

with parameters and a spin-exchange operator . Since it is a zero-range interaction, the solution of many-body problem by Hartree-Fock method can be equivalently reformulated as a density functional theory [1], where the complete density functional is

 HSk =∫d3r{b02ρ2−b′02∑qρ2q+b1(ρτ−→j2)−b′1∑q(ρqτq−→j2q)+b22(→∇ρ)2−b′22∑q(→∇ρq)2 +~b1(→s⋅→T−∑ijJ2ij)+~b′1∑q(→sq⋅→Tq−∑ijJ2q;ij)+b33ρα+2−b′33ρα∑qρ2q −b4[ρ→∇⋅→J+→s⋅(→∇×→j)]−b′4∑q[ρq→∇⋅→Jq+→sq⋅(→∇×→jq)] +~b02→s2−~b′02∑q→s2q+~b22∑ij(∇isj)2−~b′22∑q∑ij(∇isj)2q+~b33ρα→s2−~b′33ρα∑q→s2q} (2)

where the last line contains the spin (time-odd) terms, which are usually omitted. However, they have quite important contribution for magnetic excitations [11], as will be illustrated below, so we use them in all calculations. Parameters depend on the parameters from (1):

 b0=t0(2+x0)2,b′0=t0(1+2x0)2,~b0=t0x02,~b′0=t02,b1=t1(2+x1)+t2(2+x2)8,b′1=t1(1+2x1)−t2(1+2x2)8,~b1=t1x1+t2x28,~b′1=−t1+t28,b2=3t1(2+x1)−t2(2+x2)16,b′2=3t1(1+2x1)+t2(1+2x2)16,~b2=3t1x1−t2x216,~b′2=3t1+t216,b3=t3(2+x3)8,b′3=t3(1+2x3)8,~b3=t3x38,~b′3=t38,b4=b′4=t42 (3)

Most Skyrme parametrization set explicitly and this fact is denoted as exclusion of the “tensor term”. There are also parametrizations fitted with the tensor term included, e.g. SGII [12], SLy7 [13].

The operators corresponding to the densities and currents in (2) are defined as

 (4)

and they are understood as single-particle operators in a many-body system. The densities are time-even and the currents are time-odd. Spin-orbital current and current have two indices, so we will decompose them into scalar, vector and (symmetric) tensor part

 ∑ijJ2ij=13∑iJ2ii+12→J2+2∑m=−2(−1)m[Jt]m[Jt]−m=13J2s+12→J2+J2t (5)

The coefficients related to the (non-relativistic) time-reversal symmetry of the operators are defined as

 |¯α⟩=(−1)lα+jα+mα|−α⟩,u(±)αβ=uαvβ±vαuβ,^T−1^A^T=γAT^A†,γAT={+1:time-even−1:time-odd (6)
 ^A=12∑αβu(γAT)αβ⟨α|^A|¯β⟩(−^α+α^α+β+γAT^α¯α^α¯β)=12∑αβu(γAT)αβ⟨¯α|^A|β⟩(^α+¯α^α+¯β−γAT^αα^αβ) (7)

where factors define the Bogoliubov transformation of particles to quasiparticles [4]. Pairing interaction in the particle-particle channel is defined as ( is the nucleon ground-state density)

 (8)

with constant parameters and . Contribution to the density functional is then

 Hpair=⟨BCS|^Vpair|BCS⟩=14∑q=p,nVq∫κ2q(→r)d3r,optionally with ×(1−ρ(→r)ρ0); (9)

and the pairing density with the corresponding operator is

 κq(→r) =mβ>0∑β∈q2fβuβvβψ†β(→r)ψβ(→r) (10) ^κ(→r) =∑β>0fβψ†β(→r)ψβ(→r)[2uβvβ(1−^α+β^αβ−^α+¯β^α¯β)+(u2β−v2β)(^α+β^α+¯β+^α¯β^αβ)] (11)

with energy-cutoff weight [14]. We use monopole pairing, and we involve only diagonal pairs.

## 3 RPA formalism in spherical symmetry

Notation of the Clebsch-Gordan coefficients comes from the book of Varshalovich [15]. In the language of spherical tensors, a hermitian operator satisfies

 ^A†m=(−1)m^A−m, (12)

for example a vector (rank-1) operator:

 ^A1=(−^Ax−i^Ay)/√2,^A0=^Az,^A−1=(^Ax−i^Ay)/√2. (13)

The position-dependent operators (4), in general denoted as , are as well expressed in terms of quasiparticles (7) and decomposed in a manner reminiscent of Wigner-Eckart theorem

 ^Jd(→r)=12∑αβLJMJJLd;αβ(r)(−1)lβ√2J+1CJMjαmαjβmβYL∗JM(ϑ,φ)(−^α+α^α+β+γdT^α¯α^α¯β) (14)

where the symbol in bold font represents the (complex-conjugated) spherical harmonics in its scalar, vector or tensor form (depending on the rank of operator : ):

 (15)

Expression (14) defines the reduced matrix elements , e.g., a standard density matrix element:

 ρLαβ(r)=u(+)αβ√(2jα+1)(2lα+1)(2jβ+1)(2lβ+1)Rα(r)Rβ(r)(−1)jβ+12√4π{lαlβLjβjα12}CL0lα0lβ0 (16)

and radial wavefunctions are taken from

 ⟨→r|α⟩=ψα(→r)=Rα(r)∑ν,sCjα,mαlα,ν,12,sYlα,ν(ϑ,φ)χs (17)

with spinors . Remaining reduced matrix elements will be published later.

Excitations of a given multipolarity are treated as the RPA phonons. One-phonon state, , is created by action of an operator on the RPA ground state , and has an excitation energy .

 ^C+ν|RPA⟩=|ν⟩,^Cν|% RPA⟩=0 (18)

The operator is a two-quasiparticle () operator defined by real coefficients

 ^C+ν=12∑αβCλνμνjαmαjβmβ(c(ν−)αβ^α+α^α+β+c(ν+)αβ^α¯α^α¯β) (19)

and satisfies the RPA equation

 [^H,^C+ν]2qp=Eν^C+ν (20)

where the index means that we take only the two-quasiparticle portion of the commutator (after normal ordering). Although all commutators should be evaluated in the RPA ground state, it is common to evaluate them in the HF+BCS ground state (i.e., we use quasi-boson approximation), since the contribution of and higher correlations in the ground state is assumed to be low.

Formulae (14) and (19) involve duplicate pairs. To remove them consistently, we will rescale diagonal pairing factors

and will be rescaled automatically. Diagonal matrix elements contribute only to electric transitions with even.

The hamiltonian is taken as a sum of mean-field part (HF+BCS) and the second functional derivative of Skyrme density functional (+ Coulomb and pairing interaction).

 (22)

RPA equation (20) then turns into a matrix equation

 (ABBA)(c(ν−)c(ν+))=(Eν00−Eν)(c(ν−)c(ν+)) (23)

which can be recast into half-dimensional symmetric-matrix eigenvalue problem [4]. Matrices and are

 App′ (24a) Bpp′ =∑dd′LγdT(−1)lβ+lδ2λ+1∫∞0δ2HδJdδJd′JλLd;p(r)JλL∗d′;p′(r)r2dr,(p≡αβ,p′≡γδ) (24b)

The expression is symbolical, and includes the integration of the delta function, yielding . The exchange Coulomb interaction can be treated by Slater approximation as a density functional

 Hxc=−34(3π)1/3e24πϵ0∫d3rρ4/3p(→r) (25)

however, the direct Coulomb interaction gives rise to a double integral instead

 ∑L∫∞0δ2HδJdδJd′JλL∗d;αβ(r)JλLd′;γδ(r)r2dr↦ e24πϵ04π2λ+1∫∞0r2dr ∫∞0r′2dr′ρλαβ(r)ρλγδ(r′)×{rλ/r′λ+1(r

After calculation of the RPA states, yielding and , we are interested in the matrix elements of electric and magnetic transition operators and in the transition densities and currents.

 ⟨ν|^Mλμ|RPA⟩ =⟨[^Cν,^Mλμ]⟩=∑α≥β(−1)lβ+1√2λ+1Mλ;αβ(c(ν−)αβ+γMTc(ν+)αβ)∗ (27) δρq;ν(→r) =⟨[^Cν,^ρq(→r)]⟩=∑α≥β(−1)lβ+1√2λ+1ρλq;αβ(r)(c(ν−)αβ+c(ν+)αβ)∗Y∗λμ(ϑ,φ) (28) δ→jq;ν(→r) =⟨[^Cν,^→jq(→r)]⟩=∑L∑α≥β(−1)lβ+1√2λ+1jλLq;αβ(r)(c(ν−)αβ−c(ν+)αβ)∗→YL∗λμ(ϑ,φ) (29)

Besides electric () and magnetic () operators in long-wave approximation, we implement also electric vortical, toroidal and compression operators [16] (toroidal operator constitutes the next-order term in long-wave expansion of the exact electric transition operator)

 ^MEλν (30a) ^MMλν =μNc√λ(2λ+1)∑q=p,n∑i∈q([gs,q2→σ+2gl,qλ+1^→L]rλ−1→Yλ−1λμ(ϑ,φ))i (30b) ^MEvor;λμ =−i/c2λ+3√2λ+1λ+1∫d3r^→jnuc(→r)rλ+1→Yλ+1λμ(ϑ,φ)=^MEtor;λμ+^MEcom;λμ (30c) ^MEtor;λμ =−12c(2λ+3)√λλ+1∫d3r^→jnuc(→r)⋅→∇×[rλ+2→Yλλμ(ϑ,φ)] (30d) ^MEcom;λμ =i2c(2λ+3)∫d3r^→jnuc(→r)⋅→∇[rλ+2Yλμ(ϑ,φ)](=−k^ME′com;λμ) (30e) ^ME′com;λμ =∑i^ME′com;λμ(→ri)=e2(2λ+3)∑q=p,nzq∑i∈q(rλ+2Yλμ(ϑ,φ))i (30f)

where are effective charges of the nucleons, are orbital/spin g-factors (we take ; spin g-factors are reduced by a quenching factor ) and is the nuclear current composed of convective and magnetization part

 ^→jnuc(→r)=eℏmp∑q=p,n∑i∈q[zq^→ji(→r)+14gs,q→∇×^→si(→r)] (31)

where (convective) current and spin one-body operators are the same as in Skyrme functional (4).

## 4 RPA formalism in axial symmetry

Axial coordinates are

 ϱ=√x2+y2,z,φ;x=ϱcosφ, y=ϱsinφ (32)

Calculations in axially deformed nuclei don’t conserve total angular momentum (previously denoted or ), nevertheless, they conserve its -projection () and parity, so it is convenient to preserve part of the formalism from spherical symmetry, namely the convention of -components for vector and tensor operators, and their hermitian conjugation (12). Operators of differentiation are then

 ∇±1=1√2(∓∂∂x−i∂∂y)=∓e±iφ√2(∂∂ϱ±iϱ∂∂φ),∇0=∂∂z (33)

Single-particle wavefunction (and its time-reversal conjugate) is expressed as a spinor

 (34)

and the radial parts of its derivatives will be denoted by a shorthand notation

 (35)

Radial functions are real, and their spinor-wise products will be denoted by a dot to keep the expressions simple:

 Rα⋅Rβ ≡ Rα↑(ϱ,z)Rβ↑(ϱ,z)+Rα↓(ϱ,z)Rβ↓(ϱ,z) (36)

Vector currents will be decomposed in the style of rank-1 tensor operators. Vector product in the expression for spin-orbital current leads to (for vector product in -scheme see [15, (1.2.28)])

 (37)

Matrix elements of densities and currents are then

 ⟨α|^ρ|β⟩ =Rα⋅Rβei(mβ−mα)φ (38a) ⟨α|^τ|β⟩ =(R(0)α⋅R(0)β+R(+)α⋅R(+)β+R(−)α⋅R(−)β)ei(mβ−mα)φ (38b) Factor ei(mβ−mα)φ will be omitted in the following expressions. ⟨α|→J|β⟩ (38c) ⟨α|→j|β⟩ ={+1:i2eiφ(−R(−)α⋅Rβ−Rα⋅R(+)β)0:i2(R(0)α⋅Rβ−Rα⋅R(0)β)−1:i2e−iφ(−R(+)α⋅Rβ−Rα⋅R(−)β)⟨α|→s|β⟩={+1:eiφ(−√2Rα↑Rβ↓)0:Rα↑Rβ↑−Rα↓Rβ↓−1:e−iφ(√2Rα↓Rβ↑) (38d) ⟨α|→T|β⟩ ={+1:eiφ(−√2)[R(0)α↑R(0)β↓+R(+)α↑R(+)β↓+R(−)α↑R(−)β↓]0:R(0)α↑R(0)β↑−R(0)α↓R(0)β↓+R(+)α↑R(+)β↑−R(+)α↓R(+)β↓+R(−)α↑R(−)β↑−R(−)α↓R(−)β↓−1:e−iφ√2[R(0)α↓R(0)β↑+R(+)α↓R(+)β↑+R(−)α↓R(−)β↑] (38e) ⟨α|→∇×→j|β⟩ =−i(→∇ψα)†×→∇ψβ={+1:eiφ(R(−)α⋅R(0)β+R(0)α⋅R(+)β)0:R(−)α⋅R(−)β−R(+)α⋅R(+)β−1:e−iφ(−R(+)α⋅R(0)β−R(0)α⋅R(−)β) (38f) ⟨α|→∇⋅→J|β⟩ =−R(+)α↑R(+)β↑+R(+)α↓R(+)β↓+R(−)α↑R(−)β↑−R(−)α↓R(−)β↓ −√2(R(0)α↑R(−)β↓+R(−)α↓R(0)β↑+R(0)α↓R(+)β↑+R(+)α↑R(0)β↓) (38g) ⟨α|Js|β⟩ =i2[(R(0)α↑+√2R(−)α↓)Rβ↑−(R(0)α↓+√2R(+)α↑)Rβ↓ −Rα↑(R(0)β↑+√2R(−)α↓)+Rα↓(R(0)β↓+√2R(+)β↑)] (38h) ⟨α|Jt|β⟩ (38i) ^→Lψβ =−i(→r×→∇)ψβ={+1:eiφ(ϱ√2R(0)β+zR(+)β)0:ϱ√2(R(+)β+R(−)β)−1:e−iφ(ϱ√2R(0)β−zR(−)β) (38j)

In the actual calculation, it is necessary to choose projection of angular momentum and parity (together denoted also as , and ). Transition operators have the form of

 ^Mλμ=∑iMλμ(ϱi,zi)eiμφi (39)

where contains a function (or even derivatives) not dependent on . Choice of the two-quasiparticle pairs is restricted by (time-reversed states are omitted from the discussion here for simplicity; in short, it is necessary to add pairs with during the duplicate removal similar to (21)). Commutators are then evaluated in quasiparticle vacuum as

 ⟨[^A†,^B]⟩=γAT−γBT2∑αβ