A New Efficient Method of Hartree-Fock-Bogoliubov Calculation for Weakly Bound Nuclei
aa r X i v : . [ nu c l - t h ] F e b A New Efficient Method for Hartree-Fock-Bogoliubov Calculations of Weakly BoundNuclei
M. Stoitsov
Department of Physics, Graduate School of Science, Kyoto University, Kyoto 606-8502, JapanDepartment of Physics and Astronomy, University of Tennessee, Knoxville, Tennessee 37996, USAPhysics Division, Oak Ridge National Laboratory,P.O. Box 2008, Oak Ridge,Tennessee 37831, USA andInstitute of Nuclear Research and Nuclear Energy,Bulgarian Academy of Sciences, Sofia-1784, Bulgaria
N. Michel and K. Matsuyanagi
Department of Physics, Graduate School of Science, Kyoto University, Kyoto 606-8502, Japan (Dated: November 1, 2018)We propose a new method to solve the Hartree-Fock-Bogoliubov equations for weakly boundnuclei, which works for both spherical and axially deformed cases. In this approach, the quasi-particle wave functions are expanded in a complete set of analytical P¨oschl-Teller-Ginocchio andBessel/Coulomb wave functions. Correct asymptotic properties of the quasiparticle wave func-tions are endowed in the proposed algorithm. Good agreement is obtained with the results of theHartree-Fock-Bogoliubov calculation using box boundary condition for a set of benchmark sphericaland deformed nuclei.
PACS numbers: 21.10.-k,21.30.+y,21.60.JzKeywords:
I. INTRODUCTION
The study of nuclei far from stability is an increasinglyimportant part of contemporary nuclear physics. Thistopic is related to newly created radioactive beams facili-ties, allowing more experiments on nuclei beyond the sta-bility line. The new experimental opportunities on nucleiwith extreme isospin ratio and weak binding bring newphenomena which inevitably require a universal theoret-ical description of nuclear properties for all nuclei. Thecurrent approach to the problem is the nuclear densityfunctional theory which implicitly rely on Hartree-Fock-Bogoliubov (HFB) theory, unique in its ability to spanthe whole nuclear chart.The HFB equations can be solved in coordinate spaceusing box boundary condition [1, 2]. This approach (ab-breviated HFB/Box in this paper) has been used as astandard tool in the description of spherical nuclei [3].Its implementation to systems with deformed equilibriumshapes is much more difficult, however. Different ap-proaches have been developed to deal with this problem,such as the two-basis method [4, 5, 6], the canonical-basis framework [7, 8, 9] and basis-spline techniques incoordinate-space calculations developed for axially sym-metric nuclei [10, 11]. These algorithms are precise, buttime-consuming.Configuration-space HFB diagonalization is a usefulalternative to coordinate-space calculations whereby theHFB solution is expanded in a complete set of single-particle states. In this context, the harmonic oscillator(HO) basis turned out to be particularly useful. Overthe years, many configuration-space HFB codes using theHO basis (abbreviated HFB/HO) have been developed, employing either the Skyrme or the Gogny effective in-teractions [12, 13, 14, 15, 16, 17], or using a relativisticLagrangian [19] in the context of the relativistic Hartree-Bogoliubov theory. In the absence of fast coordinate-space methods to obtain deformed HFB solutions, theconfiguration-space approach has proved to be a very fastand efficient alternative allowing large-scale calculations[17, 18].Close to drip lines, however, the continuum states startplaying an increasingly important role and it becomesnecessary to treat the interplay of both continuum anddeformation effects in an appropriate manner. Unfor-tunately, none of the existing configuration-space HFBtechniques manage to incorporate continuum effects.The goal of the present work is to find an efficient nu-merical scheme to solve HFB equations for spherical andaxially deformed nuclei, which properly takes the contin-uum effects into account. We will denote this problem ascontinuum HFB (CHFB). Aiming at treating sphericaland deformed nuclei on the same footing, we rely on theconfiguration-space HFB approach.The HO basis has important numerical advantages; forexample, the use of the Gauss-Hermite quadrature al-lows fast evaluation of matrix elements. On the otherhand, its Gaussian asymptotics prevents from expand-ing systems with large spatial extension, such as halonuclear states. This problem can be successfully fixedby using the transformed HO basis (THO) [20]. Thelatter transforms the unphysical Gaussian fall-off of HOstates into a more physical exponential decay. NeitherHO nor THO bases, however, are able to provide properdiscretization of the quasiparticle continuum. This hasrepercussions already at the HFB level, for which theHO and THO bases cannot reproduce simultaneously allasymptotic properties of nuclear densities (see Sec. V).While this shortcoming is obvious for the HO basis, it alsoarises for the THO basis because the latter can provideonly one type of asymptotic form, i.e. the one insertedin the scaling function defining the THO wave functions[17]. Hence, the THO basis fails to reproduce asymptoticproperties, as asymptotic behavior is different for respec-tive channels: proton and neutron, normal and pairingdensities, different angles for the deformed case. In fact,differences between calculations using the THO and thecoordinate-space bases have been noticed in pairing prop-erties of nuclei (see Sec. V and Ref. [21]). This indicatesthat THO calculations may not always be fully accurateeven in the nuclear region and necessitate careful checkof obtained results. For the aim of carrying out quasipar-ticle random phase approximation (QRPA) calculationswith the HFB quasiparticle representation, the HO andTHO bases are very likely to be insufficient as they can-not provide accurate quasiparticle wave functions in thecontinuum region.Obviously, a more practical basis is needed. TheGamow Hartree-Fock (GHF) basis [22] would be appro-priate, as it has been demonstrated that it can providethe correct asymptotic of loosely bound nuclear states.However, it implies the use of complex symmetric matri-ces. Moreover, the presence of basis states which increaseexponentially in modulus leads to numerical divergences,unless the costly two-basis method is employed [31].As we plan to consider bound HFB ground states only,it is more advantageous numerically to employ Hermitiancompleteness relations, whose radial wave functions arereal. They are either bound, thus integrable, or oscillatewith almost constant amplitude, so that we are free fromthe numerical cancellation problems associated with theGamow states. It should be stressed that we can generatea Gamow quasiparticle basis using the HFB potentialsthus obtained. We can then describe resonant excitedstates by means of the quasiparticle random phase ap-proximation representing the QRPA matrix elements interms of the Gamow quasiparticle basis. This serves asan interesting subject for future investigation.One could expect that the employment of the sphericalHartree-Fock (HF) potential to generate the real contin-uum HF (CHF) complete basis would solve the problem.Unfortunately, the CHF basis is not numerically stabledue to the presence of resonances in the vicinity of thereal continuum. The continuum states lying close to anarrow resonance are rapidly changing, so that a verydense continuum discretization around this resonance isnecessary to accurately represent this energy region. Im-portant numerical cancellations would occur as contin-uum wave functions become very large in amplitude closeto narrow resonances.To overcome this difficulty, we adopt a technique basedon the exactly solvable P¨oschl-Teller-Ginocchio (PTG)potential [23]. The spherical HF potential, seeminglythe best candidate to generate a rapidly converging basis expansion, but providing numerically costly GHF basesor unstable CHF bases, is replaced by a PTG potentialfitted to the HF potential if the latter give rise to resonantstructure. It will be shown that the narrow resonantstates of the GHF basis will become bound in the PTGbasis, so that its scattering states will have no rapid phaseshift change, a necessary condition for numerically stablecontinuum discretization. As a result, we obtain a verygood basis for HFB calculations. We call this approachHFB/PTG.To test the feasibility of this new method, we have per-formed numerical calculations for spherical Ni isotopesnear the drip line, Ni – Ni, for a strongly deformednucleus
Zr, and two HFB solutions for Mg with dif-ferent, prolate and oblate, deformations. Good agree-ment with THO calculations is obtained.The paper is organized as it follows. The HFB/PTGalgorithm is described in Sec. II, while the method usedto generate the PTG basis is formulated in Sec. III.Asymptotic properties of the HFB quasiparticle wavefunctions are discussed in Sec. IV. Results of numeri-cal calculation are presented in Sec. V. Brief summaryand conclusions are given in Sec. VI. Some technical de-tails related to the PTG basis and calculation of matrixelements are collected in Appendices.
II. THE HFB/PTG APPROACH
Our aim is to develop an efficient method of solvingthe CHFB equation Z d r ′ X σ ′ (cid:18) h ( r σ, r ′ σ ′ ) − λ ˜ h ( r σ, r ′ σ ′ )˜ h ( r σ, r ′ σ ′ ) − h ( r σ, r ′ σ ′ ) + λ (cid:19) × (cid:18) U ( E, r ′ σ ′ ) V ( E, r ′ σ ′ ) (cid:19) = E (cid:18) U ( E, r σ ) V ( E, r σ ) (cid:19) (1)for weakly bound nuclei, which equally works both forspherical and axially deformed nuclei. In the above equa-tion, r and σ are the coordinate of the particle in nor-mal and spin space, h ( r σ, r ′ σ ′ ) and ˜ h ( r σ, r ′ σ ′ ) denote theparticle-hole and the particle-particle (hole-hole) com-ponents of the single-particle Hamiltonian, respectively, U ( r σ ) and V ( r σ ) the upper and the lower components ofthe single-quasiparticle wave function, and λ is the chem-ical potential [3]. For simplicity of notation, the isospinindex q is omitted in Eq. (1), but, of course, we solve theCHFB equation for coupled systems of protons and neu-trons. In this section, we outline the calculational schemeand details will be presented in the succeeding sections.The proposed method to solve the CHFB equations,abbreviated HFB/PTG, consists of the following steps:(1) One starts with spherical or deformed HFB cal-culations in the HO basis (HFB/HO). This provides agood approximate solution for the HF potential and theeffective mass.(2) One considers a HF potential and an effectivemass for each ℓj subspace, and fits the associated shiftedPTG potential to them when the HF potential possessesbound or narrow resonant states in this ℓj subspace (seeSec. III A). If no such states appear in the HF ℓj spec-trum, a set of Bessel/Coulomb wave functions [24] is se-lected for the ℓj partial wave basis.(3) One diagonalizes the HFB eigenvalue equations inthe basis composed of the PTG and Bessel/Coulombwave functions. This step continues until self-consistencyis achieved.The use of the Bessel/Coulomb wave functions in step(2) occurs for partial waves of high angular momentum,for which the centrifugal part becomes dominant. Asno resonant structure can appear therein in the real HFcontinuum, Bessel/Coulomb wave functions provide a nu-merically stable set of states for this partial wave. Forthe generation of Coulomb wave functions, one can usethe recently published C++ code [25] or its FORTRANalternative [26]. A complete set of wave functions is thusformed, which will be used as a basis to expand the HFBquasiparticle wave functions.The necessary truncation of the basis in step (3) im-plies that spurious effects may eventually appear at verylarge distances, where both the particle density ρ and thepairing density ˜ ρ are very small. Consequently, quasipar-ticle wave functions have to be matched to their exactasymptotics at moderate distances as it is explained fur-ther in Sec. IV. In addition, special care must be takento calculate matrix elements due to the presence of non-integrable scattering states (see Appendix B).When the HF mean-field resulting from the HFB/HOcalculation in step (1) is deformed, there are several waysto extract the HF potential for each ℓj subspace to beused in step (2). Because it is used just as a gener-ator for the complete PTG basis, its choice will havelittle effect on the final HFB solution, however. In thepresent calculation, we therefore adopt a simple proce-dure; the particle-hole part of the HFB/HO potentialand the HFB/HO effective mass are used in step (2) af-ter averaging their angular and spin degrees of freedom.The resulting HF potential is spherical and the same forall ℓj subspaces. In such a case, the effect of the spin-orbit splitting is not taken into account in the stage ofconstructing the PTG basis but it is of course taken intoaccount in step (3). This implies to consider a basis gen-erated by a spherical potential, which might seem ineffi-cient in the case of large deformation, for which deformedbases are more appropriate, as is done with the HO andTHO bases. The deformed nuclei considered in this paperare nevertheless fairly reproduced within this framework(see Sec. V). If necessary, it is possible to generate a de-formed basis by diagonalizing the deformed HF potentialwithin the PTG basis, which can then serve as a particlebasis for the HFB problem. III. GENERATION OF BASISA. PTG potentials fitting procedure
The PTG potential has four parameters Λ , s, ν and a ,which have to be determined in each ℓj subspace (seeAppendix A). For this purpose, we use the spherical HFpotential and effective mass in a given ℓj subspace. r (fm) -80-60-40-20020 V ( M e V ) neutron FIG. 1: (color online) The shifted PTG potential, the HF po-tential calculated with the SLy4-force and the unshifted PTGpotential for neutrons in Ni. HF and shited PTG potentialsto which centrifugal part is added are provided as well, andthe energies of 0 g / levels for each potential are indicated.All data respectively associated to HF, shifted and unshiftedPTG potentials are respectively shown in solid, dashed anddotted lines. The PTG mass parameter a is obtained from the re-quirement that the PTG and the HF effective masses arethe same at the origin. One first adds the centrifugalterm V ℓ ( ℓ +1) ∝ ℓ ( ℓ + 1) /r to the nuclear plus Coulombpotential, V N + V C , and determines the height E b of thecentrifugal (plus Coulomb) barrier. Then, one adds E b tothe PTG potential; the resulting potential may be calledthe shifted PTG potential. The parameters Λ and ν arefitted in such a way that the χ difference between theshifted PTG potential and the HF potential is minimal.Note that s is directly obtained from Λ and ν values dur-ing the fit, as it is determined by way of the property thatthe PTG potential of parameters Λ , s, ν and a for r → s times the PTG potential of parame-ters Λ , s = 1 , ν and a . The reason why we use the barrierheight E b in our fitting procedure will become apparentby an illustrative example presented below.To test the fitting procedure and the quality of the re-sulting PTG basis we performed GHF calculations in thecoordinate space for the spherical nucleus Ni. Let usexamine the quality of single-particle energies and wavefunctions resulting from the shifted PTG potential bycomparing them with the GHF energies and wave func-tions for bound and resonance states.Figure 1 illustrates the PTG fitting procedure andcompare the results with the GHF ones taking the neu-tron 0 g / level as an example. It is seen that the energyof the bound 0 g / state in the original (unshifted) PTGpotential (horizontal dotted line) become positive afterbeing shifted with E b (horizontal dashed line) and its po-sition agrees in a good approximation with the resonanceenergy obtained by the GHF calculation (horizontal solidline). This is due to a special feature of the PTG poten-tial, for which the centrifugal potential decreases expo-nentially and not as r − for r → + ∞ (see App. A). Thisimplies that the centrifugal + shifted PTG potential goesvery quickly to the constant value, E b , for r → + ∞ .In this way, the PTG treatment replaces the GHF res-onance with a weakly bound PTG state whose wave func-tion will be very similar in the nuclear region. Approx-imating resonant states by weakly bound states in ourframework resembles the standard two-potential methoddescribed in Ref. [27]. Thus, one can expect that thefitted PTG potential provides a rapidly converging basisfor solving the HFB equations.In fact, it is not necessary to find the PTG poten-tial that exactly minimize the χ difference with the HFpotential. As the PTG potential is used as a basis gen-erator, slight differences with the exact minimum leadonly to slightly different bases states to expand the HFBquasiparticle wave functions, preserving its rapidly con-verging properties. Thus, one can take rather large stepsfor the Λ , ν variations and few radii for the χ evaluationto save computer time, keeping the quality of the basisessentially the same. B. Single-particle energies
Single-particle energies and widths for neutrons in Niobtained by the GHF calculations are compared with thePTG energies in Table I. One can clearly see the follow-ing facts.Firstly, the overall agreement between the GHF andthe shifted PTG energies is good, which means that thePTG potential is flexible enough to reproduce the mainfeatures of the HF potential.Secondly, all narrow GHF resonances are representedas weakly bound PTG states with upward shifted PTGenergies. This is very important because the HFB upper(lower) components of quasiparticle states are likely tohave large overlaps with unoccupied (occupied) weaklybound and narrow resonance states.We note that the GHF states whose width is largerthan 1 MeV, as a rule, are not converted to bound PTGstates. This is not important, however, because scat-tering states do not exhibit rapid changes in the energyregion of broad resonances. The broad resonance regioncan indeed be well represented in terms of the continuumbasis states.
TABLE I: Neutron GHF levels in Ni calculated with theSLy4 Skyrme-force and the surface-type delta pairing interac-tion (see Sec. V for the parameters used), which are comparedwith the PTG estimates. All energies are given in MeV whilethe width Γ is given in keV.GHF PTGstates Γ e e + E b e s / s / s / p / p / p / p / d / d / d / d / f / f / g / g / h / C. PTG wave functions -0.4-0.200.20.40.6-0.4-0.200.20.40.6 u (r) proton1f proton0h proton 1d neutron0g neutron0h neutron Γ = 394 keV Γ = 10 keV Γ = 31.6 keV Γ = 26.4 keV Γ = 24.1 keV Γ = 92.9 keV r (fm) FIG. 2: (color online) The PTG (dashed lines) and GHF (solidlines) wave functions for various resonance states.
As illustrated in Fig. 2 narrow GHF resonant statesbear large overlaps with their associated PTG boundstates. Hence, the GHF resonant structure present inthe HFB quasiparticle wave functions will be sustainedby the PTG bound states, thus reducing the coupling tothe PTG scattering continuum.An example indicating the quality of the bound single-particle wave functions resulting from the fitting PTG -25 -20 -15 -10 -5 u (r) r (fm) -0.6-0.4-0.200.20.40.6 u (r) FIG. 3: (color online) The PTG (dashed lines), GHF (solidlines) and HO (dotted lines) wave functions including theasymptotic region for the bound 0 s / , s / and 2 s / neu-tron states both in normal scale (lower panel) and logarithmicscale (upper panel). procedure is shown in Fig. 3 for the bound 0 s / , s / and 2 s / neutron states. In this case, nuclear potentialhas no centrifugal barrier, so that the PTG and the HFpotentials possess the same asymptotic behavior. Verygood agreement between the PTG (dashed lines) and theGHF (solid lines) wave functions is thus not surprising.The upper panel in Fig. 3 shows the asymptotic region inlogarithmic scale where HO wave functions (dotted lines)are also given as a reference. Their Gaussian asymptoticscannot reproduce even approximately the exponential de-crease of the PTG and GHF wave functions. r (fm) -1-0.500.51 u (r) FIG. 4: (color online) The PTG (dashed lines) and GHF (solidlines) wave functions of the neutron continuum s -states cal-culated with energies of 0.118 MeV, 9.996 MeV and 66.119MeV. Neutron continuum s -states are illustrated in Fig. 4, which are properly reproduced as well by the scatteringstates for the PTG potential. In the cases when a cen- r (fm) -1-0.500.51 u (r) FIG. 5: (color online) The PTG (dashed lines) and the GHF(solid lines) wave functions for the neutron continuum d / -states calculated at the same energies as in Fig. 4 trifugal (and/or Coulomb) barrier exists, as illustrated inFig. 5 for d / states, different phase shifts develop in thePTG and GHF continuum states, as the PTG potentialbears no barrier at large distance. IV. QUASIPARTICLE WAVE FUNCTIONS INTHE ASYMPTOTIC REGION
The necessary truncation of the basis implies that spu-rious effects will eventually appear at very large radius,where both the particle density ρ and the pairing density˜ ρ are very small. Consequently, quasiparticle wave func-tions have to be matched with their exact asymptoticsat moderate distance, where the asymptotic region hasbeen attained and densities are still large enough for ba-sis expansion to be precise. Below we explain how thematching procedure is done for axially deformed nuclei.In order to deal with the asymptotics of quasiparticlewave functions, we make partial wave decomposition ofthem: U km ( r σ ) = X α U αkm Ψ α ( r ) = X ℓj U ( ℓj ) km ( r ) Y ℓjkm (Ω) ,V km ( r σ ) = X α V αkm Ψ α ( r ) = X ℓj V ( ℓj ) km ( r ) Y ℓjkm (Ω) , (2)where the subscript k specifies the quasiparticle eigen-states together with the magnetic quantum number m which is always conserved for both spherical and axiallysymmetric nuclei; Ψ α ( r ) are the PTG or Bessel/Coulombwave functions; U αkm and V αkm are the HFB expansion co-efficients; U ( ℓj ) km ( r ) and V ( ℓj ) km ( r ) are the radial amplitudeswith r = | r | for the ( ℓj ) partial wave; Y ℓjm (Ω) denotesa product wave function where the spherical harmonicswith the angular variables Ω and the orbital angular mo-mentum ℓ is coupled with spin to the total angular mo-mentum j .The partial wave amplitudes, U ( ℓj ) km ( r ) and V ( ℓj ) km ( r ), de-fined above involve summation over all quantum numbersexcept the angular momenta ℓ and j . In the sphericalcase, the sums reduce to a single element as ℓ and j aregood quantum numbers. In the asymptotic region, onlyCoulomb and centrifugal parts remain from the HFB po-tentials, so that one can continue the quasiparticle wavefunctions via their partial wave decompositions and de-cay constants k u and k v : U ( ℓj ) km ( r ) = C ( ℓj )+ km H + ℓ,η u ( k u r ) + C ( ℓj ) − km H − ℓ,η u ( k u r ) ,V ( ℓj ) km ( r ) = D ( ℓj )+ km H + ℓ,η v ( k v r ) ,k v = r m ~ ( λ − E ) , k u = r m ~ ( λ + E ) , (3)where E denotes the quasiparticle energy, λ the chem-ical potential, H ± ℓ,η the Hankel (or Coulomb) functions, η being the Sommerfeld parameter, and C ( ℓj )+ km , C ( ℓj ) − km and D ( ℓj )+ km are constants to be determined. Matching isperformed using Eq. (2) at a radius R in the asymp-totic region where the basis expansion is precise, so that C ( ℓj )+ km , C ( ℓj ) − km and D ( ℓj )+ km come forward by continuity.The value of R is typically of the order of 10 fm. V. NUMERICAL EXAMPLES
We have made a feasibility test of the HFB/PTGmethod for spherical Ni isotopes close to the neutrondrip line and for deformed neutron-rich nuclei
Zr and Mg. All calculations were done using the SLy4 den-sity functional. For the pairing interaction, we use thesurface-type delta pairing with the strength t ′ = − . for the the density-independent part and t ′ = − . t ′ MeV fm for the density-dependent part witha sharp energy cut-off at 60 MeV in the quasiparticlespace. They have been fitted to reproduce the neutronpairing gap of Sn. These values are consistent withthose given in Ref. [32]; the slight difference is due todifferent cut-off procedures, sharp cut-off in our case andsmooth cut-off in Ref. [32]. Below we discuss the majorfeatures of the result of calculation. We also make a de-tailed comparison between the HFB/PTG and HFB/Boxcalculations in the spherical case.
A. Spherical nuclei
Let us first examine how the result of calculation de-pends on the truncation of the basis. Indeed, the basishas to be truncated at a maximal linear momentum k max , and discretized with N ℓj continuum states per partialwave in the interval [0 : k max ]. Figure 6 shows that the r (fm) ρ n (r) ρ n (r) k max = 2 fm -1 k max = 3 fm -1 k max = 5 fm -1 ~ FIG. 6: (color online) Dependence on k max of the neutrondensity ρ n and the neutron pairing density ˜ ρ n calculated for Ni by the HFB/PTG method. use of values larger than k max = 3 fm − does not changethe results. Accordingly, in calculations for spherical nu-clei, we use k max = 5 fm − and discretize the continuumwith N ℓj = 60 scattering states per partial wave (seeRef. [22] for its justification).Results of the HFB/PTG calculation for a set of bench-mark Ni isotopes close to the neutron drip line are pre-sented in Table II, Fig. 7 and Fig. 8, where results ofthe HFB/Box calculation are also shown for compari-son. The Ni isotopes are spherical with pairing in theneutron channel only. We see immediately a remark-able agreement between the results of the HFB/PTG andHFB/Box calculations. The difference in total energiesis less than 85 keV and the proton rms radii agree almostperfectly, while the neutron ones are slightly different byless than 0 .
003 fm. Similarly good agreement is obtainedfor all other energy counterparts. The good agreementin the ground state characteristics evaluated by the twodifferent approaches is not surprising if one compares thedensity distributions shown in Fig. 7 and Fig. 8. In thesefigures, the neutron and proton densities, ρ n and ρ p , andthe neutron pairing density ˜ ρ n are plotted both in nor-mal (left column) and logarithmic (right column) scales.The agreement is almost perfect in the whole range of r except at the box boundary where the HFB/Box densi-ties vanish due to the boundary conditions (however notseen in Fig 8). This agreement is striking considering thesignificant impact of the continuum for these nuclei andthe fact that the HFB/PTG calculations are neverthelessperformed using the basis expansion method.Special attention has to be paid to the agreement forthe pairing quantities. Interestingly, the pairing gap ∆ n increases as one approaches the drip line, indicating theimportant role of the pairing correlations in the con- TABLE II: Results of the HFB/PTG calculation for ground state characteristics of Ni isotopes close to the neutron drip line,which are compared with results of the HFB/Box calculation. The SLy4 functional and the surface-type delta pairing[20] areused. The rms radii are in fm and all other quantities are in MeV. Proton chemical potential λ p is not provided as pairingcorrelations vanish in the proton space. Ni Ni Ni NiHFB/Box HFB/PTG HFB/Box HFB/PTG HFB/Box HFB/PTG HFB/Box HFB/PTG λ n -1.453 -1.429 -1.037 -1.029 -0.671 -0.661 -0.342 -0.329 r n r p n E pairn -30.70 -30.60 -36.52 -35.92 -41.98 -41.187 -47.158 -46.233 T n T p E son -63.379 -63.177 -61.679 61.707 -59.558 59.681 -56.898 -57.889 E Couldir E Coulexc -10.138 10.136 -10.084 -10.085 -10.033 10.033 -9.980 -9.980 E tot -654.89 -654.914 -656.933 -656.877 -658.167 -658.084 -658.665 -658.608 tinuum. This result is somehow different from that ofRef. [29] obtained by an alternative HFB calculation inthe coordinate space for the same set of nuclei but it is inagreement with the estimates from Ref. [30]. In Fig. 8,the scaling function of the THO basis is calculated withthe method described in Ref. [20], for which the quasi-exact density provided by the HFB/PTG calculationsis used, and 16 THO shells are taken into account foreach partial wave. This implies virtually optimal results,and it has been checked that densities obtained from theHFB/Box and HFB/THO methods are almost identicalup to 20 fm. On the other hand, pairing densities givenby the THO calculations are not exactly the same withthose of the HFB/PTG and HFB/Box calculations, ascan be seen from Fig. 8. While pairing densities calcu-lated with both methods for Ni and Ni are very close,those for Ni and Ni exhibit visible differences, espe-cially for Ni, for which pairing energies differ by about4 MeV. Asymptotic properties of pairing densities calcu-lated with the THO basis are also not well behaved after15-20 fm, where they saturate instead of decreasing expo-nentially. This indicates that THO basis calculations arenot always devoid of inaccuracies, even at the sphericalHFB level.
B. Axially deformed nuclei
In the case of axially deformed nuclei, few HFB/Boxcalculations are available to check the HFB/PTG re-sults. We consider the well-deformed nucleus
Zr (de-formation β ≈ . Mg. We use therein k max = 4 fm − and N ℓj = 30for all partial waves.Table III compares the three approaches with respect to ground state properties of Zr. In general theyyield similar values. The differences seen in Table III arepartially due to different structure of the model spacesadopted and the associated fitting of the pairing strength.
TABLE III: Comparison of ground state properties of
Zrcalculated with the HFB/Box, HFB/PTG and HFB/THOapproaches. The rms radii are in fm, quadrupole momentsare in barn, and all other quantities are in MeV.HFB/Box HFB/PTG HFB/THO Q tot n E pairn -1.53 -3.015 -2.05 r n r p E tot -893.93 -893.952 -893.711 Proton and neutron densities for nuclei
Zr and Mgare displayed in Fig. 9, with comparison with THO re-sults (circles) for
Zr, in normal scale (left column pan-els) and logarithmic scale (right column panels). Associ-ated pairing densities are shown in Fig. 10.While agreement between the PTG and THO densitiesfor
Zr is good in normal scale, we can notice discrep-ancies in asymptotic properties, which are visible fromthe figure in logarithmic scale (see Fig. 9). It is obviousthat all densities calculated with the THO basis eventu-ally follow the common asymptote dictated by the scalingfunction, while they are well reproduced with use of thePTG basis. This comparison also confirms the presenceof deformation effects even in the far asymptotic region.The middle and bottom panels in Figs. 9 and 10 il-lustrate the HFB/PTG normal and pairing densities for ρ (f m - ) ρ n ρ p ρ p ρ n -28 -21 -14 -7 ρ (f m - ) -28 -21 -14 -7 ρ (f m - ) -28 -21 -14 -7 r (fm) ρ (f m - ) r (fm) -35 -28 -21 -14 -7 Ni Ni Ni Ni FIG. 7: (color online) The neutron densities ρ n and proton densities ρ p both in normal (left-hand side) and logarithmic (right-hand side) scales. Results of the HFB/Box calculation are displayed by solid lines, while those of the HFB/PTG calculationsby open circles and dashed lines. The HFB/HO densities are also indicated by dotted lines in the right panels for comparison. two states with different deformations in the drip linenucleus Mg. These states possess pairing correlationsin both neutron and proton channels. The prolate andoblate states lead to asymptotic neutron densities whichare very close, as seen from the middle and bottom rightpanels in Fig. 9.
VI. CONCLUSIONS
We have proposed a new method of the CHFB calcu-lation for spherical and axially deformed nuclei, whichproperly takes the continuum into account. The methodcombines configuration-space diagonalization of the HFBHamiltonian in the complete set of analytical PTG andBessel/Coulomb wave functions with a matching proce-dure in the coordinate space which restores the correctasymptotic properties of the HFB wave functions. ThePTG potential is chosen to fit the nuclear HF potentialand effective mass. The resulting PTG wave functionsare close to the bound and continuum states of the relatedHF potential while the resonance states are substitutedby the bound PTG states with shifted single-particle en-ergies. Partial waves of high angular momentum are verywell represented by Bessel/Coulomb wave functions. The main results of the present work are twofold:First, we have obtained a new scheme (HFB/PTG) tosolve the CHFB equations as a promising tool for largescale calculation; its performance is comparable, some-times even better, to that of the HFB/THO code, forexample. It properly takes the nuclear continuum intoaccount and therefore could be used for precise densityfunctional calculations for nuclei close to the drip lines.This HFB/PTG method can also be used to provide ac-curate quasiparticle wave functions for microscopic cal-culations of dynamics beyond the nuclear mean-field ap-proximation, as for example, the QRPA calculations fordeformed nuclei.Second, the fact that the HFB/PTG calculation re-produces the results of the coordinate-space HFB calcu-lation with the box boundary condition (HFB/Box) evenfor nuclei up to the neutron drip lines is important. Thisresult indicates the validity of the HFB/Box calculationwhich is widely used, although its validity is sometimesquestioned when it is applied to the drip-line phenomenawhere continuum effects are crucially important [29].The inclusion of resonant structure in the basis is cru-cial for the success of the HFB/PTG approach. Ourtest calculations indicate significant disagreement withthe HFB/Box result if the PTG bound states represent- ρ (f m - ) ~~~~ -12 -8 -4 ρ (f m - ) -12 -8 -4 ρ (f m - ) -12 -8 -4 r (fm) ρ (f m - ) Ni Ni Ni Ni 0 5 10 15 20 25 r (fm) -12 -8 -4 FIG. 8: (color online) The neutron pairing densities ˜ ρ n in normal (left-hand side) and logarithmic (right-hand side) scales.There are no pairing correlations in the proton channel. Results of the HFB/Box and HFB/PTG calculations are displayedboth by solid lines, as they are almost indistinguishable, while HFB/THO pairing densities are represented by dashed lines. ing the resonant GHF states are removed from the basis:in their absence, the pairing densities are overestimatedin the surface region, while particle densities are slightlyunderestimated in the inner region. This means that theresonance states significantly contribute to the total en-ergy through both the particle-hole and particle-particlechannels. Their contributions to the pairing correlationenergy are evaluated to be about 2-3 MeV for the case ofNi isotopes close to the neutron drip line.A more complete investigation of the importance of theHFB resonance states could be made by a detailed com-parison with the result of the exact Gamow-HFB calcu-lation. Such an analysis is in progress for spherical nucleiand will be reported in the near future [31] . Acknowledgments
The authors acknowledge Japan Society for the Pro-motion of Science for awarding The Invitation Fellowshipfor Research in Japan (Long-term) to M. S. and TheJSPS Postdoctoral Fellowship for Foreign Researchersto N. M., which make our collaboration possible. Thiswork was supported by the JSPS Core-to-Core Pro-gram “International Research Network for Exotic Femto Systems.” This work was carried out as a part of theU.S. Department of Energy under Contracts Nos. DE-FG02-96ER40963 (University of Tennessee), DE-AC05-00OR22725 with UT-Battelle, LLC (Oak Ridge NationalLaboratory), and DE-FG05-87ER40361 (Joint Institutefor Heavy Ion Research), the UNEDF SciDAC Collabora-tion supported by the U.S. Department of Energy undergrant No. DE-FC02-07ER41457.
APPENDIX A: PTG BASIS1. PTG potential
The one-body Hamiltonian for the exactly solvablePTG model reads: H P T G = ~ m (cid:18) − ddr µ ( r ) ddr + ℓ ( ℓ + 1) r µ ( r ) (cid:19) + V P T G ( r ) (A1)with m the particle free mass, r is the radial coordinate(in fm), µ ( r ) its dimensionless effective mass (the full ef-fective mass is m µ ( r )), ℓ its orbital angular momentumand V P T G is the PTG potential. The potential V P T G ( r )0 ρ (r) (f m - ) -30 -25 -20 -15 -10 -5 ρ (r) (f m - ) -30 -25 -20 -15 -10 -5 r (fm) ρ (r) (f m - ) r (fm) -30 -25 -20 -15 -10 -5 Zr Mg (prolate) Mg (oblate)
FIG. 9: (color online) The neutron and proton densities of the prolately deformed nucleus
Zr ( β = 0 . Mg calculated by the HFB/PTG method for two states with different deformations (oblate β = − .
09 and prolate β = 0 .
26) in normal (middle and bottom left) and logarithmic (middle and bottom right) scale are alsoprovided with the same line convention. and the effective mass µ ( r ) are written: µ ( r ) = 1 − a (1 − y ) , (A2) V P T G ( r ) = ~ s m µ ( r ) × ( V µ ( r ) + V ℓ ( r ) + V c ( r )) , (A3)where s is the scaling parameter, V µ the potential partissued from the effective mass, V ℓ its ℓ -dependent part and V c its main central part, defined by V µ ( r ) = (cid:2) − a + (cid:2) a (4 − ) − − Λ ) (cid:3) y − (Λ − − a ) + 2 ay ) y (cid:3) × aµ ( r ) (1 − y ) (cid:2) − y (cid:3) , (A4) V ℓ ( r ) = ℓ ( ℓ + 1) (cid:20) (1 − y )(1 + (Λ − y ) y − s r (cid:21) , r > , (A5) V c ( r ) = (1 − y ) (cid:2) − Λ ν ( ν + 1) − Λ − (cid:0) − (7 − Λ ) y − − y (cid:1)(cid:3) . (A6)The quantities V P T G ( r ) and µ ( r ) depend on an implicitfunction y = y ( r ) defined in the following way:Λ s r = arctanh( y ) + √ Λ − √ Λ − y )(A7)1 ρ (r) (f m - ) -7 -6 -5 -4 -3 -2 -1 ρ (r) (f m - ) -7 -6 -5 -4 -3 -2 r (fm) ρ (r) (f m - ) r (fm) -7 -6 -5 -4 -3 -2 ~~~ Zr Mg Mg(oblate)(prolate)
FIG. 10: (color online) Same as in Fig. 9 but for pairing densities and without HFB/THO results. Proton pairing density isnot represented for
Zr as it is negligible therein. so that 0 ≤ y < ≤ r < ∞ .The numerical solution of Eq. (A7) by way of New-ton/bisection methods is stable but one should take spe-cial care at large distances when y becomes closely equalto one. For example, this can be done by rewritingEq. (A7), introducing the new variable x = arctanh( y ):Λ s r = x + √ Λ − √ Λ − x )) , (A8)It is solved with respect to x with a fixed-point algorithm.In this region, 1 − y should be calculated in terms of theexpression 1 − y = 4 e − x / (1 + e − x ) to avoid numeri-cal cancellations.One has to mention that, in the calculation of V P T G ( r ), V ℓ ( r ) is finite for all r ≥ r →
0. Thus, to be precise in thisregion, Eq. (A7) must be rewritten as a power series in y , so that the main diverging terms of Eq. (A5) cancelanalytically.As seen from the equations above, there are four pa-rameters in the PTG model: the effective mass parameter a , the scaling parameter s , the parameter Λ determiningthe shape of the potential and the parameter ν associatedwith the depth of the potential. They can take differentvalues for different angular momenta ℓ . We can use thisfreedom in order to approximate the nuclear potential for each ℓj -subspace, as described in Sec. II.
2. PTG states
The PTG wave functions and eigen-energies are deter-mined by the Schr¨odinger equation for the Hamiltonian(A1) H P T G Ψ k ( r ) = E Ψ k ( r ) (A9)with energies E = ~ k m , (A10)where k stands for the complex linear momentum asso-ciated with E .For bound states, if they exist, the parameter ν de-termines the maximal value n max of the radial quantumnumber n = 0 , , , ..., n max as the largest integer inferiorto (cid:26) (cid:18) ν − ℓ − (cid:19)(cid:27) , (A11)and defines the complex momentum k nl = is − A nl + √ ∆ nl − a , (A12)2with A nl = 2 n + ℓ + 32 , (A13)∆ nl = Λ (cid:18) ν + 12 (cid:19) (1 − a ) − (cid:2) (1 − a )Λ − (cid:3) A nl . (A14)For continuum states, k can take any real positive valuesfrom zero to infinity.
3. PTG wave functions
In order to express the PTG wave function Φ k ( r ) = r Ψ k ( r ) in a closed analytical form, let us introduce thefollowing three functions f k ( r ) = F (cid:18) ν − , ν + , ℓ + 32 , x − (cid:19) (cid:0) x + (cid:1) ¯ β/ , (A15) f + k ( r ) = F (cid:0) ν − , ν + , ¯ β + 1 , x + (cid:1) (cid:0) x + (cid:1) ¯ β/ , (A16) f − k ( r ) = F (cid:0) µ − , µ + , − ¯ β + 1 , x + (cid:1) (cid:0) x + (cid:1) − ¯ β/ (A17)and χ k ( r ) = s x − + Λ (1 − a ) x + √ x − + Λ x + ( x − ) ℓ + 322 , (A18)where x = 1 − (Λ + 1) y − y , x − = 1 − x , x + = 1 + x , (A19) ν + = ℓ + + ¯ β + ¯ ν , ν − = ℓ + + ¯ β − ¯ ν , (A20) µ + = ℓ + − ¯ β + ¯ ν , µ − = ℓ + − ¯ β − ¯ ν , (A21)¯ β = − ik Λ s , (A22)¯ ν = q ( ν + 1 / + ¯ β (1 − Λ (1 − a )) , (A23)and F ( a, b, c, z ) is the Gauss hypergeometric function[24].In the case of bound states, k nl determines the mo-menta k which are pure imaginary (see Eq. (A12)), whilethey are real positive numbers in the case of scattering states. This defines all other quantities entering the equa-tions above. For both cases, the PTG wave functions canbe written either asΦ k ( r ) = N χ k ( r ) f k ( r ) (A24)or as Φ k ( r ) = N χ k ( r ) (cid:0) A + f + k ( r ) + A − f − k ( r ) (cid:1) . (A25)Equation (A24) is suitable for numerical work for smalldistances since x − → r → x − = 1. Similarly, Eq. (A25) is applicable for largedistances since x + → r → + ∞ and the pole x + = 1 of the hypergeometric function in Eqs. (A16)and (A17) is avoided.In the case of bound states, the quantum numbers { nℓ } are the principal quantum number n and the angular mo-mentum ℓ . The constants N , A + , A − entering Eqs. (A24)and (A25) are given by: N = s s ¯ β ( ℓ + + ¯ β + 2 n )( ℓ + + ¯ β Λ (1 − a ) + 2 n ) × s Γ( ℓ + + ¯ β + n )Γ( ℓ + + n )Γ( n + 1)Γ( ¯ β + n + 1)Γ( ℓ + ) ,A + = Γ( ℓ + )Γ( − ¯ β )Γ( µ + )Γ( µ − ) , A − = 0 , (A26)where Γ( z ) is the Gamma function [24].In the case of scattering states, the quantum numbers { kℓ } include the momentum k and the angular momen-tum ℓ while the associated constants N , A + , A − read: N = s Γ( ν + )Γ( ν − )Γ( µ + )Γ( µ − )2 π Γ( ¯ β )Γ( − ¯ β )Γ( ℓ + ) A + = Γ( ℓ + )Γ( − ¯ β )Γ( µ + )Γ( µ − ) A − = Γ( ℓ + )Γ( ¯ β )Γ( ν + )Γ( ν − ) . (A27)The normalization constant N is determined from thenormalization condition Z ∞ Φ nl ( r )Φ n ′ l ( r ) dr = δ nn ′ (A28)for bound states and from the Dirac delta function nor-malization for scattering states: Z ∞ Φ kl ( r )Φ k ′ l ( r ) dr = δ ( k − k ′ ) (A29)All bound and scattering wave functions are orthogo-nal to each other Z ∞ Φ k ( r )Φ k ′ ( r ) dr = 0 , k = k ′ (A30)3and they form a complete basis X nl Φ nl ( r )Φ nl ( r ′ )+ X l Z ∞ Φ kl ( r )Φ kl ( r ′ ) dk = δ ( r − r ′ ) . (A31)One can check that at large distances x → − e − s ( r − r ) , r → + ∞ , (A32)whereΛ s r = p Λ − p Λ − − log (cid:18) Λ2 (cid:19) . (A33)Substituting this into Eq. (A25) one obtains the asymp-totic form of the PTG wave functionsΦ k ( r ) C + e ikr + C − e − ikr (A34)where C + = N A + e − ikr and C − = N A − e ikr , definedby Eqs. (A26) and (A27).The PTG wave functions are numerically stable andaccurate when using Eq. (A24) up to y ≤ .
99 thenapplying the form (A25). They accurately land ontotheir asymptotic representation of Eq. (A34) at large dis-tances.
APPENDIX B: MATRIX ELEMENTS
Let us deal with numerical integration in r and k space.The integration in the r space is performed in terms of N r Gauss-Legendre integration points x i and weights w i within the interval [0 , R max ], Z ∞ O ( r ) Φ k ( r ) Φ k ′ ( r ) dr ≃ N r X i =1 O ( r i ) Φ k ( r i ) Φ k ′ ( r i ) w i , (B1)where O ( r ) is an arbitrary function of r and R max is apoint where nuclear potential disappears. Usually a value R max = 15 fm is used. In the same way, integrationin the k space is done in terms of N k Gauss-Legendreintegration points k i and weights w k i within the interval[0 , k max ], Z k max O ( k ) Φ k ( r ) Φ k ( r ′ ) dk ≃ N k X i =1 O ( k i ) Φ k i ( r ) Φ k i ( r ′ ) w k i , (B2)where O ( k ) is an arbitrary function of k .Radial integrals must be calculated cautiously due tothe presence of non-integrable scattering states in thebasis. When the radial operator represents the nuclear potential or explicitly depends on nuclear densities orcurrents, one can safely integrate the matrix elements tosome large but finite distance R max . Beyond R max , thecontribution of the integral becomes negligible due to thepresence of the densities or currents. However, it is notthe case for the kinetic + Coulomb part of the Hamilto-nian. This operator is infinite-ranged and induces Diracdelta functions in the matrix elements, which have to beregularized directly. For this, one separates the matrix el-ement in two integrals, defined on the intervals [0 : R max ]and [ R max : + ∞ [. The first part is finite and treated withstandard methods. For the second part, if one deals withBessel/Coulomb wave functions, one can assume that thenuclear part is negligible after R max so that they are so-lutions of the asymptotic HF equations. Hence, one ob-tains: Z + ∞ R max u α ( r ) h ( r ) u β ( r ) dr = k α δ αβ − Z R max u α ( r ) u β ( r ) dr ! (bound)= k α δ ( k α − k β ) − Z R max u α ( r ) u β ( r ) dr ! (scat)= − k α Z R max u α ( r ) u β ( r ) dr (mixed) (B3)where h ( r ) is the HF potential which reduces to thekinetic + Coulomb Hamiltonian asymptotically. Here,“bound” (“scat”) means that both α and β states arebound (scattering) and “mixed” means that α is boundand β scattering or vice-versa. The Dirac delta witha discretized basis becomes δ αβ /w k α with w k α beingthe Gauss-Legendre weight associated to the discretizedvalue k α , so that its implementation is immediate; sinceall integrals are finite, they pose no problem. When thePTG basis states are used instead of the Bessel/Coulombwave functions, it turned out that it is numerically preciseto disregard the Coulomb/centrifugal part of the Hamil-tonian after R max , so that Eq. (B3) is the same for boththe PTG and Bessel/Coulomb wave functions. Indeed,Eqs. (A32) and (A34) imply that the PTG wave func-tions behave asymptotically like neutron waves functionsof angular momentum ℓ = 0. The above seemingly crudeapproximation can, in fact, be mathematically justified.The HFB matrix evaluated using such a procedure con-verges weakly to the exact HFB matrix for R max → + ∞ [33]. This means that the HFB matrix elements dependon R max asymptotically, some of them even divergingwith R max → + ∞ , whereas its eigenvalues and eigenvec-tors converge to a finite value.4 [1] A. Bulgac, Preprint FT-194-1980, Central Institute ofPhysics, Bucharest, 1980; nucl-th/9907088.[2] J. Dobaczewski, H. Flocard, and J. Treiner, Nucl. Phys. A422 , 103 (1984).[3] J. Dobaczewski, W. Nazarewicz, T.R. Werner,J.F. Berger, C.R. Chinn and J. Decharg´e, Phys. Rev. C , 2809 (1996).[4] B. Gall, P. Bonche, J. Dobaczewski, H. Flocard and P.-H. Heenen, Z. Phys. A348 , 183 (1994).[5] J. Terasaki, H. Flocard, P.H. Heenen and P. Bonche,Nucl. Phys.
A621 , 706 (1997).[6] M. Yamagami, K. Matsuyanagi and M. Matsuo, Nucl.Phys.
A693 , 579 (2001).[7] P.-G. Reinhard, M. Bender, K. Rut, and J.A. Maruhn,Z. Phys.
A358 , 277 (1997).[8] N. Tajima, RIKEN Review , 29 (1998).[9] N. Tajima, Phys. Rev. C , 034305 (2004).[10] E. Ter´an, V.E. Oberacker and A.S. Umar, Phys. Rev. C , 064314 (2003).[11] V.E. Oberacker, A.S. Umar, E. Ter´an andA. Blazkiewicz, Phys. Rev. C , 064302 (2003).[12] D. Gogny, Nucl. Phys. A237 , 399 (1975).[13] M. Girod and B. Grammaticos, Phys. Rev. C , 2317(1983).[14] J.L. Egido, H.-J. Mang, and P. Ring, Nucl. Phys. A334 ,1 (1980).[15] J.L. Egido, J. Lessing, V. Martin and L.M. Robledo,Nucl. Phys.
A594 , 70 (1995).[16] J. Dobaczewski and P. Olbratowski, Comput. Phys.Commun. , 158 (2004).[17] M.V. Stoitsov, J. Dobaczewski, W. Nazarewicz and P. Ring, Comput. Phys. Commun. 167, 43 (2005).[18] S. Goriely, M. Samyn, P.-H. Heenen, J. M. Pearson andF. Tondeur, Phys. Rev. C , 024326 (2002).[19] P. Ring, Prog. Part. Nucl. Phys. , 193 (1996).[20] M.V. Stoitsov , J. Dobaczewski, W. Nazarewicz, S. Pitteland D.J. Dean, Phys. Rev. C , 054312 (2003).[21] A. Blazkiewicz, V.E. Oberacker, A.S. Umar andM. Stoitsov, Phys. Rev. C , 054321 (2005).[22] N. Michel, W. Nazarewicz and M. Ploszajczak, Phys.Rev. C , 064313 (2004).[23] J. Ginocchio, Ann. of Phys. , 467 (1985).[24] M. Abramowitz and I.A. Stegun, Handbook of Mathemat-ical Functions (Dover, New York, 1970).[25] N. Michel, Comp. Phys. Comm., , 232 (2007).[26] I.J. Thomson and A.R. Barnett, Comp. Phys. Comm., , 363 (1986).[27] S.A. Gurvitz, P. B. Semmes, W. Nazarewicz andT. Vertse, Phys. Rev. A , 042705 (2004).[28] E. Chabanat, P. Bonche, P. Haensel, J. Meyer andF. Schaeffer, Nucl. Phys. A635 , 231 (1998).[29] M. Grasso, N. Sandulescu, Nguen Van Giai and R.J. Li-otta, Phys. Rev. C , 064308 (2005).[31] N. Michel, M. Stoitsov and K. Matsuyanagi, in prepara-tion.[32] K. Bennaceur and J. Dobaczewski, Comp. Phys. Comm.,168