Canonical-basis solution of the Hartree-Fock-Bogoliubov equation on three-dimensional Cartesian mesh

1.1 MB·English

by  Naoki Tajima

Canonical-basis solution of the Hartree-Fock-Bogoliubov equation on
three-dimensional Cartesian mesh
Naoki Tajima∗
Department of Applied Physics, Fukui University
(Dated: February 8, 2008)
A method is presented to obtain the canonical-form solutions of the HFB equation for atomic
nuclei with zero-range interactions like the Skyrme force. It is appropriate to describe pairing
correlationsinthecontinuumincoordinate-spacerepresentations. Animprovedgradientmethodis
used for faster convergences under constraint of orthogonality between orbitals. To prevent high-
lying orbitals to shrink into a spatial point, a repulsive momentum dependent force is introduced,
whichturnsouttounveilthenatureofhigh-lyingcanonical-basisorbitals. Theasymptoticproperties
at large radius and the relation with quasiparticle states are discussed for the obtained canonical
basis.
4
0 PACSnumbers: 21.10.Pc,21.30.Fe,21.60.Jz,27.30.+t
0
2
I. INTRODUCTION takeintoaccounttheseeffectsleadstotheHartree-Fock-
n
a Bogoliubov (HFB) equation, with which the density is
J localized whenever the Fermi levels are negative [3, 4].
Pairingcorrelationsplay anessentialrole inthe deter-

