UUser guide for the hfbcs-qrpa (v1) code
G. Col`o ∗ and X. Roca-Maza † Dipartimento di Fisica, Universit`a degli Studi di Milano and INFN,Sezione di Milano, Via Celoria 16, 20133 Milano, Italy. (Dated: February 15, 2021)We briefly explain the structure of the Hartree-Fock-Bardeen-Cooper-Schrieffer (HFBCS) andQuasi-particle Random Phase Approximaiton (QRPA) code hfbcs-qrpa , highlighting only the dif-ferences with the code skyrme rpa that has been previously published [1]. The code deals with openshell spherical systems and non-spin flip, non-charge-exchange and natural parity ( − l excitations. I. DIFFERENCES WITH THE PUBLISHEDCODE
SKYRME RPA
The present code is an extension of the Hartree-Fockplus Random Phase Approximation (HF plus RPA) codethat has already been published in Ref. [1]. HF is ex-tended to HF-Bardeen-Cooper-Schrieffer (HF-BCS) and,accordingly, RPA to Quasiparticle RPA (QRPA). Thescheme remains self-consistent, once the effective pairingforce is introduced.Along the code we use (cid:126) c = 197 . mc = 938 . m p + m n ) /
2, forboth neutrons and protons.
A. HF-BCS
HF equations are solved in coordinate space, assumingspherical symmetry, as explained in [1], and employingthe standard form of a Skyrme interaction . AvailableSkyrme interactions are: SkM*, SkP, SLy4, SLy5, SkX,KDE0-J33 (KDE33) and SAMi. The HF equations arecoupled to the standard BCS equations, that in sphericalsymmetry with α = ( n, l, j ) read (cid:88) α (2 j α + 1) v α = N, (1)∆ α = − (cid:88) α (cid:48) ∆ α (cid:48) E α (cid:48) V α ˜ αα (cid:48) ˜ α (cid:48) , (2)where E , u and v are the usual quasi-particle energiesand BCS amplitudes as defined in the literature [3], andthe tilde denotes the time-reversal state. N is the par-ticle number. Either proton pairing, or neutron pairing,or both, can be active. No neutron-proton pairing isincluded and, thus, the particles retain their quantumnumber t z .To keep u and v positive, we use the phase conventionin which the single-particle (s.p.) wave functions include ∗ [email protected] † [email protected] In the current version, tensor terms [2] are implemented at theHF-BCS level but not in the residual QRPA interaction. a factor i l [Eq. (7) of Ref. [1] should be modified ac-cordingly]. The number and gap equations, (1) and (2),are solved in a limited window which is specified by theuser if the pairing is run in manual mode. In that case,an upper limit on the corresponding s.p. energy ε α anda fixed number of maximum s.p. levels above the Fermienergy to be considered should be specified. Note thatonly one of those two pairing cutoffs, namely the low-est, will act as the real one. These two possibilities givemore flexibility in the usage of the code. Both types ofcutoffs can be different for neutrons and protons. Thereis no lower limit imposed for this pairing window. Forthose users not conversant with pairing calculations ornot interested in exploring pairing effects, we stronglyrecommend to use pairing in automatic mode, in whichcase, a volume pairing force (see below) has been fitted toreproduce the experimental neutron pairing gap in Sn(∆ = 1 .
245 MeV) with a neutron pairing window fixedby taking into account only bound states above the neu-tron Fermi energy and to a maximum of 6 neutron s.p.levels. The default pairing window in automatic mode isset in an analogous way.The matrix elements V α ˜ αα (cid:48) ˜ α (cid:48) are calculated using azero-range, density-dependent pairing force of the type V ( r , r ) = V (cid:20) x (cid:18) ρ ( r + r ) ρ (cid:19)(cid:21) δ ( r − r ) ≡ V pair ( r ) δ ( r − r ) . (3)In the above equation, ρ = 0 .
16 fm − while the othertwo parameters x and V have to be adjusted. x = 0 , . V is usually adjusted to reproduce some given pair-ing gap but it has to be noted that the code does notmake this adjustement by itself unless automatic modeis choseen where only the use of volume pairing ( x = 0)is allowed.The total HF-BCS energy can be calculated, and is cal-culated by the code, in two different manners. It can becalculated from the force, or energy functional, directly.In this case E = E kin + E Skyrme + E Coul + E pair . (4)The first two terms are the kinetic and Skyrme contri-butions to the energy and are written as the volume in-tegrals of the corresponding energy densities, E kin and a r X i v : . [ nu c l - t h ] F e b E Skyrme , as in Eqs. (2.1) and (2.2) of Ref. [4], or in Eqs.(24) and (38) of Ref. [5]. The Coulomb energy has adirect part and an exchange part that is calculated inthe usual Slater approximation. The pairing energy iswritten as in Eq. (16) of [5]. As an alternative, E = 12 (cid:88) α v α ε α + 12 E kin + E rearr + E pair , (5)where the first term includes the s.p. energies ε α , and E rearr is the rearrangement term associated with thedensity-dependent term in the Skyrme force (as in Eq.(109) of [5]), as well as with the Coulomb exchange en-ergy within the Slater approximation. The equality be-tween the two expressions (4) and (5) for the total energyis, among the rest, a useful test for the convergence of theHF-BCS equations.Regarding the calculation of charge radii, we take intoaccount the electromagnetic finite size of the proton andthe neutron. That is, we assume the proton and neutronelectric sizes as (cid:104) r p (cid:105) = 0 .
64 fm and (cid:104) r n (cid:105) = − .
11 fm ,respectively, as well as the anomalous proton magneticmoment κ p = µ p − µ p = 2 .
79 and neutron mag-netic moment κ n = µ n = − .
91. The expression for thecharge radius (cid:104) r (cid:105) ch used is derived in the non-relativisticlimit and can be written as [6] (see also [7]), (cid:104) r (cid:105) ch = (cid:104) r (cid:105) p + (cid:104) r p (cid:105) + NZ (cid:104) r n (cid:105) + 1 Z (cid:18) (cid:126) mc (cid:19) (cid:88) ατ v ατ (2 j α + 1) κ τ (cid:104) σ · l (cid:105) α , (6)where τ = p, n indicates if the nucleon is a proton or aneutron, respectively, and (cid:104) r (cid:105) p is the mean square radius of the density distribution of protons, (cid:104) r (cid:105) p = (cid:90) d r r ρ p ( r ) . (7)We also evaluate the isospin mixing ε in the obtaineds.p. wave-functions as discussed in Sec. IIA of Ref. [8].Specifically, ε ≡ N − Z + 2 (cid:88) n p ,n n l p = l n j p = j n (2 j p + 1) v p u n O np (8)where we have omitted the α label (previously defined)for the sake of simplicity and clarity. The overlap factor O np between the neutron ( n ) and proton ( p ) radial part R n,l ( r ) of the considered s.p. wave function is defined as O np ≡ (cid:90) ∞ drr R n p ,l p ( r ) R n n ,l n ( r ) . (9) B. QRPA
The QRPA equations are (cid:18) A αβ,γδ B αβ,γδ − B ∗ αβ,γδ − A ∗ αβ,γδ (cid:19) (cid:18) X νγδ Y νγδ (cid:19) = E ν (cid:18) X ναβ Y ναβ (cid:19) . (10)The basis on which these equations are written is madeup with two quasiparticle states, like ( αβ ) and ( γδ ).These respect the usual selection rules in that they arecoupled to good angular momentum J and parity π . Theupper limit on the two quasiparticle energy can be spec-ified (see below). On the other hand, the code automat-ically discards configurations for which both quasiparti-cle states are fully occupied or both have an occupationprobability lower than 10 − .The matrix elements A and B , in the angular momentum-coupled representation, written on the HF-BCS basistwo-quasiparticle basis, have the form A αβ,γδ = 1 (cid:112) δ αβ (cid:112) δ γδ × [( E α + E β ) δ αγ δ βδ + G αβγδ ( u α u β u γ u δ + v α v β v γ v δ )+ F αβγδ ( u α v β u γ v δ + v α u β v γ u δ ) − ( − j γ + j δ − J (cid:48) F αβδγ ( u α v β v γ u δ + v α u β u γ v δ ) (cid:105) , (11) B αβ,γδ = 1 (cid:112) δ αβ (cid:112) δ γδ × (cid:104) − G αβδγ ( u α u β v γ v δ + v α v β u γ u δ ) − ( − j δ + j γ − J (cid:48) F αβδγ ( u α v β u γ v δ + v α u β v γ u δ )+( − j α + j β + j γ + j δ − J − J (cid:48) F αβγδ ( u α v β v γ u δ + v α u β u γ v δ ) (cid:105) , (12)with G αβγδ = (cid:88) m α m β m γ m δ (cid:104) j α m α j β m β | JM (cid:105)(cid:104) j γ m γ j δ m δ | J (cid:48) M (cid:48) (cid:105) V ppαβ,γδ , (13) F αβγδ = (cid:88) m α m β m γ m δ (cid:104) j α m α j β m β | JM (cid:105)(cid:104) j γ m γ j δ m δ | J (cid:48) M (cid:48) (cid:105) V phα ¯ δ ¯ βγ . (14) V phα ¯ δ ¯ βγ are the (uncoupled) matrix elements of theparticle-hole effective interaction, and V ppαβ,γδ representthe (uncoupled) matrix elements of the particle-particleeffective interaction. The p-h effective interaction is de- scribed in detail in Ref. [1]. In fact, one can immediatelyspot that the matrix element in Eq. (14) is the same thathas been defined in Eq. (16) of Ref. [1]; the different no-tation here has the purpose of distinguishing it clearlyfrom the matrix element G that is the one of the pairingforce and reads G αβγδ = R (cid:88) λ (cid:18) ( − ) j β + j γ + J (cid:26) j α j γ λj δ j β J (cid:27) (cid:104) α || Y λ || γ (cid:105)(cid:104) β || Y λ || δ (cid:105) + ( − ) j β + j γ (cid:26) j α j δ λj γ j β J (cid:27) (cid:104) α || Y λ || δ (cid:105)(cid:104) β || Y λ || γ (cid:105) (cid:19) , (15)where R is the radial integral R = (cid:90) drr u α ( r ) u β ( r ) u γ ( r ) u δ ( r ) V pair ( r ) . (16)The excitation operators and associated sum rules usedin the code are those discussed in Ref. [1]. Regardingspurious states: the HF-BCS is known to break differ-ent symmetries present in the original Hamiltonian. Inprinciple QRPA should restore those symmetries exactly.However, numerical implementations are not exact, and spurious states do not appear necessarily at zero energyas expected. As a rule, the sppurious state should bethe first excited state, close to zero energy, in the cal-culated spectrum. Specifically, within the present code:i) for 0 + excitations, we have not implemented any pro-jection technique to avoid the spurious state originatedfrom the violation of the particle number; and ii) for1 − excitations, we use the operators defined in [1] to ap-proximately remove the spurous state that arise from theviolation of translational invariance. II. INPUT OF THE CODE
In this section we briefly describe the input file hfbcs-qrpa.in needed to run the code. The format of the inputfile should be as follows:
Line Input Explanation---- ----- -----------1 FORCE Name of the Skyrme interaction to be used[options: SkM*, SkP, SLy4, SLy5, SkX, KDE33, SAMi]2 MESH STEP Number of mesh points and step (fm)3 A Z Mass number (A) and proton number (Z)4 EPS Required accuracy for the HFBCS equations. Specifically, it correspondsto the maximum accepted difference between consecutive iterations and itis used for: single-particle energy, gap and chemical potential.5 FLAG BCSP BCSN Automatic (0) or manual(1) BCS calculaitons. Active (1) non-active (0)for protons. Active (1) non-active (0) for neutrons.Only if manual mode has been chosen next three lines (6-8) should be added:6 LEVP LEVN Number of levels taken above the Fermi level in the BCS calculation forprotons and neutrons.7 EMAXP EMAXN Maximum proton and neutron energies for the BCS calculation.8 -V_0 x Parameters of the pairing interaction to be used defined inthe previous section.9 CALEX Type of calculations [options: - UNP: stands for unperturbed calculation, residual interaction set to zero- TDA: Tamm-Dancoff Approximation is used- RPA: Random Pahse Approximation is used].10 J P Total angular momentum and pairty (+1 or -1). At the moment only naturalparity is considered.11 ECP ECN Minimum and maximum energy for the considered two bodyneutron-neutron and proton-proton basis states.
III. OUTPUT OF THE CODE
The code produces different output files. Definitions of the different quantities such as operators or transitiondensities are as in skyrme rpa [1]. Below a brief description:
File name Description--------- -----------hfbcs-qrpa.out main output of the codequasiparticle_states.out full list of quasi-particle states used in the calculation.Information on the Fermi energies is also givenbasis_states.out full list of two-body neutron-neutron and proton-protonbasis states used in the calculation.d.out proton and neutron densitiestd.out proton and neutron transition densities for all excited statesstrength_*.out three files contain isoscalar, isovector and electromagneticstrength functions (reduced transition probabilities convolutedby a Lorenzian function with a energy step of 0.1 MeV and widthof 1 MeV)
In the main output file ( hfbcs-qrpa.out ) details on the parameters of the used functional as well as active terms areprinted. One-body center of mass correction is active for all the Skyrme interaction assuming the same prescriptionas in skyrme rpa . Convergence information and a list of quasi-particle levels is given. Total binding energy aswell as its different contributions are printed in two different forms (see the discussion above). Firstly, the energy iscalculated as the integral of the energy density. Then, the energy is calculated via the sum of the eigenvalues. In thislatter case the rearrangement energy is needed. The calculation of the binding energy in two different ways servesas a strong test for the convergence the code. Neutron and proton root mean square radii are also printed and theroot mean square radius is estimated. Average pairing gaps, and isospin mixing in the groud state wave function arealso given in the output. Finally, different sum rules and isoscalar, isovector and electromagnetic reduced transitionprobabilities are listed.
IV. EXAMPLE
In this section, we show an input ( hfbcs-qrpa.in ) example of the code which uses the manual pairing option(FLAG=1):
SLy5200 0.1120 501.d-61 0 10 60. 0.870.6 1.RPA2 +10. 100.
The main output ( hfbcs-qrpa.out ) of the code corresponding to this input is given below. . . .>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>READING SKYRME INTERACTION TO BE USED>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>NAME: SLy5SKYRME FORCE PARAMETERS IN STANDARD FORM:t0 = -2484.88000x0 = 0.77800t1 = 483.13000x1 = -0.32800t2 = -549.40000x2 = -1.00000t3 = 13763.00000x3 = 1.26700g = 0.16667W0 = 126.00000W0p = 126.00000PARAMETERS OF THE SKYRME ENERGY DENSITY FUNCTIONAL:as defined in J. Dobaczewski and J. Dudek,Phys. Rev. C52, 1827 (1995)C_0^\rho = -931.83000 + 860.18750 *\rho^ 0.16667C_1^\rho = 793.91916 + -1013.30087 *\rho^ 0.16667C_0^\delta\rho = -76.52453C_1^\delta\rho = 16.37485C_0^\tau = 56.24937C_1^\tau = 23.95020C_0^s = -172.69916 + 439.84254 *\rho^ 0.16667C_1^s = 310.61000 + -286.72917 *\rho^ 0.16667C_0^\nabla J = -94.50000C_1^\nabla J = -31.50000C_0^\Delta s = 46.08734C_1^\Delta s = 14.06234C_0^T = -15.66645C_1^T = -64.53312>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>HF+BCS CALCULATION IN A BOX>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>NUCLEUS:-Mass number 120.-Proton number 50.BOX PARAMETERS:-N. of points 200-Step size 0.10E+00CONVERGENCE PARAMETERS:-Error spe 0.10E-05-Error gap eq. 0.10E-05-Error num. eq. 0.10E-05
POTENTIAL INCLUDED IN THE CALCULATIONS:-SKYRME STANDARD WITH|___ J2 TERMS... YES|___ TENSOR..... NO-COULOMB DIRECT....... YES-COULOMB EXCHANGE..... YES-SPIN-ORBIT........... YES-DEFAULT PAIRING...... NO-PROT. PAIRING (BCS).. NO-NEUT. PAIRING (BCS).. YESENERGY CUT..... 0.000000000000000E+000LEV. ABOVE F... 6Type: delta forceVp= V0(1-x(rho/rho0))^gV0 = 0.87060E+03x = 0.10000E+01rho0= 0.16000E+00g = 0.10000E+01-CENTER OF MASS CORR.. YESCONVERGENCE:-At iteration number: 77-Maximum error in spe: 0.99830E-05>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>SINGLE PARTICLE ENERGIES, OCCUP., V^2 AND GAPS>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>PROTON STATES1S1/2 E=-0.4664E+02 DEG= 0.2000E+01 V^2= 0.1000E+01 D= 0.00001P3/2 E=-0.3873E+02 DEG= 0.4000E+01 V^2= 0.1000E+01 D= 0.00001P1/2 E=-0.3778E+02 DEG= 0.2000E+01 V^2= 0.1000E+01 D= 0.00001D5/2 E=-0.2980E+02 DEG= 0.6000E+01 V^2= 0.1000E+01 D= 0.00001D3/2 E=-0.2756E+02 DEG= 0.4000E+01 V^2= 0.1000E+01 D= 0.00002S1/2 E=-0.2559E+02 DEG= 0.2000E+01 V^2= 0.1000E+01 D= 0.00001F7/2 E=-0.2040E+02 DEG= 0.8000E+01 V^2= 0.1000E+01 D= 0.00001F5/2 E=-0.1634E+02 DEG= 0.6000E+01 V^2= 0.1000E+01 D= 0.00002P3/2 E=-0.1522E+02 DEG= 0.4000E+01 V^2= 0.1000E+01 D= 0.00002P1/2 E=-0.1361E+02 DEG= 0.2000E+01 V^2= 0.1000E+01 D= 0.00001G9/2 E=-0.1083E+02 DEG= 0.1000E+02 V^2= 0.1000E+01 D= 0.0000NEUTRON STATES1S1/2 E=-0.5616E+02 DEG= 0.2000E+01 V^2= 0.1000E+01 D= 0.20171P3/2 E=-0.4748E+02 DEG= 0.4000E+01 V^2= 0.1000E+01 D= 0.46781P1/2 E=-0.4615E+02 DEG= 0.2000E+01 V^2= 0.1000E+01 D= 0.41911D5/2 E=-0.3787E+02 DEG= 0.5999E+01 V^2= 0.9998E+00 D= 0.79631D3/2 E=-0.3481E+02 DEG= 0.3999E+01 V^2= 0.9998E+00 D= 0.70452S1/2 E=-0.3334E+02 DEG= 0.1999E+01 V^2= 0.9997E+00 D= 0.85191F7/2 E=-0.2775E+02 DEG= 0.7993E+01 V^2= 0.9991E+00 D= 1.15551F5/2 E=-0.2250E+02 DEG= 0.5992E+01 V^2= 0.9987E+00 D= 1.04692P3/2 E=-0.2194E+02 DEG= 0.3993E+01 V^2= 0.9982E+00 D= 1.18192P1/2 E=-0.2013E+02 DEG= 0.1995E+01 V^2= 0.9975E+00 D= 1.18661G9/2 E=-0.1747E+02 DEG= 0.9937E+01 V^2= 0.9937E+00 D= 1.48132D5/2 E=-0.1137E+02 DEG= 0.5787E+01 V^2= 0.9644E+00 D= 1.25431G7/2 E=-0.9976E+01 DEG= 0.7133E+01 V^2= 0.8916E+00 D= 1.39243S1/2 E=-0.9077E+01 DEG= 0.1618E+01 V^2= 0.8089E+00 D= 1.0877 |HF> = sqrt(1-e^2)|T0, T0> + e|T0+1, T0>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>e^2 (%) = 0.31506E+00>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>LINEAR RESPONSE EVALUATED INRANDOM PHASE APPROXIMATION>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>CONSIDERED TRANSITIONS WITH:J 2PARITY 1>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>MOMENTS OF THE STRENGTH FUNCTION>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>-NON-ENERGY WEIGHTED SUM RULE:m(0)(IS) = 0.26417E+05m(0)(IV) = 0.11253E+05-INVERSE ENERGY WEIGHTED SUM RULE:m(-1)(IS) = 0.10068E+05m(-1)(IV) = 0.14632E+04-ENERGY WEIGHTED SUM RULE:m( 1)(IS) = 0.20875E+06m( 1)(IV) = 0.23841E+06m(1)[D.C.]IS = 0.21407E+06Exhaus. IS(%)= 97.51375IV DOUBLE COMMUTATOR EWSR NOT AVAILABLE-ENERGY CUBE WEIGHTED SUM RULE:m( 3)(IS) = 0.42103E+08m( 3)(IV) = 0.16041E+09>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>AVERAGE EXCITATION ENERGIES>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>E_CONSTR. (IS) = 0.45533E+01E_CENTR. (IS) = 0.79021E+01E_SCAL. (IS) = 0.14202E+02E_CONSTR. (IV) = 0.12765E+02E_CENTR. (IV) = 0.21187E+02E_SCAL. (IV) = 0.25939E+02>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>REDUCED TRANSITION PROBABILITIES>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>E BEL_is FRAC_NEWSR BEL_em BEL_iv FRAC_NEWSR------------------------------------------------------------------------------1 0.13051E+01 0.11155E+05 0.42229E+02 0.12264E+04 0.12658E+04 0.11249E+022 0.27263E+01 0.29313E+02 0.11096E+00 0.41912E+00 0.16969E+02 0.15080E+003 0.30799E+01 0.57219E+03 0.21660E+01 0.87101E+02 0.27614E+02 0.24539E+00
Finally, for this example, we show in Fig. 1 the reduced transition probabilities contained in hfbcs-qrpa.out andtheir convolution with a Lorenzian function with a width of 1 MeV. The latter information is also given by thecode and can be found in the file strength is.out for the isoscalar response and strength iv.out for the isovectorresponse. − x B ( E , I S ) [f m ] E [MeV]0246810 − x B ( E , I S ) [f m ] SLy5
Sn PairingNo pairing
FIG. 1. Reduced transition probabilities (bars) for the text example contained in hfbcs-qrpa.out and their convolution with aLorenzian function (solid lines) with a width of 1 MeV. The latter information is also given by the code and can be found in thefile strength is.out for the isoscalar response (upper panel) and strength iv.out for the isovector response (lower panel). [1] G. Col`o, L. Cao, N. Van Giai, and L. Capelli, Computer Physics Communications , 142 (2013).[2] F. Stancu, D. Brink, and H. Flocard, Physics Letters B , 108 (1977).[3] P. Ring and P. Schuck, “The nuclear many-body problem,” (Springer, 1980).[4] E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer, Nuclear Physics A , 231 (1998).[5] W. Ryssens, V. Hellemans, M. Bender, and P.-H. Heenen, Computer Physics Communications , 175 (2015).[6] C. J. Horowitz and J. Piekarewicz, Phys. Rev. C , 045503 (2012).[7] E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer, Nuclear Physics A , 710 (1997).[8] X. Roca-Maza, H. Sagawa, and G. Col`o, Phys. Rev. C102