SELF-CONSISTENT METHODS IN NUCLEAR STRUCTURE PHYSICS 8 JACEK DOBACZEWSKI 9 Institute of Theoretical Physics, Warsaw University 9 ul. Hoz˙a 69, PL-00681 Warsaw, Poland 1 E-mail: [email protected] n a Joint Institute for Heavy Ion Research, Oak Ridge National Laboratory, J P.O. Box 2008, Oak Ridge, TN 37831, U.S.A. 7 2 Department of Physics, University of Tennessee, Knoxville, TN 37996, U.S.A. 1 v WepresentaverybriefdescriptionoftheHartree-Fockmethodinnuclearstructure 6 physics,discussthenumericalmethodsusedtosolvetheself-consistentequations, 5 and analyze the precision and convergence properties of solutions. As an appli- 0 cation we present results pertaining to quadrupole moments and single-particle quadrupolepolarizationsinsuperdeformednucleiwithA∼60. 1 0 8 1 Introduction 9 / h Self-consistent methods have been used in the low-energy nuclear structure t studies over many years, and represent a mature field with numerous suc- - l cessful applications.1,2,3,4,5 A number of computer codes solving the nuclear c u Hartree-Fock (HF) problem have already been developed. Two types of effec- n tive nucleon-nucleon interactions have been mainly employed. Starting with : the work of Vautherin and Brink6 many authors have applied the nuclear v i HF theory with the Skyrme effective interaction, while the work of Gogny7,8 X initiated numerous studies with the force which carries his name. r The methods employed to solve the HF equations depend mainly on the a effective force used and on the assumed symmetries of the many-body wave functions. For the solutions which allow at least triaxial deformations, two different methods have been applied for the two above mentioned effective in- teractions. The first one, used in conjuncture with the Skyrme interaction, is formulated in the spatial coordinates and makes use of the finite-difference,9 orFourier,10 or spline-collocation11,12 methods to approximatedifferentialop- erators. The solution is then obtained by using the imaginary time evolution operator.13 The second one, used for the finite-range Gogny interaction, employs a truncatedharmonicoscillator(HO)basis14,15 andsolvestheproblemeitherby aniterativediagonalizationofthemean-fieldHamiltonian,orbythegradient16 1 ortheconjugategradient17methods. Recently,themethodwhichincorporates the advantages of both existing approaches, and combines the robustness of theCartesianHObasiswiththesimplicityoftheSkyrmeinteraction,hasbeen implemented.18 The methods using spatial coordinates have several advantages. First of all,variousnuclearshapes canbe easilytreatedonthe samefooting; the same cubic lattice of points in three spatial dimensions is suitable to accommodate wave functions with, in principle, arbitrary deformations restricted only by a specific symmetrization of the lattice. This allows easy studies of systems for which the deformation is not a priori known, or is ill defined because of deformation instabilities or a shape coexistence. Secondly, using spatial coordinates allows studies which address the asymptotic form of nucleonic wave functions at large distances. This is particularly important for a precise description of weakly bound nuclei, where the use of spatial coordinates is a necessity.19 Third, for the Skyrme zero-range interaction, the mean fields are local(apartfroma velocitydependence) andcanbe easilyprogrammedin the spatialcoordinates. Lastbutnotleast,thetreatmentofwavefunctionsonlarge lattices(12 12 12isatypicalexample)is easilyamenabletovectorizationor × × parallelization of the algorithm. Methods using the HO basis have other advantages. Firstly, the basis provides a natural cut-off for many operators which otherwise are unbound and require particularly delicate treatment in the spatial coordinates. This concerns in particular the multipole moment and the angular momentum op- erators which are often used as constraining operators. For the corresponding constraints the solutions can become unstable when the non-zero probability amplitudes (wave functions) move towards large distances as it is the case for e.g. weakly bound nucleons. Secondly, much smaller spaces are usually requiredtodescribethenuclearwavefunctionswithinagivenprecision. Typi- callyabasisofabout300HOwavefunctionsissufficientformostapplications. Third,the iterativediagonalizationofthe meanfield Hamiltoniancanbe used to find the self-consistent solutions, which provides a rapidly converging al- gorithm, and, last but not least, scalar or superscalar computers can also be used, because the typical sizes of the information handled is smaller and the performance is less dependent on the use of a vector processor. Detaileddiscussionofthe stability,convergence,andefficiencyofmethods using the expansion on the Cartesian HO basis has recently been presented together with the description of the hfodd code.18 In the present communi- cation we give a brief r´esum´e of the methods employed (Sec. 2), present tests of the precision and discuss the CPU times required (Sec. 3), and illustrate possible applications by presenting a few recentresults obtained for the A 60 ∼ 2 nuclei (Sec. 4). 2 Hartree-Fock method ′ The eigenequations for the HF single-particle Routhians h are called the HF α equations, h′ ψ (rσ)=e′ ψ (rσ), (1) α i,α i,α i,α where the neutron (α=n) and proton (α=p) Routhians read ¯h2 h′ = ∆+Γ +δ UCoul+Umult ω Jˆ. (2) α −2m α αp − y y Detailed expressions for the nuclear mean field operators Γ , Coulomb poten- α tial UCoul, and multipole constraint potential Umult are given in Ref.18 After solving Eq. (1) one calculates the total energy of the system as a sum of the kinetic, Skyrme, Coulomb, and pairing terms = kin+ Skyrme+ Coul+ pair. (3) E E E E E ′ Thesameenergycanbe calculatedbyusingtheRouthianeigenvaluese and i,α occupation probabilities v2 as i,α ˜= 1 v2 e′ + 1 kin+ pair rear+ 1 Coul mult cran. (4) E 2X i,α i,α 2E E −E 3Eexch −Ecorr −Ecorr i,α Here, the rearrangement energy rear results from the density dependence of E the Skyrme interaction. Similarly, the Coulombexchange energy Coul can be Eexch considered as resulting from a zero-order interaction term depending on the density as ρ−2/3. The remaining terms correct for the fact that Routhians p contain constraint terms which are not present in the total energy. The difference between the energies ˜and is exactly equal to zero when E E thedensitiesandfieldsdonotchangefromoneiterationtothenextone. Hence their difference δ = ˜ (5) E E −E provides a useful measure of the quality of convergence, and is called the sta- bility of the HF energy. The HF equations(1)canbe solvedbyexpanding the single-particlewave functions ψ (rσ) onto the deformed HO wave functions ψ (rσ) in the i nxnynz,sz Cartesian coordinates, i.e., Nx Ny Nz ψ (rσ)= Anxnynz,szψ (rσ). (6) i X X X X i nxnynz,sz nx=0ny=0nz=0sz=−21,12 3 Here N , N , and N are the maximum numbers of the HO quanta corre- x y z sponding to the three Cartesian directions. However, the sums over n , n , x y andn areperformedoverthegridofpointswhichformapyramidratherthan z a cube. The HO wave functions have the standard form ψ (rσ)=ψ (x)ψ (y)ψ (z)δ , (7) nxnynz,sz nx ny nz szσ where ψnµ(xµ)=bµ12Hn(0µ)(ξµ)e−21ξµ2, (8) and ξ =b x are dimensionless variables scaled by the oscillator constants µ µ µ b = mω /¯h. (9) µ q µ PolynomialsH(0)(ξ)areproportionaltothestandardHermiteorthogonalpoly- n nomials H (ξ),20 n −1 H(0)(ξ)= √π2nn! 2 H (ξ). (10) n (cid:0) (cid:1) n Accuracy of the solution of the HF equations with the wave functions ex- pandedontotheCartesianHObasis,Eq.(6),dependsonthethreeparameters ¯hω , h¯ω , h¯ω defining the HO frequencies in three Cartesian directions, and x y z on the number M of the HO states included in the basis. In the code hfodd we use the standard prescription21,22 to chose the HO states included in the basis, namely, the M states with the lowest HO single-particle energies, ǫ =h¯ω (n + 1)+h¯ω (n + 1)+h¯ω (n + 1), (11) nxnynz x x 2 y y 2 z z 2 areselectedamongthosewhichhaven N ,n N ,andn N ,whereN is x 0 y 0 z 0 0 ≤ ≤ ≤ the fixed maximum number of HO quanta. It should be noted that in general both M and N have to be specified to define the basis. Only for large N , 0 0 the basis is defined solely by M and does not depend on N . In this case, 0 the grid of points (n ,n ,n ) defining the states included in the basis forms a x y z pyramidinthreedimensions,withthe inclinedfacedelimitedbythe condition ǫ const. On the other hand, only for small values of N the basis is nxnynz≤ 0 defined solely by N and does not depend on the energy cut-off. In this case 0 the corresponding grid of points n n n forms a cube of the size N . In all x y z 0 intermediate cases the shape of the basis corresponds to a pyramid with the corners cut off, or to a cube with the corners cut off. Usually N is chosen 0 large enough so that all the states allowed by the energy cut-off are included in the basis. 4 3 Precision and convergence TheHFequations(1)canbesolvedbyusingstandarditerativemethods1which arebasedontheself-consistencyprinciple. Namely,theself-consistentsolution ′ is attained when the mean fields Γ appearing in the mean-field Routhian h α α equal to the mean fields calculated from the occupied single-particle states ψ (rσ). Therefore,by iteratively calculating the meanfields and using them i,α toobtainnewsingle-particlestatesonemayeventuallyreachtheself-consistent solution. A convergence of such a procedure is not guaranteed and, in fact, in nuclear applications the method such as formulated above almost always diverges. Thereasonforthatisthefactthattheconsecutivecorrectionsofthe mean fields turn out to be too violent and lead to too strong modifications of the states. Since early years of practical implementations of the HF method in nuclear structure physics the simple method to remedy this situation has been devised. Namely, only the fraction f of the mean fields calculated from the single-particlestatesisusedineveryiteration,anditiscombinedwiththe fraction 1 f of the old mean fields. − In practicalcalculations,the optimal value off is determinedby checking the convergence rate. An example of such an analysis is presented in Fig. 1, where the stability energy (5) is shown for the values of f ranging from 30% to 90%. It can be seen that the rate of convergence increases very fast when f increases from 30% to about 80%, but then it decreases again dramatically. Unfortunately, the optimal value of f depends critically on the shell structure of the nucleus; it is rather high for magic nuclei (as superdeformed 152Dy ) 66 86 which have large shell gaps, and has to be lowered when there are several single-particlestatesneartheFermisurface. Inpracticalcalculationsitismore efficienttouseratherlowvaluesoff,andensureconvergenceindependentlyof the shell structure, than to risk a divergence. Usually, f=50% can be a safely recommended value. The total CPU times required to solve the HF equations depend not only on the convergence rate governed by f, but also (of course) on the overall speedofthe computer and(veryimportantly)onthe performancesofthe spe- cific compiler optimization procedures. This is illustrated in Fig. 2, where the relative times required to perform specific tasks are shown in percent of the totaltime. ThetotalCPUtimesrequiredtocalculatethecompletesuperdefor- medbandin152Dy rangefromone houronthe CrayC-90supercomputerto 66 86 several hours on fast workstations (SG Power Challenge L or IBM RS/6000), and to about 36 hours ona typicalworkstation(SUN-Hyper). However,these times result from averaging very different performances pertaining to specific tasks. Typically, three tasks take most of the time, namely, the calculation of 5 102 30% ---> 25 iter/decade 40% 19 ) V 50% 15 Me 100 60% 12 y ( 70% 10 g 80% 9 er 10-2 90% 23 n E f o y 10-4 bilit a St 10-6 10-8 0 20 40 60 80 100 120 140 160 Number of iterations Figure1: Stabilityenergyasfunctionofnumberofiterationsforthecalculationsperformed for the superdeformed state of 16562Dy86. Different curves correspond to different values of thefractionf (seetext), andleadtodifferentnumbers ofiterations whicharenecessaryto improvethestabilitybyanorderofmagnitude. the HO matrix elements of the mean fields, the diagonalizationof the Routhi- ans, and the calculation of the Coulomb field. On a vector computer, the first of these tasks takes relatively long time (40%) because it is composed of mul- tiple short loops,18 whereas the last one contributes very little (5%) due to rather long loops it requires. On superscalarcomputers these three tasks take almostthe sametimesofabout20%each. The exampleoftypicalworkstation shown in Fig. 2 illustrates the importance of the optimization options; here the compiler failed to optimize the subroutine calculating the mean fields and therefore this task takes now almost 30% of the CPU time. 4 Superdeformed bands in the A 60 nuclei ∼ Very recently there have been several reports of superdeformed rotational bandsdiscoveredinnucleiaroundA 60.23,24,25,26 Specific nucleiinthisregion ∼ havealreadybeenaddressedfromthepointofviewoftheNilsson-Strutinsky27 and Relativistic Mean Field28 descriptions. The single-particle structure of 6 Execution times in % Matrix elements Diagonalisation Coulomb Densities HFODD S K S.p. averages S A T Mean fields CRAY-C90 (x1) I/O SG-PCL (x3) Start & End RISC/6000 (x4) SUN (x36) Other 0 10 20 30 40 Figure2: Percentages oftheCPUtimespentonperformingseparatetasks requiredbythe hfoddcodetosolvetheHFequations. the light superdeformed nuclei is very simple and thus amenable to a system- atic description in terms of a microscopic self-consistent mean-field approach. Using the code hfodd (v1.75),29 hundreds of bands have recently been calcu- lated in light Ni-Ga nuclei.30 Based on these results one can find and analyze the generic features exhibited by rapidly rotating nuclei. In the present com- munication we report on quadrupole polarizabilities characterizing individual single-particle orbitals in this region. Following the method introduced in the A 150 nuclei31 we aim at a de- ∼ scription of the proton quadrupole moments Q in terms of contributions 20 givenby the particle (p) and hole (h) states defined with respect to the magic superdeformed nucleus 60Zn , i.e., 30 30 Q =Q 60Zn + δQp n + δQh n (12) 20 20(cid:2)30 30(cid:3) X 20 p X 20 h p h 7 Deviations from additivity 100 hw w =1MeV 80 s 283 bands d n a 60 b f o o. 40 N 20 0 -0.1 -0.05 0 0.05 0.1 d d Q (b) 20 Figure3: DeviationsofcalculatedHFprotonquadrupolemomentsinA∼60nucleifromthe sumofindividualsingle-particlecontributions (12). Thecore60Zn correspondstoalln =0andalln =0,whiletheothersuperde- 30 30 p h formedstates in this regioncanbe obtained by creatingup to 6 particles with n =1 and/or up to 6 holes with n =1, both for protons and for neutrons.30 p h Hence, the sum in Eq. (12) contains altogether 24 unknown coefficients δQ 20 whicharefitted by the least-squaremethod to reproduceallcalculatedproton quadrupole moments Q , separately at every rotational frequency h¯ω. 20 The quality of describing the proton quadrupole moments by the simple additive contributions is illustrated in Fig. 3 where we show deviations of the calculatedHFvaluesfromthesumgiveninEq.(12). Onecanseethatformost of the bands the additivity hypothesis is correctup to 0.05b. This constitutes aremarkableagreementwiththe valuesofQ whichareoftheorderof2–3b. 20 InFig.4thevaluesofindividualcontributionsδQ areshownasfunctions 20 of the rotational frequency for the 24 orbitals considered. It seen that both proton and neutron orbitals significantly contribute to the proton quadrupole moment. This is so because the direct contributionsof protons are ofthe sim- ilarorderas the polarizationcontributionsexisting forboth kinds ofparticles. The positive-parity orbitals [440]1/2 and [431]3/2, originating from the 1g 9/2 intruder orbital, have a larger influence on the quadrupole moment than the negative-parityorbitals(notethatthe scaleistwiceexpandedinthetoppanel of Fig. 4). Deformations gradually increase when consecutive 1g orbitals 9/2 are occupied. The influence of proton positive-parity orbitals is about twice larger than that of their neutron counterparts. The negative-parity [303]7/2 8 0.6 n[431]3/2+ (p) 0.4 n[431]3/2- 0.2 (n) n[440]1/2+ n[440]1/2- 0 p[431]3/2+ (n) -0.2 p[431]3/2- -0.4 (p) p[440]1/2+ -0.6 p[440]1/2- 0.1 n[310]1/2+ n[303]7/2+ 0 n[321]1/2+ ) b ( n[312]5/2+ 0 -0.1 2 n[310]1/2- Q d -0.2 n[303]7/2- d n[321]1/2- -0.3 n[312]5/2- 0.1 p[310]1/2+ p[303]7/2+ 0 p[321]1/2+ p[312]5/2+ -0.1 p[310]1/2- -0.2 p[303]7/2- p[321]1/2- -0.3 p[312]5/2- 0 0.4 0.8 1.2 1.6 hw w (MeV) Figure4: ContributionsδQ20 ofindividualsingle-particleorbitalstotheprotonquadrupole moments of nuclei in the A∼60 region. The orbitals are denoted by the dominant Nilsson labels. Solid (dashed) lines correspond to particle (hole) states, while the closed and open symbols correspond to positive (r=i) and negative (r=−i) signatures, respectively. Error barsaredeterminedbytheleast-squareprocedure. 9 orbitals,originatingfromthe1f extruderorbital,significantlycontributeto 7/2 thequadrupolemomenttoo. Inthiscasetheneutronandprotoncontributions arealmostequal,i.e.,thedirectprotoncontributionisverysmall. Ontheother hand, the proton contribution from the [312]5/2 orbital is much larger than thatofthe neutroncounterpart;hence this orbitalcontributesmostly through the direct quadrupole moment. Contributionscomingfromthesignature-partnerorbitalsarealmostequal, except from those related to the [431]3/2 and [310]1/2 orbitals where the positive-signature partners contribute significantly more. These two orbitals, and also [303]7/2,exhibit contributions slightly dependent on the angular fre- quency. However, these variations are relatively much smaller (do not exceed 0.2b) than the variation of the core quadrupole moment Q 60Zn . In- 20(cid:2)30 30(cid:3) deed, the latter decreases quite substantially with increasing h¯ω, from 3.28b at h¯ω=0.2MeV to 2.56b at h¯ω=1.6MeV. This is at variance with the results calculated for 152Dy , where the core quadrupole moment is fairly indepen- 66 86 dent of spin. Acknowledgments ThisworkissupportedinpartbythePolishCommitteeforScientificResearch (KBN) and by the computational grants from the Interdisciplinary Center for Mathematical and Computational Modeling (ICM) of the Warsaw University andfromtheInstitut du D´eveloppement et de Ressources en Informatique Sci- entifique (IDRIS) of CNRS, France (Project No. 960333). References 1. P. Ring and P. Schuck, The Nuclear Many-Body Problem (Springer- Verlag, Berlin, 1980). 2. P.QuentinandH.Flocard,Annu. Rev. Nucl. Part. Sci. 28,523(1978). 3. S. ˚Aberg, H. Flocard and W. Nazarewicz, Ann. Rev. Nucl. Part. Sci. 40, 439 (1990). 4. C. Baktash, B. Haas, and W. Nazarewicz, Annu. Rev. Nucl. Part. Sci. 45, 485 (1995). 5. P. Butler and W. Nazarewicz, Rev. Mod. Phys. 68, 349 (1996). 6. D. Vautherin and D.M. Brink, Phys. Rev. C5, 626 (1972). 7. D. Gogny, in: Proc. Int. Conf. on Nuclear Physics, eds. J. De Boer and H. Mang (North-Holland, Amsterdam, 1973). 8. D. Gogny, in Nuclear Self-Consistent Fields, ed. G. Ripka and M. Porneuf (North-Holland, Amsterdam, 1975) p. 333. 10