7 The HFB method in the coordinate space was first
mination of the ground-state structure of the vast ma-
2 formulated using the quasiparticle states and solved for
jority of atomic nuclei. Among its treatments, an easy
but still enough accurate one is to consider only single- sphericallysymmetric statesinRef.[3]. Althoughspher-
2
ical solutions can be obtained easily with present com-
v particlestatesneartheFermilevelwhiletheeffectsofthe
puters(forzero-rangeforces),deformedsolutionsarestill
5 other states are assumed to be absorbed in the strength
7 of an effective pairing interaction. In such a treatment, difficult to obtain because there are quite a large num-
0 our experiences in Hartree-Fock(HF)+BCS calculations berofquasiparticlestatesevenforamoderatesizeofthe
7 normalizationbox(i.e.,thecavitytoconfinethenucleons
suggestthatatleasthalfofamajorshellabovetheFermi
0 to discretize the positive-energy single-particle states).
level must be considered for meaningful estimations. As
3 An orthodox approach to face this difficulty is the two-
a consequence, one has to explicitly consider positive-
0
basis method [5, 6, 7, 8, 9], in which the quasiparticle
/ energystatesiftheFermilevelofneutronsishigherthan
h thenegativeofhalfofthemajorshellspacing( 41A−1/3 states are expanded in bound and unbound HF orbitals.
-t MeV).Thisconditionappliestoabouthalfoft∼he2104 nu- This method requires heavy numerical calculations be-
l cause there are a pile of unbound HF orbitals below an
c clides whichexist between protonandneutron driplines
u in the nuclear chart, not only to nuclei near the neu- energy cut-off of even only a few MeV.
n tron drip line or outside the s-process path. Thus there An alternative approach is the canonical-basis HFB
v: areplentyofnecessitiesfortheoreticalframeworkswhich method. According to the Bloch-Messiah theorem [10],
i candescribepairingcorrelationsinvolvingthecontinuum every HFB solution obtained as the vacuum of a set of
X
states. Bogoliubovquasiparticleshasanequivalentexpressionin
r terms of a BCS-type wave function. The single-particle
a Ifthepairingcorrelationsignificantlyinvolvesthecon-
states appearing in this expression are called the HFB
tinuum states, the HF+BCS approximation becomes in-
canonicalbasis. Inthecanonical-basisHFBmethod,one
adequatebecausetheoccupationofunboundHForbitals
obtains the solution in the canonicalform without using
leadstothe unwanteddislocalizationofthe nucleonden-
thequasiparticlestates. Thismethodappearedoriginally
sity. This is a serious problem in coordinate-spacetreat-
in Ref. [11] to obtainsphericalsolutions. However,there
ments,whichismorefavorabletodescribelooselybound
are no sever difficulties in obtaining HFB solutions for
systemslikedrip-linenucleithanexpansionsinharmonic-
sphericalnucleiusinganyothermethods. Applicationsto
oscillator eigenstates (except the transformed oscillator
deformednucleihavebeendonebyus[8,12,13,14]using
basis of Ref. [1]). The reason why a dislocalized solu-
a three-dimensional Cartesian mesh representation[15].
tion becomes the variational minimum is that, in order
Incidentally, a different line of application is also found
to separate the variational equations into HF and BCS,
in literature [18].
onehastoneglectthe effectsthatthe matrixelementsof
Let us explain the advantage of the canonical-basis
pair-scattering processes are affected by the changes in
method over the two-basis method concerning the treat-
the wave functions of the orbitals involved [2]. To fully
ment of the pairing in the continuum using Fig. 1. On
both sides of the figure, the ordinate represents the ex-
pectation value of the mean-field (HF) energy, while the
∗[email protected]; abscissa stands for the radius r from the center of the
http://serv.apphy.fukui-u.ac.jp/~tajima nucleus. ǫc and ǫF mean a cut-off energy and the Fermi
2
ported on the discovery of the point-collapse problem
and the tests of the state-dependent cut-off factors and
pairing-density-dependentforcesastheremedies. Inthis
paper we have rewritten those contents to include the
momentum-dependent interaction and, by including the
rewritten contents, organized the formalism for consis-
tency and selfcontainedness. All the numerical calcula-
FIG. 1: Comparison between the HF orbitals and the HFB tions shownin this paper havebeenperformedusing the
canonical-basis orbitals. See text for explanations. renewed form of the interaction.
Incidentally,someimportantpartsofthosereportsare
notrepeatedinthispaper. Theexamplesareananalysis
level, respectively. On the left-hand side, the wavy hori- of the spherical two-basis method concerning the preci-
zontal lines stand for energy levels and their spatial ex- sion of density tails [8] and numerical demonstrations of
tent for the HF potential denoted by a solid curve. The the point collapses of high-lying orbitals [14].
wavefunctionsofpositiveenergystatesextendtothewall
oftheboxandtheirleveldensityismuchlargerthanthat
ofnegativeenergystates. Inthe two-basisHFB method, II. HFB IN THE CANONICAL
onehastomixthesepositive-energyorbitalstoconstruct REPRESENTATION
localized canonical-basis orbitals, which is a numerically
demanding task. On the right-hand side, the wavy lines To begin with, let us formulate HF and HFB methods
represent the HFB canonical-basis orbitals. Unlike HF in the coordinate-space representation in order to eluci-
orbitals,theyarespatiallylocalizedforbothnegativeand datethedifficultyofthequasiparticle-basedHFBmethod
positive energies. Because of this localization, the level comparedwith the HF method and to propose its possi-
density is much smaller than the unbound HF orbitals. blesolutionintermsofthecanonical-basisHFBmethod.
Therefore, one needs much less orbitals. More specifi- For the sake of simplicity, we consider only one kind of
cally, the number of necessary single-particle states to nucleons in this section. The number of that kind of nu-
describe the HFB ground state of a nucleus is propor- cleonsaredesignatedbyN. Thez-componentofthespin
tionaltothevolumeofthenucleusinthecanonical-basis of a nucleon is represented by s (= 1).
method while it is related to the volume of the normal- In HF, one should minimize ΨH±2Ψ among all the
ization box in the other methods. Incidentally, with the N-body single Slater-determinanht s|tat|esi,
dash curve, we suggest the existence of some potential
which binds the high-lying canonical-basis orbitals. The N
identity of this potential is unveiled in Sec. VIB. Ψ = a† 0 , (1)
| i i| i
In this paper we formulate the canonical-basis HFB i=1
Y
methodonacubicmeshanddevelopanefficientgradient
a† = d3r ψ (r,s)a†(r,s), (2)
method to obtain its solutions. In order to decrease un- i i
s Z
physical influences of high-momentum components due X
tozero-rangeinteractions,weintroduceamomentumde-
by varying ψ under orthonormalityconditions
i i=1,···,N
pendent term to the pairing interactionand show how it ψ ψ = δ{ .}The operator a† creates a nucleon with
suppresses a problematic behavior of wave functions pe- ha wi|avjei funcitjion ψ (r,s). Theiket state 0 stands for
culiar to the canonical-basisHFB method. As a byprod- i | i
the state in which no nucleons exist. The distribution
uct, we show a clear way to understand the nature of
function of the number density of the nucleons is related
high-lying canonical-basis orbitals.
to the single-particle wave functions as
The contents of this paper is as follows. In sections
II – V, we give the formulation and general considera- N
tions. Discussionsusing resultsofnumericalcalculations ρ(r)= ψ (r,s)2. (3)
i
are collected in section VI. Conclusions are given in sec- | |
i=1 s
XX
tion VII. In Appendix A we examine the errors due to
themeshrepresentationandtheboxboundarycondition. There is arbitrariness in the selection of ψi as for uni-
{ }
Appendix B gives a model analysis of the problempecu- tary transformations among them. ψi are called HF or-
liar to the method, i.e., a phenomenon that canonical- bitals when they diagonalize the HF Hamiltonian.
basis orbitals collapse to a spatial point once their occu- In HFB, the solution takes the following form,
pation probabilities fall below some critical value when
pure delta-function forces are used. NB
Ψ = b 0 , (4)
Some of the contents of this paper have already been | i i| i
discussed in our earlier reports [8, 12, 13, 14]. For ex- iY=1
ample, Ref.[8]includes the firstreportonthe canonical- b = d3r φ∗(r,s)a(r,s)+ϕ∗(r,s)a†(r,s)(5,)
i i i
basis HFB method on a cubic mesh. Refs. [13, 14] re-
s Z
X (cid:8) (cid:9)
3
where b is the annihilationoperatorofapositive-energy a condition on the expectation value of the number of
i
Bogoliubov quasiparticle with two types of amplitudes nucleons,
φ (r,s) and ϕ (r,s) for particle-like and hole-like parts
i i
of the excitation. NB is the number of the basis states Nw
of the representation used and is by far larger than N. 2 vi2 =N, (13)
Oneshouldvary φ ,ϕ underanorthonormal- i=1
{ i i}i=1,···,NB X
ity condition
and the normalization of the u-v factors, u2+v2 =1.
i i
φ∗(r,s)φ (r,s)+ϕ∗(r,s)ϕ (r,s) d3r =δ , Reinhard et al. wrote that the advantage of the repre-
{ i j i j } ij sentation (9) over (4) is that one has to consider only a
s Z
X (6) singlesetofwavefunctions ψi ,notadoubleset φi,ϕi
{ } { }
and a constrainton the expectation value of the number [11]. In our opinion, however, one may expect much
of nucleons, greater benefit from the canonical-basis representation.
Namely, i may be truncated as i N = (N) 1N
≤ w O ≪ 2 B
ρ(r)d3r=N, (7) to a very good approximation.
The reasonis the localizationof the density. The den-
Z
where the nucleon density for state (4) is expressed as sitydistributionofHFBsolutionscanbeshowntobehave
asymptotically for large r as
NB
ρ(r)= |ϕi(r,s)|2. (8) e−κr 2 √ 2mǫF
i=1 s ρ(r) , κ − , (14)
XX ∼ r ≥ ¯h
The essential difference between HF and HFB lies in (cid:18) (cid:19)
the number of necessary single-particle wave functions. asfarastheFermilevelǫ isnegative[3],wheremisthe
F
It is N in the former and 2NB in the latter. Obviously, nucleon mass. Consequently, ψi appearing on the right-
the latter is much larger. handsideofEq.(11)mustbealocalizedfunctionasρ(r)
OwingtotheBloch-Messiahtheorem[10],thestate(4) ontheleft-handside. Ontheotherhand,theorthogonal-
can be transformed(as for the ordinaryground states of ity relation (12) restricts the number of wave functions
even-even nuclei) into the following form, which can exist in the vicinity of the nucleus. Therefore
thenumberofcanonicalorbitalscannotbeverylarge. In
Nw
Ψ = u +v a† a† 0 , (9) mesh representations, NB is proportional to the volume
| i i i i ¯ı | i of the box while N is proportionalto the volume of the
iY=1(cid:16) (cid:17) nucleus. The latterwis 10−1 10−2 times as small as the
a† = d3r ψ (r,s)a†(r,s), (10) former in typical calculation−s.
i(¯ı) i(¯ı)
s Z Situation is quite different in the quasiparticle formal-
X
where a† and a† create a nucleon with wave functions ism: The number of quasiparticle states is proportional
i ¯ı tothevolumeofthebox. Thisisbecausethelocalization
ψ (r,s) and ψ (r,s), respectively, which are called the
i ¯ı of the density demands only that of ϕ through Eq. (8)
canonical basis of the HFB vacuum Ψ . They may i
be called the natural orbitals[11] instea|d,ibecause they while φi does not have to be localized in general. The
orthogonalitycondition(6)involvesbothϕ andφ . This
are the eigenstates of the one-body density matrix. In i i
enablesmanyquasiparticlestateshavingsimilarϕ tobe
this paper, we often callthem canonical-basisorbitalsor i
orthogonalto each other by differing their φ .
canonicalorbitals. ItisrequiredN = 1N fortheexact i
w 2 B
equivalence between Eqs. (4) and (9) to hold in general.
Inthispaperweassumethat Ψ isinvariantunderthe
| i III. MEAN FIELDS FOR ZERO-RANGE
time reversal operation. In this case, ψ and ψ are the
i ¯ı INTERACTIONS
time-reversalstatetoeachotherandonlyoneofthemhas
to be considered explicitly in the variational procedure.
Inthispaper,weemploytheSkyrmeinteraction[2,16],
The nucleon density can be expressed as
whichisadensity-andmomentum-dependentzero-range
Nw interaction:
ρ(r)=2 v2 ψ (r,s)2, (11)
i | i |
Xi=1Xs vˆ = t0(1+x0Pσ)δ+ 21t1(1+x1Pσ)(k2δ+δk2)
in the time-reversalinvariant case. + t (1+x P )k δk+ 1t (1+x P )ραδ, (15)
In orderto obtainsolutions expressedin the canonical 2 2 σ · 6 3 3 σ
form without using the information of the quasiparticle where P is the spin exchange operator, δ = δ(r r )
σ 1 2
states, one should vary {ψi,ui,vi}i=1,···,Nw under three withr1 andr2 thepositionvectorsoftheinteractin−gtwo
kinds of constraints, i.e, the orthonormality condition, nucleons, ρ is the nucleon density at r (=r ), and
1 2
ψ ψ = ψ∗(r,s)ψ (r,s)d3r =δ , (12) 1
h i| ji i j ij k= ( r r ). (16)
s Z 2i ∇ 1 −∇ 2
X
4
¯hk is the relative momentum between the two nucleons. reason[20].) Thedependenceonthepairingdensityρ˜(to
Since it is hermitian under the box (vanishing and peri- be defined by Eq. (26)) has been introduced in Ref. [13]
odic)boundaryconditions,wehavenotspecifiedinwhich to prevent unphysical behaviors of high-lying canonical
way(tobraorketstates)itoperatesinEq.(15)unlikein, orbitals discussed in Appendix B. This term is squared
e.g.,Ref.[2]. Zero-rangenessmakesthemean-fieldpoten- because ρ˜ can become negative in principle. Except in
tialslocal,whichisanessentialadvantageforcoordinate- Appendix B, we use ρ˜ = . The momentum depen-
c
∞
space solutions. dent term acts on s-wave relative motions and quenches
Among the terms of the Skyrme force, those with interactionsbetweennucleonsin high-momentumstates.
strength parameters t , t , and t act only on s-wave Theinteractionvanishesatrelativemomentumk atlow-
0 1 3 c
relative motions. The term with t serves to take into density points. If ρ˜ = k = , this pairing force is re-
1 c c
∞
account the effects of the finite-rangeness in terms of duced to that introduced in [19]. The typical values of
the lowest order momentum dependence, while the term theparametersarev = 1000MeVfm3 andk =2fm−1
p c
−
with t expresses a density dependence. The term with but different values are also used.
3
strength t2 acts on relative p-wave states through a TheS=0andT=1partoftheSkyrmeforceofEq.(15)
different form of momentum dependence. Note that maybe rewrittenasEq.(19)withoutthepairingdensity
the isospin channel (T=0 or 1) is uniquely determined dependent term. Using different parameterizations be-
through the requirement of the two-body antisymmetry. tween the two forces has an advantage of avoiding con-
Owing to the p-wave as well as the s-wave interaction fusions.
terms, the Skyrme force can act on all the four kinds of
When one considers both protons and neutrons, the
spin-isospin channels.
state of the nucleus is usually assumed to be a product
The completeSkyrmeinteractionincludes aspin-orbit
of two BCS-type states of Eq. (9), one for the protons
termiW(σ1+σ2) k δk,whichwehaveexcludedfrom and the other for the neutrons:
· ×
Eq. (15). This elimination decreases the size of the nu-
merical calculations by a factor of four (two from spin n Nw(q)
up and down components, another two from real and Ψ = u +v a† a† 0 , (20)
imaginary parts of the spatial wave function). Due to | i qi qi qi q¯ı | i
this size reduction, we can perform very severe tests of qY=p iY=1 (cid:16) (cid:17)
the canonical-basisHFB method by, e.g., taking into ac-
wherethe index q distinguishes betweenprotons(p) and
count unusually many canonical orbitals or expressing
neutrons(n). a† createsaprotonhavingawavefunction
wave functions on a very large mesh. pi
As the parameters of the force, we choose the SIII set ψpi(r,s)whilea†nicreatesaneutronwithawavefunction
[17] except for the spin-orbit term, i.e., t0 = 1128.75 ψni(r,s). This product form matches our choice of the
MeVfm3,t =14000MeVfm3+3α,α=1,andW−=0. As pairing interaction of Eq. (19) which acts only on T =1
3
for the parameters appearing in the momentum depen- pairs.
dent terms,the following twocombinations aresufficient In this paper we treat only N=Z nuclei without
to treat N =Z even-even nuclei, Coulomb interaction for the sake of simplicity. In this
case, the wave functions are the same between protons
f = 1 (3t +5t +4t x ), (17)
1 16 1 2 2 2 and neutrons. Moreover, because the potentials are in-
f2 = 614(−9t1+5t2+4t2x2), (18) dependent of the spin, the wave functions ψi(r,s) can
be factorized into a product of a spin wave function and
whosevaluesaref =44.375MeVfm5 andf = 62.969 a real-number function of the position, which we write
1 2
MeV fm5. − ψ (r) hereafter. More specifically, we assume ψ (r,s)
i pi
Asisusuallydone,weemploydifferentparametersbe- = ψ (r,s) = ψ (r)δ and ψ (r,s) = ψ (r,s) =
ni i s,1 p¯ı n¯ı
tween the mean-field and the pairing interactions. We ψ (r)δ . 2
i s,−1
use the following interaction for the pairing channel, 2
The isoscalarpairingmaytakeoverthe isovectorpair-
2 ing in N Z nuclei [21]. Nevertheless, we assume
vˆ = v 1−Pσ 1 ρ ρ˜ δ Eq. (20) b∼ecause the aim of this paper is to examine
p p
2 "( − ρc −(cid:18)ρ˜c(cid:19) ) a method whose most important applications are to the
1 nuclei near the neutron drip line.
k2δ+δk2 , (19)
−2k2 Our Hamiltonian H consists of a kinetic energy term
c (cid:21)
(cid:0) (cid:1) andtwo-bodyinteractionterms. Itshouldbeunderstood
where v is the strength and 1(1 P ) is the projector thatbothareexpressedinthesecondquantizedformbe-
p 2 − σ
into spin-singlet states. Terms in the braces represent causeHFBstatesdonothaveafixednumberofparticles.
dependences on particle and pairing densities. The in- The kinetic energy term is ¯h2 2/2m in the one-body
− ∇
teraction becomes repulsive where ρ>ρ or ρ˜>ρ˜ . For expression, for which we assume the mass of a nucleon
c c
theparticle-densitydependence,weuseρ =0.32fm−3to m to be the average of the proton and the neutron bare
c
vanishthe volume-changingeffect[19](This choice,twice masses(m andm )withanapproximatecorrectionfac-
p n
ofthematterdensity,isalsorecommendedforadifferent tor for the center-of-mass motion of a nucleus of mass
5
number A: where h and h˜ are called the mean-field and pairing
Hamiltonians. They are given by
−1
1 m +m
p n
m= 1 . (21) h = ∇ B∇+V, (30)
− A 2
(cid:18) (cid:19) − ·
h˜ = ∇ B˜∇+V˜, (31)
Forthe interactionterms,we use the interaction(15)for − ·
themean-fieldtypecontractionsandtheinteraction(19) with
for the pairing type contractions in evaluating the ma-
¯h2
trix elements of the two-body interactions using Wick’s B = +f ρ, (32)
1
2m
theorem.
The totalenergyforthe state (20)canbe expressedin V = 3t ρ+ 2+αt ρ1+α+f τ +2f 2ρ vpρ˜2,(33)
terms of a single space integral as, 4 0 16 3 1 2∇ − 8ρc
v
B˜ = p ρ˜, (34)
Etot = ΨH Ψ = (r)d3r, (22) −8kc2
h | | i H
Z v ρ ρ˜ 2 τ˜ 2ρ˜
V˜ = p 1 2 ρ˜ + ∇ .(35)
where the integrand is given by 4"( − ρc − (cid:18)ρ˜c(cid:19) ) − 2kc2 kc2 #
= ¯h2 τ + 3t ρ2+f ρτ +f ρ 2ρ+ 1 t ρ2+α We call V and V˜ the mean-field and the pairing poten-
0 1 2 3
H 2m 8 ∇ 16 tials, respectively.
2 When there are no pairing correlations, the mean-
v ρ ρ˜ ρ˜
+ p 1 ρ˜2 τ˜ 2ρ˜ .(23) field Hamiltonian h is identical to the HF single-particle
8"( −ρc−(cid:18)ρ˜c(cid:19) ) − kc2 −∇ # Hamiltonian. Therefore, its expectation value,
(cid:0) (cid:1)
On the right-hand side of the above equation, functions ǫ = ψ hψ , (36)
i i i
of position r are h | | i
maybecalledthesingle-particleenergyoftheithcanon-
Nw icalorbital. The negativeof the expectationvalue of the
ρ(r) = 4 vi2|ψi(r)|2, (24) pairing Hamiltonian,
i=1
X ∆ = ψ h˜ ψ , (37)
Nw i −h i| | ii
τ(r) = 4 v2 ∇ψ (r)2, (25)
i| i | has a meaning of the pairing gap of the orbital. The
i=1
X pairing gap should be positive or zero. Indeed, ∆ could
i
Nw takenegativevaluesonlyifthewavefunctionwouldcon-
ρ˜(r) = 4 u v ψ (r)2, (26)
i i| i | sist of very high-momentum components [22]. A related
Xi=1 discussion is given in Sec. VIB.
Nw The quasiparticle states of Eq. (5) are the eigenstates
τ˜(r) = 4 u v ∇ψ (r)2, (27)
i i| i | of their own Hamiltonian, i.e., so-called HFB super ma-
Xi=1 trix composed of h and ˜h:
which are called [3] the (particle) density, the kinetic- h ǫ h˜ φ φ
energy density, the pairing density, and the pairing − F i =E(i) i , (38)
kinetic-energy density, respectively. (cid:18) ˜h −h+ǫF (cid:19)(cid:18)ϕi (cid:19) qp (cid:18)ϕi (cid:19)
There is a set of single-particle operators (Hi, h, and where ǫF is the Fermi level. In small subspaces like sev-
h˜) which are useful to obtain mean-field solutions. They eral major shells, the easiest way to obtain an HFB so-
can be obtained by taking the functional derivative of lution is to solve Eq. (38) as an eigenvalue problem of a
Etot with respect to a canonical wave function ψi. This hermitian matrix. In large spaces like mesh representa-
turnsouttobeequivalenttooperatingastate-dependent tions,however,itiseasiertoobtainthecanonicalorbitals
single-particle Hamiltonian i to the wave function : directly, rather than via quasiparticle states. Such a di-
H
rect method is explained in the next section.
δE δE
δψtoit =4Hiψi∗, δψti∗ot =4Hiψi. (28) funInctciiodne,ntoanlley,mifayonneoctoicnesitdhearts stehreiofuusnlycttiohnaatlψdieirsivaarteivael
should be doubled since
Here,wehavetakenthefamiliarpointofviewthatψ and
i
ψi∗ can be formally regarded as independent variables, δEtot = δEtot δψibra + δEtot δψiket =2 4 ψ , (39)
although ψi is actually a real function. δψi δψibra δψi δψiket δψi · Hi i
One can express , in turn, as a linear combination
i
of state-independentHsingle-particle Hamiltonians as, as ψibra = ψi and ψiket = ψi. Eq. (39) seems to be a
more appropriate expression of what our numerical cal-
=v2h+u v h˜, (29) culationprogramactuallydoes. However,thefactor2at
Hi i i i
6
the beginning ofthe right-handside as wellas the factor with
4 which reflects the spin-isospin degeneracy do not mat-
µ2
ter because they can be conveniently absorbed into the f =exp k2 , k2 = ψ∗(r) 2ψ (r)d3r. (43)
step-size parameter of the gradientmethod (see Sec. V). i − 4 i i − i ∇ i
(cid:18) (cid:19) Z
In Eq. (43), we assume a dependence on the kinetic
IV. CUT-OFF SCHEMES FOR THE PAIRING energy (∝ ki2), not on the mean-field energy ǫi as in
INTERACTION Eq.(41), to avoida highly complicatedexpressionof the
gradientfor the lattercase. (In BCS,suchcomplications
are simply neglected.) The gradient is given for k =
In this section we discuss on the cut-off schemes c
∞
by
for zero-range pairing interactions in relation to the
canonical-basis HFB method.
1δE
Delta function forces without a cut-offlead to a diver- tot = v2h+u v h˜ f(k2)+f′(k2) 2 ψ . (44)
4 δψ∗ i i i i i ∇ i
gence of the strength of the pairing correlation. In order i
h (cid:8) (cid:9)i
to prevent the divergence, one usually takes only those
This methodworkswellto preventthe pointcollapsesin
quasiparticleswhose excitation energyEq(ip) is lowerthan numerical tests using µ = 1.2 fm taken from the Gogny
some cut-off energy Ec in constructing the ground state force[23]. However,this cut-offschemehasandisadvan-
[3]. Namely Eq. (4) is modified to tage that the HFB super matrix in Eq. (38) cannot be
defined,becauseinEq.(44)thepairingHamiltonian(the
NB
Ψ = θ E E(i) b 0 , (40) partproportionaltouiviinsquarebrackets)isdependent
| i c− qp i| i onthecanonical-orbitalindexi. Thismeansthatthereis
iY=1 (cid:16) (cid:17) nocommonpairingHamiltoniantoallthequasiparticles.
where θ is the Heaviside (step) function.
Consequentlyone does nothave well-definedBogoliubov
Inthispaperwedonotdiscusshowoneshouldchange
quasiparticlestatesand cannotutilize them to refine the
the strength as a function of Ec for renormalization. In- HFB solution, which is a ratherserious drawbackas dis-
stead, we regardE as one of the constants to define the
c cussed in the next section. If one can use N in place
w
force. of E , however, one can recover a unique quasiparticle
c
In order to implement the cut-off, one should pay at-
Hamiltonian.
tentiontothefollowingtwopeculiaritiesofthecanonical-
Incidentally, a similar kind of ambiguity arises from
basis HFB method. i) The introduction of E causes a
c the cut-off scheme in terms of the quasiparticle energies
serious inconvenience to this method, i.e., the absence
accordingtoEq.(40). Inthiscase,thecanonicalorbitals
of the quasiparticle Hamiltonian. For this reason, it is
do not have well-defined Hamiltonians. It may not mat-
preferabletousethenumberofcanonicalorbitalsN in-
w ter,however,becauseonedoesnotneedHamiltoniansto
stead of E . ii) For delta-function pairing interactions,
c define canonical orbitals when they can be obtained by
substitutionofE withN leadstoaphysicallymeaning-
c w transforming the quasiparticle states.
lesssituation,inwhichcanonicalorbitalsarecollapsedto
ii)Letusdiscussthephenomenonofthepointcollapse.
aspatialpoint(apointcollapsephenomenon,hereafter).
In quasiparticle (or two-basis) HFB method, specifying
To overcome this difficulty, we modify the pairing inter-
thenumberofquasiparticle(HFsingle-particle)statesto
actionbyaddingarepulsivemomentum-dependentterm.
be consideredhas the same consequencesas imposingan
In the restofthis sectionwe discuss these two problems.
energycut-offbecausethestatestobechosenareanyway
i) Let us explain how the introduction of the cut-off
the eigenstatesof the lowestenergiesofthe quasiparticle
causes the absence of the quasiparticle Hamiltonian. An
(HF) Hamiltonian.
analogywiththe BCSapproximationindicatesanatural
On the other hand, the canonical-basis HFB method
wayto introduceanenergycut-offto the canonical-basis
is very peculiar concerning this point. The canonical or-
method. Namely,thesmoothcut-offmethod[15]modifies
bitalsarenottheeigenstatesofasingleoperatorbuteach
the seniority interaction as of them has its own Hamiltonian as a function of the
i
H
occupation probability given by Eq. (29). becomes h
i
H
vˆpair =−G fi a†ia¯†ı! fj a¯aj, (41) tainvdely˜h.inIttihseotnwlyo ewxhterenmoenesitcuoantsiiodnesrsvi2a=s1maanlldn0u,mrebsepreoc-f
i>0 j>0
X X orbitals that all of them are the lowest-energy approxi-
 
where fi = f(ǫi) is a function of single-particle energy mate eigenstates of h. Otherwise, some of the orbitals
ǫi and takes on 1 well below a chosen cut-off energy maybecometheeigenstatesofh˜ ratherthanh: Weshow
∼
and smoothly becomes zero above it. The smoothness is in Appendix B that, if the occupation probability of an
indispensable to the stability ofthe solution. Inanalogy, orbital becomes sufficiently small, it spontaneously de-
we may modify the pairing density as creases the probability further toward zero (to decrease
the total energy of the nucleus) and change its Hamilto-
ρ˜=4 Nw u v ψ 2 4 Nw f u v ψ 2 (42) nian to h˜. This occurs when the pairing force is a delta
i i i i i i i
| | → | | function (with or without particle-density dependence),
i=1 i=1
X X
7
for which ˜h is a local potential and thus its eigenstates V. GRADIENT METHOD FOR
aredeltafunctionsinthespatialcoordinate. Thismeans CANONICAL-BASIS HFB
that their h is infinitely high. On the other hand, the
h i
rest of the orbitals are left in low-energy subspace of h. Weexplaininthissectionthe detailsofaprocedureto
ThereforetherestrictiononNwresultsinutterlydifferent obtain the canonical-basis solution of the HFB equation
solutions from those obtained by using a cut-off energy directly,notbywayofquasiparticlestates. Whatshould
explicitly. bedoneistominimize E ofEq.(22)underconstraints
tot
ofEqs.(12)and(13). Equivalently,onemay introduce a
Thispeculiarityofthecanonical-basismethodisessen-
Routhian,
tial to the implementation of the method. In the quasi-
particleandthetwo-basismethods,theunphysicaldiver- Nw
genceofthepairingcorrelationstrengthoccuronlyinthe R = E 4ǫ v2
tot− F i
limit of Nw → ∞ (supposing Nw is used instead of Ec). Xi=1
In the canonical-basis method, in addition to this diver- Nw Nw
gence, a different type of unphysical behavior (the point 4 λ ψ ψ δ , (45)
ij i j ij
− {h | i− }
collapse)canalsohappenforN fixedatfinite numbers.
w i=1j=1
XX
One has to take special care of the latter problem if one
tries to implement the canonical-basis method. and minimize it without constraints. ǫF is probably the
mostfamiliarLagrangemultiplier,whosephysicalmean-
Now, if one can suppress the point collapse, one is al- ingistheFermilevel. Inthedefinition(45),N2 Lagrange
w
lowed to use N as the control parameter of the cut-off multipliers λ obeying hermiticity,
w ij
energy. Appendix B shows that this can be achieved,
e.g., by introducing a sufficiently strong repulsive mo- λij =λ∗ji, (46)
mentum dependence to the delta-function pairing inter-
action. With this term applied to finite (small or large) areintroducedinsteadof 21Nw(Nw+1)independentmul-
numbers of canonical orbitals, unphysical situations in- tipliers. The hermiticity ensures the equality between
cluding the point collapse never happen. (If it is applied the number of constraints and the number of indepen-
to infinite number of orbitals, not the point collapse but dent multipliers. Moreover, it makes R real so that two
the divergenceoccurs[22].) Inthis paperweemploythis conditions, δR/δψi =0 and δR/δψi∗ =0, become equiv-
momentum-dependent force. alent and thus one has to consider only one of them.
Note that δ is subtracted from ψ ψ , in contrast to
ij i j
h | i
In addition to the role to enable a cut-off in terms of Ref.[11]. Thissubtractionisnecessarytotreatλ notas
ij
N , equallyimportantis the byproductthatit alsoclar- constantslikeǫ butasfunctionalsofthewavefunctions.
w F
ifies the nature of high-lying canonical orbitals by pro- The stationary condition of R results in two kinds of
viding a kinetic energy term to the pairing Hamiltonian equations. One is ∂R/∂v =0, which concerns the occu-
i
h˜. Further explanations using a numerical example are pation amplitudes v and is fulfilled by [11]
i
given in Sec. VIB.
1 1 ǫ ǫ
v2 = i− F . (47)
Incidentally, Appendix B also gives a discussion on an i 2 ± 2 (ǫ ǫ )2+∆2
alternative way in terms of a pairing-density dependent i− F i
force to enable cut-off in terms of N . Among the double sign,pthe minus sign corresponds to
w
the minimum. The relative sign between v and u =
i i
Before ending this section, we mention the prospect 1 v2 should be plus, i.e., u v 0 (unless ∆ <0).
about the renormalization of the pairing force strength. ±The st−atiionarycondition for a waivie≥function ψ (r)i leads
i
Considering the recent quantitative success of the regu- topthe following equation:
larization method of delta-function forces based on the
Thomas-Fermiapproximation[24,25,26,27],wethinkit 1 δR Nw
is worth trying in the canonical-basis method in future. 4δψ∗ = Hiψi− λijψj
i j=1
X
However, it will require a special treatment. The
Nw Nw δλ
Thomas-Fermi approximation makes the interaction jk
ψ ψ δ =0, (48)
strength a decreasing function of the particle density. − δψ∗ {h j| ki− jk}
j=1k=1 i
One might expect this modification of the interaction XX
would hinder the point collapse. Actually, dependences Owing to the state dependence of , the orthogonal-
i
H
on particle density cannot prevent it because the col- ity condition becomes nontrivial to fulfill. In HF, the
lapsemakesonlythepairingdensitydivergentbutleaves orthogonality is automatically satisfied because ψ are
i
the particle density unchanged (because v w3/2 and eigenstates of the same hermite operator h, which are
ψ (0)2 w−3 asshowninAppendix B). T∝hereforeone orthogonaltooneanother. The orthogonalizationproce-
w
| | ∝
has to keep the momentum-dependent term. This will dureisneededonlybecauseoftheinstabilityfordecaying
lead to a modification of the regularizationprocedure. into Pauli-forbidden configurations. On the other hand,
8
in the canonical-basis HFB method, the orthogonaliza- In order to justify the negligence of the higher-order
tion is indispensable and the explicit functional form of terms in Eq. (55), it must hold[28]
λ is the most important secret of the method.
ij
Reinhard et al. have proposed [11] fevo
f <2 for ∆t = , (56)
evo
T
1 max
λ = ψ ( + ) ψ . (49)
ij j i j i
2h | H H | i where T is the maximum kinetic energy (used as
max
Let us give our reasoning on how the above form can the upper bound of the single-particle energy). For the
be deduced. This consideration plays the crucial role in Cartesian mesh representation with a mesh spacing a,
order to modify the form for the sake of a faster conver-
π 2
gence. The requirement that Eq. (48) must hold at the T =3 B(ρ), (57)
max
a
solution (where ψ ψ =δ ) means,
h i| ji ij (cid:16) (cid:17)
where B(ρ) is given by Eq. (32). By choosing ρ between
Nw
ψ = λ ψ (at the solution). (50) ρ=0(freespace)andρ=0.16fm−3 (nuclearmatterden-
i i ik k
H | i | i sity), one obtains T = 1254 MeV for the SIII force
k=1 max
X with a = 0.8 fm−1. For the calculations shown in this
By taking the overlap with hψj|, one obtains, paper,we haveusedfevo =1.0,i.e., ∆t =5×10−25 sec.
Kinetic energy also arises from the pairing Hamiltonian
λ = ψ ψ (at the solution). (51)
ij h j|Hi| ii h˜. However,an estimation by replacing B(ρ) of Eq. (57)
Eqs.(49) and(51)areequivalentatthe solutionbecause with B˜(ρ˜) of Eq. (34) results in an order of magnitude
λ is defined to be hermite by Eq. (46). However, one smaller value.
ij
shouldemployEq.(49)ratherthanEq.(51)becausethe Whatispresentedabovemaybecalledana¨ıveversion
former,butnotthelatter,ishermiteatpointsotherthan ofthegradientmethod. TheconvergencetotheHFBso-
the solution. lution turns out to be very slow with this na¨ıve method.
OnecanutilizethegradientmethodtoobtaintheHFB Anumericalexamplewillbe giveninSec.VIA. Theori-
solutions in the canonical-basis formalism. Let us de- ginof this slownesscanbe tracedback to the factor u v
i i
scribe this methodbriefly: Smallvariationsinψ andψ∗ in , whichis justa linearcombinationofh andh˜ with
change the Routhian as, i i coeHffiicients v2 and u v . The effects of can be very
i i i Hi
weak for canonical orbitals having small v , which leads
δR δR i
δR= δψ + δψ∗, (52) to the smallness of their changes in a gradient-method
δψ i δψ∗ i
i i step.
Now, steepest-descent paths depend on the choice of
which indicates the steepest descent direction to be,
thevariables. Forexample,Eq.(54)isobtainedwhenone
δψ δR , δψ∗ δR. (53) uses ψi as the variables. If one uses norm-resized wave
i ∝−δψi∗ i ∝−δψi functions χi = ψi/√αi, where αi is a scaling factor, a
gradient-methodstep becomes δχ δR/δχ∗, which is
One takes a small distance movement toward this direc- i ∝− i
equivalent to
tion, i.e., ψ ψ +δψ with
i i i
→
Nw δR
δψ = ∆tα = α ∆t ψ λ ψ .
δψi =−∆t Hiψi− λijψj. (54) i − iδψi∗ − i Hi i− j ij j
 Xj=1   X (58)
By repeatingthe evaluationof δψi andthe movement, Thus the path is changed. By choosing αi to cancel
{ }
one eventually reaches the minimum of R. Incidentally, the small factors in , one can accelerate the otherwise
i
H
in the braces on the right-hand side, the term including slow convergence,which we call the acceleratedgradient
δλ /δψ∗ in the second member of Eqs. (48) is dropped method.
jk i
because we know empirically that the orthogonality is In this paper, we take
stable without this term.
In Eq. (54), ∆t is a parameter introduced to control 1 facc
α min , , (59)
thesizeofamovement. Oneshouldregardthatthepref- i ≃ v2 u v
(cid:18) i i i(cid:19)
actors on the right-hand side of Eq. (39) are included in
this parameter. One may call ∆t the imaginary-time- wherethefactorfacc maybechosenempiricallytomaxi-
step size because the first order approximationin ∆t of mizethespeedofconvergence. Wechoosefacc =5. This
theimaginary-timeevolutionmethod[29]becomesessen- factor seems to take care of the fact that h˜ is smaller
tially the same procedure as the gradient method : than h by an order of magnitude. With this choice for
the scaling factor, α h for deeply bound orbitals
δψi = e−∆tHiψi ψi and α 5h˜ for hiigHhi-l∼ying orbitals. Thus, all the or-
− i i
= ∆t ψ + ∆t2 . (55) bitals eHvol∼ve at roughly the same pace.
i i
− H O
(cid:0) (cid:1)
9
Incidentally,byusingtheapproximateequalitysymbol renewedset of canonicalorbitals as the eigenvectorsand
in Eq. (59), we mean the utilization of some common v2 as the eigenvalues. If the initial set of canonical or-
i
recipes for numerical calculations like taking a moving bitalsaretakenfromtheexactsolution,therenewedand
average and restricting on the maximum values. the initial sets are identical. However, if the initial set
It should be noticed here that one should also mod- is taken before the state converges to the solution, the
ify the form of multipliers λ when one introduces the renewedsetcorrespondstoabettersolutionthantheini-
ij
acceleration factors α = 1. In principle one may reach tialset. Itisbecausetheaboveproceduregivestheexact
i
6
thesolutionbyusingEq.(49)becausetransformationsof variationalminimum[21,30]inthesubspacespannedby
variables do not change the location of the minimum of the initial canonical orbitals. The gradient-step method
R. Nevertheless,oneshouldmodifyEq.(49)forapracti- takes care of both the variation inside this subspace and
calreason. Incalculatingthegradientvectorwhoseexact the optimization of the subspace itself. The diagonal-
formisgivenbythesecondmemberofEqs.(48),thelast ization performs only the former part but perfectly in
termtakesmuchmorecomputingtimethanthefirsttwo a single step. This diagonalization method is not indis-
terms due to δλ /δψ∗. One can drop the last term if pensable but useful to obtain the solutions. It makes
jk i
the orthogonality relation (12) is fulfilled along the path the convergence more robust and somewhat quicker if
of the evolution. Let’s suppose that the relation is sat- it is performed after every 100 gradient-method steps.
∼
isfied before a gradient-methodstep is takenand require Its effect seems to saturate with this interval since using
thatitisconservedwithoutthedroppedtermtothefirst shorterperiodsdoesnotleadtonoticeableimprovements.
order in ∆t after the step, i.e., A numerical example is given in Sec. VIA.
Incidentally, the idea of this diagonalization method
ψ ψ =δ , ψ′ ψ′ =δ + (∆t)2 , (60) originatesinthetwo-basisformalismofHFB[5],inwhich
h i| ji ij h i| ji ij O
the quasiparticleHamiltonianis diagonalizedin the sub-
(cid:16) (cid:17)
withψ′ =ψ +δψ andδψ givenbyEq.(58). Substitut- space spanned by low-lying HF orbitals (eigenstates of
i i i i
ing ψ and ψ′ in Eq. (60) and requiring the hermiticity h), not the canonical orbitals as in this paper.
i i
(46) result in Therelationbetweenthese quasiparticlesandthe true
quasiparticles defined in the full space is discussed in
1 Sec. VID.
λ = ψ (α +α ) ψ . (61)
ij j i i j j i
α +α h | H H | i
i j
This new form of λ , as well as the original form of
ij VI. RESULTS OF NUMERICAL
Eq. (49), fulfills the requirement that it should agree
CALCULATIONS
with the expression (51) at the solution. The original
form differs from the new form, however, before reach-
In this section, we show the results of numerical cal-
ing the solution if α = α . Consequently, only the new
i j
6 culations using a newly developed computer program
formconservestheorthogonalityduringthecourseofthe
[8, 12, 13, 14] of the canonical-basis HFB method in a
evolution.
three-dimensional Cartesian mesh representation. The
We have indeed suffered from large errors of orthog-
parameters of the model are as follows: For the mean-
onality by using the original form. On the other hand,
field interaction, Skyrme SIII force is employed but its
by using the new form, we have observed that the er-
spin-orbit term and the Coulomb force are not included.
ror does not grow but decreases without performing ex-
For the pairing interaction, ρ = 0.32 fm−3 and ρ˜ =
plicit orthogonalizations. The reason for this stability c c ∞
areused. The value ofv is changeddepending onk (=
will probably be found in the second order terms in ∆t p c
1.7 – 5 fm−1) and N (= 14 – 210) to make the average
neglected in Eq. (60). w
pairinggaptobe2–2.5MeV.Unlessthevaluesarespec-
There is an auxiliary method to speedup the conver-
ified,k =2fm−1,v = 1050MeVfm3,andN =21are
gence. It is to diagonalize the HFB super matrix given c p − w
used. With N =21, the number of orbitals amounts to
by Eq. (38) in the subspace spanned by N canonical w
w
three times as large as the number of the nucleons,since
orbitals. Then one obtains N quasiparticle states with
w eachorbitalcanbe occupied by four nucleons due to the
positiveenergy,whosetwo-componentwavefunctionsde-
spin-isospin degeneracy.
fined in Eq. (5) are expressed as,
The calculated state is the axially symmetric prolate
Nw Nw solutionof 2184Si14. The single-particlewave functions are
φ (r)= ψ (r)U , ϕ (r)= ψ (r)V , (62) expressedonaCartesianmeshhaving(in the defaultset
i j ji i j ji
up) 29 29 29 mesh points with mesh spacing a=0.8
j=1 j=1
X X × ×
fm. As the boundary condition, we assume a cubic box
where (U , ,U ,V , ,V )T is the ith normal- whoseinfinitely highwallsarelocatedatthe0thandthe
1i ··· Nwi 1i ··· Nwi
ized eigenvector of Eq. (38). The HFB ground state is 30th mesh points. The center of mass of the nucleus is
constructed as the vacuum of these quasiparticles as in placed at the center of the box. The 17-point formu-
Eq. (4). By diagonalizing the one-body density matrix lae [31] are employed to calculate the first and second
of this vacuum in this subspace, V∗VT, one obtains a derivatives. In applying these formulae, wave functions
10
are assumed to be antisymmetrically extended beyond
the walls. Spatial integrals are done using the mid-point
(=trapezoidal) rule.
A. Convergence to HFB solutions
InFig.2,examplesoftheconvergencecurvesareshown
for severalversionsofthe gradientmethod. The calcula-
tions are done using the default parameters. The initial
wave functions of the canonical orbitals are taken from
theeigenstatesofthestandardharmonicoscillatorpoten-
tialwithaquadrupledeformationofβ =0.3andγ =15◦.
The D symmetryis conservedtoa verygoodaccuracy
2h
during the course of the gradient-method evolutions.
For the size of an imaginary time step, we use f =
evo
1.0, defined by Eq. (56), for all the curves. This value
of f roughly optimizes the convergence speed for the
evo
fastestmethod. Forslowermethods,f maybeslightly
evo
largerbutsuchasmallchangedoesnotaffectourdiscus-
sions below.
In connection with the step size, we also adopt a
recipe from the HF+BCS method of Ref. [15], accord-
ing to which the changes of h and ˜h after each gradient-
method step should be further damped compared with
the changes of ψ (by a factor of 0.4 in our case). It is a
i
precautionagainstoccasionalinstabilities like oscillating
behaviors. This additional damping does not affect our
conclusion on the comparison of the methods, either.
The top portion of Fig. 2 shows the total energy E
tot
of Eq. (22). The middle portion shows the maximum
over i = 1, ,N of the error of the second equality of
w
···
Eqs. (48) (neglecting the contribution from the error of
orthogonality),
FIG. 2: Convergence to HF and HFB solutions for 28Si. See
Nw text for explanations.
∆ǫHFB = ψ λ ψ , (63)
i Hi i− ij j
(cid:12) j=1 (cid:12)
(cid:12) X (cid:12)
(cid:12) (cid:12)
which is an indicator of(cid:12) the accuracy of(cid:12)HFB solutions. for other quantities like the axial asymmetry γ and the
AsforHF,thisquantityisreducedtothemaximumover r.m.s.radius. Thereisaseriousproblem,however,ofthe
i=1, ,1A of obvious slowness of the convergence.
··· 4 Dot curves show the convergence to an HF solution
∆ǫHF = hψ ǫ ψ = ψ h2 ψ ψ hψ 2, (64) withpairinginteractionsturnedoff. (With“w/D”inthe
i | i− i i| h i| | ii−h i| | ii legend in the top portion of the figure we mean that di-
which is just the energypwidth of ψ for h. The bot- agonalization of h is done in the subspace of occupied
i
| i
tomportionshowsthe massquadrupoledeformationpa- 1A = 7 HF orbitals to renew the orbitals after every
4
rameter β defined as the general definition [21] times 100 gradient-method steps. Without this diagonaliza-
45/16π 0.95 disregarding nucleon form factors. HF tion,ψi arenotHForbitalsanddonotsatisfy∆ǫHi F =0.
≃
resultforβ,convergingto0.44,isomitted. Theabscissae This diagonalization also makes the convergence some-
p
arecommontoalltheportionsanddesignatethenumber what quicker.) One can see that the HF energy reaches
of gradient-method steps. a plateau by three orders of magnitude faster than the
Thin solid lines in all the three portions are the re- na¨ıveHFBmethod. Inthemiddle portion,whileHFcan
sult of the na¨ıve (i.e., with α = 1) gradient method. achieve precision of 0.1 keV with only 1400 steps, HFB
i
These curvesdemonstratethat onecanindeedobtainan requires 55000 steps for 10 keV, 3 105 steps for 1 keV,
×
HFB solution with the canonical-basis HFB method on and 1.5 106 steps for 0.1 keV precisions.
×
a mesh, because after a million steps, E and β appear A method to improve the convergence speed is the di-
tot
to reach a plateau and ∆ǫHFB becomes as small as 0.1 agonalization of the HFB super matrix in the subspace
i
keV. We have also obtained similar convergence curves spanned by N canonical orbitals. Dot-dash curves are
w