MICROSCOPIC CALCULATION OF COLLECTIVE HAMILTONIAN FOR COMPLEX
NUCLEI
By
Liyuan Jia

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Physics
2012

ABSTRACT
MICROSCOPIC CALCULATION OF COLLECTIVE HAMILTONIAN FOR
COMPLEX NUCLEI
By
Liyuan Jia
A formalism is proposed for microscopic calculation of the collective nuclear Hamiltonian
from the underlying nucleonic Hamiltonian, based on the generalized density matrix. It goes
beyond the well-known random phase approximation, which calculates only the harmonic
terms in the collective Hamiltonian. Many nuclei, including those at the γ-unstable and
rotational limits, have a vanishing or negative harmonic potential; there the anharmonic
terms calculated by the generalized density matrix method are indispensable restoring the
stability of the system. The proposed formalism is tested in several models by comparing
with the exact results. Realistic applications on tin isotopes are discussed.

ACKNOWLEDGMENT

I gratefully express many thanks to my advisor Vladimir Zelevinsky. He helped me in
his best to grow as a physicist, starting from an engineering undergraduate major. His
understanding for the student, especially during the difficult times, makes my PhD years an
enjoyable journey.
I am grateful to people at the cyclotron and the physics department. It is a very good
working environment here. I received numerous help from professors, staff members, and
fellow students here. Without them, it would be impossible to finish this work.
Support from friends is always important. It is fortunate for me to meet some new friends
here at MSU; I believe our friendship will continue.
At last, thanks to my parents for their support, as always.

iii

TABLE OF CONTENTS

List of Tables . . . . . . . . . . . . . . .
List of Figures
1 Introduction

. . . . . . . . . . . . . . . . . . .

vi

. . . . . . . . . . . . . . . . . .

vii

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

. . . . . . . . . . . . . .

2 Generalized Density Matrix Method
. . . . . . . .
2.1 Development . . . . . . . . . . . . . . . . . . . . . .
2.2 Hartree-Fock Equation . . . . . . . . . . . . . . . . .
2.3 Random-Phase Approximation . . . . . . . . . . . . .
2.4 Higher Orders . . . . . . . . . . . . . . . . . . . . .
2.5 Self-consistent Hamiltonian Condition . . . . . . . . .
2.6 Equivalent Representations of Collective Hamiltonian
2.7 Transition Rates . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.

4
4
8
9
10
13
14
17

3 Validity of the Formalism
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1 Lipkin Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Models with Factorizable Force . . . . . . . . . . . . . . . . . . . . . . . . .

19
19
26

4 Formalism for Realistic Nuclei
. . . . . . . . . . . .
4.1 Formalism . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Equations of Motion for Quadrupole Phonon . . . . .
4.2.1 Hartree-Fock-Bogoliubov Equation . . . . . .
4.2.2 Quasi-Particle Random Phase Approximation
4.2.3 Cubic Order . . . . . . . . . . . . . . . . . . .
4.2.4 Quartic Order . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

38
39
42
44
45
45
46

5 Quadrupole-plus-Pairing Model . . . . . . . . . . . .
5.1 The GDM Method . . . . . . . . . . . . . . . . . . .
5.1.1 BCS Theory . . . . . . . . . . . . . . . . . . .
5.1.2 Quasi-Particle Random Phase Approximation
5.1.3 Cubic Anharmonicity . . . . . . . . . . . . . .
5.1.4 Quartic Anharmonicity . . . . . . . . . . . . .

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

49
51
51
53
54
55

iv

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

5.2
5.3

Comparison with Exact Results . . . . . . . . . . . . . . . . . . . . . . . . .
Tin isotopes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6 Conclusion

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

A Values of g L and γ L

65

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

72

l,l

B Conventions

55
60

l,l

C Systems Near Critical Point . .
C.1 Harmonic Order . . . . . . . . .
C.2 Cubic Anharmonicity . . . . . .
C.3 Quartic Anharmonicity . . . . .
D Mode Coupling
E Odd-Mass Nuclei

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

74
75
75
77

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

78

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

80

Bibliography . . . . . . . . . . . . . . .

v

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

. . . . . . . . . . . . . . . . . . .

83

LIST OF TABLES

Table 5.1

Results of the quadrupole plus pairing model for different values of
the pairing strength G. All quantities are in units of MeV. ∆ is the
solution of Eq. (5.12). The chemical potential µ is the solution of Eq.
(5.16). κc is the critical κ such that ω 2 in Eq. (5.18) becomes zero.
Λ(30) is given by Eq. (5.20). Λ(40) is given by Eq. (5.22) setting
Λ(12) = 0. “GDM E2+ ” is the excitation energy of the first 2+ state
(22)

from diagonalizing Eq. (4.18) for ω 2 = Λ(12) = Λ(04) = ΛL = 0.
“NuShellX E2+ ” is the exact excitation energy of the first 2+ state
from diagonalizing Eq. (5.5), in which G and κ are given by G and
κc in the table. Similarly, “NuShellX E4+ ” is the exact excitation
energy of the first 4+ state. . . . . . . . . . . . . . . . . . . . . . .
Table 5.2

57

Neutron single-particle energies for 100 Sn and 132 Sn in units of MeV.

60

vi

LIST OF FIGURES

Figure 3.1

The excitation energies En − E0 of the first five excited states as a
function of κ in the model (3.1) with Ω = 30. The red dashed-dotted
line is the RPA calculation (diagonalizing the harmonic Hamiltonian
ω 2 α2 + 1 π 2 ). The blue dashed line is obtained by diagonalizing
2
2
ω 2 α2 + 1 π 2 + Λ(40) α4 in the infinite phonon space {|0 ≤ n < +∞ },
2
2
4!
where ω 2 and Λ(40) are given by Eq. (3.10). The black solid line shows

the exact results from diagonalizing the original fermionic Hamiltonian (3.1) directly. For interpretation of the references to color in this
and all other figures, the reader is referred to the electronic version
of this dissertation. . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.2

23

The errors of the first four excitation energies E −Eexact as a function
of κ in a set of calculations. Panels (a), (b), (c), and (d) are for the
first, second, third, and fourth excitation energy, respectively. Four
lines on each panel are obtained by diagonalizing different collective
Hamiltonians in different phonon spaces. 1 (green dashed line): H =
ω 2 α2 + 1 π 2 + Λ(40) α4 in {|0 ≤ n < +∞ }. 2 (red dotted line):
2
2
4!
ω 2 α2 + 1 π 2 + Λ(40) α4 in {|0 ≤ n ≤ Ω }. 3 (blue dashedH = 2
2
4!
ω 2 α2 + 1 π 2 + Λ(40) α4 in {|0 ≤ n ≤ Ω }, but with
dotted line): H = 2
2
4!
ω 2 = 1 − κ2 Ω(Ω + 2) replacing that in Eq. (3.10). 4 (black solid line):
2

(04)

(40)

H = ω α2 + 1 π 2 + Λ 4! α4 + Λ 4! π 4 in {|0 ≤ n ≤ Ω }, also with
2
2
2 = 1−κ2 Ω(Ω+2). 1 and 2 closely overlap and are indistinguishable
ω
in the figure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3.3

24

Single-particle level scheme. ˜ is the time-reversal level of 1. Each
1
level has a quantum number m. . . . . . . . . . . . . . . . . . . . .

26

vii

Figure 3.4

Figure 3.5

Figure 3.6

Figure 3.7

Figure 5.1

Figure 5.2

Excitation energies En − E0 and transition matrix elements m|Q|n
in model 1 as a function of κ. The black solid lines are exact results of
shell-model diagonalization. The black dashed line is the beginning
of the single-particle continuum. The red dashed-dotted lines are
the RPA results. The blue symbols are the GDM results. The stars
are energies; the squares, circles, up-triangles and down-triangles are
matrix elements of Q between different states. m|Q|n that are not
shown vanish in both the shell model and the GDM calculations. . .

29

Excitation energies En − E0 and transition matrix elements 0|Q|1
in model 2 as a function of κ. The black solid lines are exact results of
shell-model diagonalization. The black dashed line is the beginning
of the single-particle continuum. The red dashed-dotted lines are the
RPA results. The blue symbols (stars, up-triangles and circles) are
results of three different sets of GDM calculations, as labeled in the
legend. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

Excitation energies En − E0 and transition matrix elements m|Q|n
in model 3 as a function of κ. The black solid lines are exact results of
shell-model diagonalization. The black dashed line is the beginning
of the single-particle continuum. The red dashed-dotted lines are the
RPA results. The blue symbols are the GDM results. The stars are
energies; the squares, circles and up-triangles are matrix elements of
Q between different states. . . . . . . . . . . . . . . . . . . . . . . .

34

Excitation energies En − E0 and transition matrix elements m|Q|n
in model 4 as a function of κ. The black solid lines are exact results of
shell-model diagonalization. The black dashed line is the beginning
of the single-particle continuum. The red dashed-dotted lines are the
RPA results. The blue symbols are the GDM results. The stars are
energies; the squares, circles and up-triangles are matrix elements of
Q between different states. . . . . . . . . . . . . . . . . . . . . . . .

35

Excitation energies (from Table 5.1) in the quadrupole plus pairing
model as a function of the pairing strength G, at the critical point
ω = 0 (κ = κc ). The black squares and red circles show the exact
excitation energy of the first 2+ and 4+ state, respectively, “NuShellX
E2+ ” and “NuShellX E4+ ” in Table 5.1. The blue triangles give
“GDM E2+ ” from Table 5.1. . . . . . . . . . . . . . . . . . . . . .

56

Even-odd mass staggering in the tin isotopes. The black squares show
the experimental data MA −(MA−1 +MA+1 )/2. The red circles show
the BCS pairing gap ∆ as the solution of Eq. (5.12). . . . . . . . .

61

viii

Figure 5.3

Excitation energy E(2+ ) and transition rates B(E2; 0+ → 2+ ) for the
first excited state 2+ in the tin isotopes. The black squares show the
experimental data taken from Ref. [26]. The blue triangles are the
GDM results. The red circles show the QRPA results. On the upper
panel, negative values for the QRPA E(2+ ) actually mean imaginary
frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

63

Chapter 1
Introduction
The nuclear many-body problem aims at describing a nucleus in terms of nucleonic degrees of
freedom. It is very difficult in at least two aspects. First, the interaction between nucleons, as
the residual force among their quark components, is still not firmly established. Second, exact
solution is impossible except for very light nuclei, as is common in many-body physics. The
non-perturbative nature of nucleon force makes the nuclear problem particularly difficult.
Until now the most effective treatment starts from a self-consistent mean field, generated
by the participating nucleons themselves. In the first step, nucleons move independently in
the mean field obeying the Fermi statistics. In the next step, the nucleons interact with
each other through the residual force that is not included into the mean field. The standard
treatment is given by the nuclear shell model (configuration interaction). There the residual
interaction is diagonalized in a huge many-body space that consists of Slater determinants
built of single-particle levels. However, the dimension of the space makes it impractical in
cases of many active nucleons in lots of relevant single-particle levels, which is typical for
medium and heavy nuclei.

1

The collective motion is an important type of dynamics in nuclei. Intuitively, it can be
thought of as the motion of the nucleus as a whole, such as vibration and rotation, where
many nucleons move coherently (in phase). Typical features of the collective states are their
low excitation energies and large transition rates. For these collective states, another type
of description using an effective bosonic Hamiltonian is successful, including the well-known
geometric Bohr Hamiltonian [1, 2] and interacting boson model [3]. Broad sets of nuclear
data are described with a few parameters that change smoothly across the nuclear chart. This
shows that, out of the huge Slater-determinant space, there exist a few collective degrees of
freedom, which are usually sufficient for describing the low-lying collective states. To put the
phenomenological theory on solid grounds, serious efforts were devoted to calculating these
parameters microscopically from the underlying nucleonic Hamiltonian. But the complete
theory is still missing after several decades.

The purpose of this thesis is the microscopic calculation of the collective/bosonic Hamiltonian by the generalized density matrix (GDM) method. This method was proposed long
ago [4, 5, 6] and was applied to nuclear rotation [5, 7, 8] and large-amplitude collective motion [9, 10, 11]. Here the construction is considered in the most general way. As a result,
the whole nuclear chart, from vibrational, through γ-unstable, to rotational regions, can be
described on the same footing. From another perspective, the approach greatly generalizes
the well-known random phase approximation (RPA). The lowest orders of the GDM method
give naturally Hartree-Fock (HF) equations and RPA; higher orders fix the anharmonic terms
in the collective Hamiltonian. Only near the closed shells, where the RPA frequency ω 2 is
large, the harmonic potential term dominates and the RPA is a good approximation. Going
away from closed shells, ω 2 becomes small and finally negative. There the anharmonic terms

2

are indispensable restoring the stability of the system.
The thesis is organized as following. In chapter 2 we present the GDM formalism in the
simplest scenario, without complications owing to angular-momentum vector coupling. The
validity of the GDM method is then examined in chapter 3 by comparing with the exact
shell-model results in several models. Chapter 4 is devoted to applications for a realistic
nuclear system with rotational symmetry and pairing. Finally, we consider in chapter 5 a
specific Hamiltonian – the quadrupole plus pairing Hamiltonian, with realistic calculations
for tin isotopes. Some details of the calculations are given in Appendices. The results are
partly published [12, 13, 14].

3

Chapter 2
Generalized Density Matrix Method
In this chapter we review the essence of the GDM method in simple systems without complications owing to rotational symmetry and pairing correlations. A single collective mode
is assumed; the case of multiple modes is discussed briefly in Appendix D.

2.1

Development

The starting point is the microscopic fermionic Hamiltonian
1
†
12 a1 a2 +

H=
12

4

Here the operators of single-particle energy

† †

V1234 a1 a2 a3 a4 .

(2.1)

1234

and residual interaction V are expressed using

an arbitrary basis. We find it convenient for H, 12 and V1234 to be dimensionless; in other
words H is measured in some unit of energy. We have assumed in Eq. (2.1) a two-body force,
inclusion of three-body forces is discussed in Appendix A of Ref. [12]. In accordance with
the discussion in Introduction, we consider a typical case when the system described by the

4

Hamiltonian (2.1) has a band of collective states {|Ci }, characterized by low energies and
large transition rates. We assume that there exists a reference state |Φ , a collective mode
√
operator A† = (uα + ivπ)/ 2 (uv = −1, α and π are collective coordinate and momentum),
such that approximately

[α, π] = i,
|Ci = [ ci0 + ci1 A† + ci2 (A† )2 + . . . ] |Φ ,
ω 2 2 1 2 Λ(30) 3 Λ(12)
α + π +
α +
{α, π 2 }
2
2
3!
4
Λ(40) 4 Λ(22)
Λ(04) 4
+
α +
{α2 , π 2 } +
π + . . . |Cj .
4!
8
4!

(2.2)
(2.3)

Ci | H |Cj = Ci | E0 +

(2.4)

Eq. (2.2) says that A† is effectively a boson operator. Eq. (2.3) says that the collective
band {|Ci } can be built by repeated action of A† on the reference state |Φ . Later |Φ
will be identified as the HF ground state. Eq. (2.4) says that within the band, the effect
of the fermionic Hamiltonian (2.1) can be approximated by an expansion over the bosonic
operators, where all time-even terms are kept (α is time-even, π is time-odd). The curly
brackets represent anti-commutator, {A, B} ≡ AB + BA.

The generalized density matrix operator is defined as

†

R12 ≡ a2 a1 .

(2.5)

Similarly to Eq. (2.4), we assume that within the band the effect of R can be approximated

5

by a boson expansion

(10)

(20) α

(01)

Ci |R12 |Cj = Ci | ρ12 + r12 α + r12 π + r12
(30)

+r12

α3
3!

(03)

+ r12

π3
3!

(21)

+ r12

{α2 , π}
4

2

2

(02) π

+ r12

2

(11) {α, π}

+ r12

2
2
2}
(12) {α, π
+ r12
+ . . . |Cj .
4

(2.6)

The expansion coefficients r(mn) are called generalized density matrices. They are matrices
in the single-particle space and do not act in the band of collective states. Later the constant
term ρ in Eq. (2.6) will be identified as the usual density matrix on the HF ground state.
Terms with operators α and π generate interaction within the band.

The exact equation of motion for the density matrix operator in the full many-body
Hilbert space is

[R12 , H] = [ , R]12 −

where [ , R]12 =

1
2

† †

V5432 a5 a4 a3 a1 +
345

1
2

† †

V1345 a2 a3 a4 a5 ,

(2.7)

345

3 ( 13 R32 − R13 32 ). Now we map Eq. (2.7) onto the collective subspace,

namely, we take matrix elements of Eq. (2.7) between two collective states on both sides.
Thus R12 and H in Eq. (2.7) can be substituted by their boson expansions (2.6) and (2.4).
For the two-body density matrix operator, we assume

† †

†

†

†

†

Ci | a4 a3 a2 a1 |Cj = Ci | a4 a1 · a3 a2 − a4 a2 · a3 a1 |Cj ,

(2.8)

that is, it factorizes into antisymmetrized products of one-body density matrix operators
within the band. The constant terms (without operators α and π) in Eq. (2.8), ρ1234 =
ρ14 ρ23 − ρ24 ρ13 , is routinely used to derive the time-dependent Hartree-Fock equation. The

6

validity of Eq. (2.8) was studied in detail in Ref. [12]. Substituting Eq. (2.8) into Eq. (2.7)
we get the equation of motion in the collective band

[R12 , H]

[ + W {R}, R]12 ,

(2.9)

where “ ” means projecting onto the collective subspace (taking matrix elements Ci |...|Cj ).
The generalized self-consistent field is defined as

W {R}12 ≡

V1432 R34 .

(2.10)

34

On the left-hand side of Eq. (2.9), the intermediate states (between R12 and H) are restricted
to those of the collective subspace, assuming large transition amplitudes.
For conciseness we rewrite Eq. (2.4) in compact form,
m+2l≥2

H

Λ(m,2l)

m≥0,l≥0

1 {αm , π 2l }
,
2 m!(2l)!

(2.11)

where Λ(20) = ω 2 , Λ(02) = 1. Similarly Eq. (2.6) is rewritten as
m n
(mn) 1 {α , π }

R12

r12
m≥0,n≥0

2

m! n!

.

(2.12)

The substitution of Eq. (2.12) into Eq. (2.10) gives the mapping of the latter,

w(mn)

W {R}
mn

where w(mn) ≡ W {r(mn) }.

7

1 {αm , π n }
,
2 m! n!

(2.13)

Substituting Eqs. (2.11) and (2.12) into Eq. (2.9), on the left-hand side we calculate
commutators such as [α2 , π] = 2iα. Then, comparing coefficients with the same operator
structure, we come to the final set of GDM equations with different r ≥ 0, s ≥ 0,
p+2l≥2
0≤p(≤r+1),0≤2l(≤s+1)

=

2.2

1
[ , r(rs) ] +
r! s!

2l(r + 1 − p) − (s + 1 − 2l)p
· i Λ(p,2l) r(r+1−p,s+1−2l)
(r + 1 − p)! (s + 1 − 2l)! p! (2l)!

0≤p(≤r),0≤q(≤s)

1
[w(r−p,s−q) , r(pq) ].
(r − p)! (s − q)! p! q!

(2.14)

Hartree-Fock Equation

Equation (2.14) with (rs) = (00) gives

0 = [ + W {ρ}, ρ].

(2.15)

Thus f ≡ + W {ρ} and ρ can be diagonalized simultaneously in some single-particle basis:

f12 = δ12 e1 ,

ρ12 = δ12 n1 ,

(2.16)

providing mean-field single-particle energies e1 and occupation numbers n1 . In what follows,
we will always use this single-particle basis. If we restrict the reference state |Φ to be a
Slater determinant, then the occupation numbers n1 can be only 0 or 1; in this case Eq.
(2.15) is the usual HF equation, |Φ is the HF ground state. More general choices, such as
the thermal ensemble, are also possible. For future convenience we define

e12 ≡ e1 − e2 ,

n12 ≡ n1 − n2 .

8

(2.17)

We assume that degenerate single-particle levels have the same occupancies,

e1 = e2 ⇒ n1 = n2 ,

(2.18)

but the reverse is not necessarily true.

2.3

Random-Phase Approximation

Equation (2.14) with (rs) = (10), (01) gives

ir(10) = [f, r(01) ] + [w(01) , ρ],

(2.19)

−iω 2 r(01) = [f, r(10) ] + [w(10) , ρ],

(2.20)

where w(10) = W {r(10) }, and w(01) = W {r(01) } are the corresponding components of the
self-consistent field. This is the set of RPA equations. The formal solution is

(10)

n12
(01)
(10)
[−iω 2 w12 + e12 w12 ],
2 − ω2
(e12 )
n12
(01)
(10)
(01)
r12 =
[i w12 + e12 w12 ].
2 − ω2
(e12 )

r12 =

(2.21)
(2.22)

Note that r(10) and r(01) have only n1 = n2 matrix elements. From Eqs. (2.10), (2.21) and
(2.22) we obtain a linear homogenous set of equations for w(10) and w(01) :

(10)

w34 =

n12
(10)
(01)
[−iω 2 w12 + e12 w12 ],
2 − ω2
(e12 )
n12
(01)
(10)
V3214
[iw12 + e12 w12 ].
2 − ω2
(e12 )

V3214

12
(01)
w34 =
12

9

(2.23)
(2.24)

A non-zero solution requires a zero determinant, from which we solve for the harmonic
frequency ω 2 . The lowest solution should be taken corresponding to the collective motion.

2.4

Higher Orders

For K = r + s ≥ 1 in Eq. (2.14), we solve a linear set of coupled equations for r(rs) |r+s=K
in terms of lower-order quantities. The formal solution can be written as
p+2l≥3

r(mn) |m+n=K = −
r+s=K 0≤p(≤r+1),0≤2l(≤s+1)

2l(r + 1 − p) − (s + 1 − 2l)p
(r + 1)(s + 1)
(mn)

p

2l
· Cr+1 Cs+1 i Λ(p,2l) · η(rs) : r(r+1−p,s+1−2l)
p+q≤r+s−1

p q

+
r+s=K 0≤p(≤r),0≤q(≤s)

(mn)

Cr Cs · η(rs) : [w(r−p,s−q) , r(pq) ],

(2.25)

(mn)

q

where Cp = p!/[q!(p − q)!], and we have introduced the “weight” matrix η(rs) so that
(mn)

(mn)

(η(rs) : r)12 = η(rs)12 r12 . The matrix η is given by
(mn)

(rs)

−1
η(rs) = (DK )(mn) ,

(2.26)

where DK is a tridiagonal matrix of dimension K + 1, with non-vanishing matrix elements
(mn)

(DK )(mn) = −e,
(m+1,n−1)

= i · n,

(m−1,n+1)

= −i · ω 2 m.

(DK )(mn)
(DK )(mn)

10

(2.27)

We give as an example the D and η matrices in the lowest two orders. For K = 1,




 (mn) \ (rs) (01) (10) 


(mn)

2 ,
D(rs) = 
(01)
−e −iω 




(10)
i
−e

(2.28)

from which we solve by Eq. (2.26)




 (mn) \ (rs) (01) (10) 


1
(mn)


η(rs) = 2
·
(01)
−e −i  .
2 


e −ω


2
(10)
iω
−e

(2.29)

Hence the RPA solutions (2.21) and (2.22) are rederived substituting the above expression
into Eq. (2.25). For K = 2,




(20) 
 (mn) \ (rs) (02) (11)






(02)
−e −iω 2
0
(mn)


D(rs) = 
,

2 

2i
−e −2iω 
(11)




(20)
0
i
−e

(2.30)

from which we solve by Eq. (2.26)



(02)
(11)
(20)
 (mn) \ (rs)



(02)
e2 − 2ω 2
2ie
−2
1
(mn)

η(rs) =
·
−e(e2 − 4ω 2 ) 

(11)
−iω 2 e
e2
ie


−2ω 4
−2iω 2 e e2 − 2ω 2
(20)

11






.





(2.31)

(20)

From Eq. (2.31) we read, for example, η(02)12 = {−2ω 4 }/{−e12 [(e12 )2 − 4ω 2 ]}, where 1
and 2 are single-particle indices. All η matrices with an even K have the factor 1/e [see for
(mn)

example Eq. (2.31)], thus the e1 = e2 matrix elements of r12

cannot be directly calculated

from the solutions (2.25). However, if we set n1 = n2 in Eq. (2.25) and simplify, the 1/e12
(mn)

divergence is canceled; the resulting expression is used to calculate r12
(mn)

matrix elements of r12

. The e1 = e2

are expressed in terms of lower-order quantities.
(mn)

In each order K, substituting r(mn) from Eq. (2.25) into w12

=

(mn)
results
34 V1432 r34

in a linear set of coupled equations for the latter, from which the e1 = e2 matrix elements of
(mn)

w12

are solved in terms of lower-order quantities, which in turn gives the e1 = e2 matrix
(mn)

elements of r12

by Eq. (2.25). However, if K = 2L + 1 is odd, the determinant for w(mn)

is zero. This can be proved as following. Summing Eq. (2.14) with proper weights we get

i x = [W {y}, ρ] + [f, y] + . . . ,
(2.32)
−i ω 2 y = [W {x}, ρ] + [f, x] + . . . ,
where

x=
0≤t≤L

νt
r(2t+1,2L−2t) ,
(2t + 1)! (2L − 2t)!

y=
0≤t≤L

µt
r(2t,2L+1−2t) ,
(2t)! (2L + 1 − 2t)!

in which µt and νt are solved from (0 ≤ t ≤ L, µL+1 = ν−1 = 0)
2L − 2t 2
1
ω µt+1 + µt =
νt ,
2t + 1
2t + 1
1
2t
νt−1 =
ω 2 µt .
ω 2 νt −
2L − 2t + 1
2L + 1 − 2t
−

12

The “. . .” in Eq. (2.32) are lower-order quantities. It is seen that the explicitly shown parts of
Eq. (2.32) have the same structure as the RPA equations (2.19) and (2.20). This finishes the
proof. The zero determinant means that the set of equations for w(mn) is linearly dependent,
and there is a constraint in each order of odd K, entering as a solvability condition. These
constraints are the main results of the GDM formalism, from which the parameters Λ(pq) of
the collective Hamiltonian are calculated. Then w(mn) is solved from this zero-determinant
set, with only a factor undetermined. This factor is fixed by the constraint indicated in the
next section below Eq. (2.33).

2.5

Self-consistent Hamiltonian Condition

If the GDM formalism is self-consistent, the substitution of the solutions (2.25) into Eq.
(2.1) should reproduce the assumed Hamiltonian (2.11),
1
Λ(mn) = Tr [ r(mn) ] +
2

p

q

Cm Cn Tr[r(pq) w(m−p,n−q) ].

(2.33)

0≤p(≤m),0≤q(≤n)

In an order of odd K = m + n, all the parameters Λ(mn) are checked correctly. In an order
of even K = m + n, all but one Λ(mn) are checked correctly; this one leftover degree-offreedom/constraint is used to fix the remaining “undetermined factor” mentioned at the end
of the last section.

13

2.6

Equivalent Representations of Collective Hamiltonian

We have shown in Sec. 2.4 that in each even order in the collective Hamiltonian the GDM
method gives one constraint (in the harmonic order it is the RPA secular equation). In this
section we show that these constraints fix the collective Hamiltonian (2.4) completely. In Eq.
(2.4) all terms with the right symmetry are kept thus the expansion is complete. However, in
fact it is over-complete. Different expansions are equivalent if they are related by canonical
transformations of collective variables α and π. The number of independent parameters of
H should be equal to the number of possible combinations as in Eq. (2.4) minus the number
of allowed transformations.
Microscopic estimates of quite general type [15] show that Λ(mn) ∼ Ω−(m+n−2)/2 , where
Ω is the collectivity factor, the effective number of simple quasiparticle excitations contributing to the collective mode. The solvable Lipkin model and quadrupole plus pairing model
confirm these estimates [12]. In the case of strong collectivity, Ω

1.

Let us count the number of transformations (α, π) → (¯ , π ),
α ¯

x(mn)

α=
m≥0,n≥0

1 {¯ m , π n }
α ¯
,
2 m! n!

αm ¯ n
(mn) 1 {¯ , π } ,
y

π=

2

m≥0,n≥0

(2.34)

m! n!

which preserve the commutator algebra,

[α, π] = [¯ , π ] = i.
α ¯

14

(2.35)

The transformation parameters x(mn) , y (mn) being of the order Ω−(m+n−1)/2 will not change
the dependence of Λ(mn) in Eq. (2.4) on Ω. The constant terms, x(00) and y (00) , corresponding to a trivial translation of origin are not needed as this choice is already made by
selecting Λ(10) = Λ(01) = 0 in the collective Hamiltonian (2.4); in the case of a multipole
collective mode such terms would violate rotational symmetry. In the linear terms we can
set x(10) = y (01) = 1, which would correspond to a rescaling of α or π [Λ(02) = 1 in Eq.
(2.4)]. The parameters x(mn) and y (mn) vanish for odd and even n, respectively, because of
the wrong time-reversal symmetry.
Using Eqs. (2.34) and (2.35) we have

[α, π] =

[m(s − n) − n(r − m)]
i
x(mn) y (r−m,s−n) ·
{¯ r−1 , π s−1 },
α
¯
2 rsmn
m! n! (r − m)! (s − n)!

(2.36)

where in the coefficient of {¯ r−1 , π s−1 } we keep only the leading terms in 1/Ω, that is,
α
¯
terms ∼ Ω−(r+s−2)/2 . The summation runs over r ≥ m and s ≥ n, and as seen from the
numerator, r ≥ 1 and s ≥ 1. In addition, s is odd, otherwise x(mn) y (r−m,s−n) vanishes. The
starting term, r = s = 1, gives correctly i. The terms with r + s ≥ 3 and an odd s should
vanish,

x(mn) y (r−m,s−n)

0=
mn

m(s − n) − n(r − m)
.
m! n! (r − m)! (s − n)!

(2.37)

These relations constrain x(mn) and y (mn) in the transformations (2.34).
Let us identify the independent parameters in the collective Hamiltonian (2.4), removing
the redundant degrees of freedom related to the transformations (2.34). In the quadratic
2

order, the transformations (2.34) do not change the harmonic terms ω α2 + 1 π 2 ; thus,
2
2

15

there is one independent parameter ω 2 . In the cubic order, the transformations (2.34)
with nonzero x(20) , x(02) , and y (11) influence the Hamiltonian parameters Λ(30) and Λ(12)
through the harmonic terms; there is one constraint (2.37) with (rs) = (21). Thus, the
renormalization of the collective variables removes the skew terms: Λ(30) and Λ(12) can be
set to zero and there remains no independent parameter in this order. In the quartic order,
the transformations (2.34) with nonzero x(30) , x(12) , y (21) , and y (03) influence Λ(40) , Λ(22) ,
and Λ(04) ; and there are two constraints (2.37) with (rs) = (31) and (13). Thus, there is
one independent parameter; we can, for example, choose it to be Λ(40) , and set Λ(22) and
Λ(04) to zero.
This process continues to anharmonic terms of higher orders. There is one independent parameter in each even order (we can choose it to be Λ(n0) excluding all momentumdependent high-order terms), and there are no independent parameters in odd orders. In
summary, the independent parameters in the collective Hamiltonian (2.4) can be identified
in the following form:
1
H = π 2 + V (α2 ) ,
2
α2
α4
α6
α8
V (α2 ) = ω 2
+ Λ(40)
+ Λ(60)
+ Λ(80)
+ ....
2
4!
6!
8!

(2.38)

We have shown in Sec. 2.4 that the GDM method gives one constraint in each even order
of anharmonicity, thus fixes all the independent parameters in Eq. (2.38). In this sense the
collective Hamiltonian is completely determined.
In practical applications, Eq. (2.38) may not be the most convenient choice for the
independent parameters of the collective Hamiltonian, which means solving the equations of
motion in the GDM method to infinitely high orders. Alternatively, we can pick up a certain

16

number (labeled N ) of terms in Eq. (2.4), putting other terms to zero; in other words, we
assume that the original fermionic Hamiltonian (2.1) can be sufficiently accurately mapped
onto a collective Hamiltonian with these N terms. Then in the GDM method we need to solve
the equations of motion up to the (2N )th order, in order to get N constraints. The quality
of the assumption of mapping can be checked self-consistently within the GDM method:
if the assumption is good, constraints from the orders higher than 2N should be satisfied
automatically. For the realistic nucleonic Hamiltonian, mapping onto a bosonic Hamiltonian
is guaranteed by the success of numerous phenomenological studies. The mapped quadrupole
phonon αµ is not necessarily the RPA phonon that is “proportional” to the real quadrupole
operator Qµ ; rather αµ is such a renormalized operator (2.34) that the mapping onto a given
form (the selected N terms) of the bosonic Hamiltonian is the “best”. The possibly infinite
series of the bosonic Hamiltonian expanded in the RPA phonon is “pushed”/resummed
into the selected finite-N terms by the renormalization (2.34). The expansion of the real
quadrupole operator Qµ in terms of αµ and πµ is obtained by substituting the solution (2.12)
into Qµ = Tr[qµ R] as shown in the next section.

2.7

Transition Rates

Let us consider the calculation of the transition rates for a one-body operator

†

Q=

q12 a1 a2 = Tr[qR].
12

17

(2.39)

The boson expansion of Q is calculated by substituting the solution (2.12) into Eq. (2.39),

Tr[qr(mn) ]

Q
m≥0,n≥0

1 {αm , π n }
≡
2 m! n!

Q(mn)
m≥0,n≥0

1 {αm , π n }
.
2 m! n!

(2.40)

The transition rates of Eq. (2.40) between eigenstates of the bosonic Hamiltonian (2.11)
should reproduce that of Eq. (2.39) between eigenstates of the nucleonic Hamiltonian (2.1).

18

Chapter 3
Validity of the Formalism
In this chapter we examine the validity of the GDM formalism by comparing its results with
that of the exact shell-model diagonalization. First we consider the well-known Lipkin model
that is analytically solvable. Next we consider models with factorizable forces numerically.
As will be shown, the GDM method reproduces the exact results quite well, for both energies
and transition rates, through the whole parameter space going from vibrational, γ-unstable
to deformed regions.

3.1

Lipkin Model

The Lipkin model [16] is widely used as a test ground for methods of collective motion. In
this model, there are two single-particle levels with energies ±1/2 (the spacing is the energy
unit), each with degeneracy Ω. The model Hamiltonian contains only “vertical” transitions

19

(σ = ±1; l = 1, 2, ..., Ω):

H=
σ,l

κ
σ †
aσ,l aσ,l +
2
2

†

aσ,l a

†
a
a
.
σ,l −σ,l −σ,l

(3.1)

σ,l,l

There are Ω particles; thus without interaction the lower levels (σ = −1) are completely
filled and the upper levels (σ = 1) are empty.
We follow Ref. [17] to get the analytical solution. It is easy to show that the quasi-spin
operators,

†

†

J+ = J− = Jx + iJy =

a+1,l a−1,l ,

Jz =

l

1
2

†

σaσ,l aσ,l ,

(3.2)

σ,l

satisfy the angular momentum algebra. Using Eq. (3.2) the Hamiltonian (3.1) is written as
1
2
2
H = Jz + κ(J+ + J− ).
2

(3.3)

Thus the total quasi-spin J is a good quantum number. With the Holstein-Primakoff transformation,

†

J+ = J− = −i A†

2J − A† A,

Jz = −J + A† A,

(3.4)

where A† and A are bosonic creation and annihilation operators with commutation relation
[A, A† ] = 1, the Hamiltonian (3.3) is written as an expansion over A† and A; or collective
coordinate α and momentum π by the canonical transformation
1
1
A = √ (uα − ivπ), A† = √ (uα + ivπ), uv = −1.
2
2

20

(3.5)

Assuming strong collectivity J

1, we keep only the leading order in 1/J. Under the

choice

u≈

√

1 + 2κJ,

1
v=− ,
u

(3.6)

the Hamiltonian becomes

H=

ω 2 2 1 2 Λ(40) 4 Λ(04) 4
α + π +
α +
π ,
2
2
4
4

(3.7)

with

ω 2 ≈ 1 − 4κ2 J 2 ,

Λ(40) ≈ κu4 ,

Λ(04) ≈ −κv 4 .

(3.8)

Other Λ(mn) vanishes in their leading order of 1/J.
For the ground-state band, the quantum number J is found from Eq. (3.2):

J = |Jz |max =

Ω
.
2

(3.9)

Using Eqs. (3.6) and (3.9) the result (3.8) can be written as

ω 2 ≈ 1 − (κΩ)2 ,

Λ(40) ≈ 6κ(1 + κΩ)2 ,

Λ(04) ≈

−6κ
.
(1 + κΩ)2

(3.10)

Equations (3.10) are accurate in the leading order of 1/Ω.

Now we apply the GDM formalism to the Hamiltonian (3.1). By going up to the sixth

21

order in the equation of motion, we get three constraints to fix ω 2 , Λ(40) , and Λ(04) ; they
agree with the exact results (3.10). Constraints from the eighth order and higher are satisfied
identically. Hence, the GDM method solves the Lipkin model completely in the leading order
of 1/Ω.

In order to compare directly the GDM results with the exact ones, we perform a numerical
example at Ω = 30. The results are shown in Fig. 3.1. As one can see, the first few excited
states of the anharmonic Hamiltonian agree very well with the exact results, while the RPA
fails very soon as κ increases to the critical point (RPA ω 2 = 0). For realistic medium
and heavy even-even nuclei, only in the vicinity of magic numbers the low-lying collective
excitations can be sufficiently described by the QRPA; in very many cases, including the softspherical, γ-unstable, and rotational dynamics, the collective modes lie near or beyond the
critical point of QRPA, so that the higher-order anharmonicities in the collective Hamiltonian
are indispensable.

In Fig. 3.1, as conventionally done, the bosonic Hamiltonian (3.10) is diagonalized in the
infinite phonon space {|0 ≤ n < +∞ } (|n is the state with n phonons, A† A|n = n|n ),
dropping the “divergent” Λ(04) < 0 term. However, as discussed in Ref. [12], the Hamiltonian
(3.10) should really be diagonalized in the finite physical space {|0 ≤ n ≤ Ω }. Acting A†
more than Ω times on the ground state runs out of valence particles thus gives zero. Similar
phenomena should exist in realistic nuclei, when the relevant particles (holes) near the Fermi
surface are exhausted by repeated actions of the phonon creation operator A† . In the GDM
method we need to determine the parameter u in A† (3.5). It is fixed by minimizing A|Φ
in its one-particle-one-hole components, where |Φ is the Hartree-Fock ground state that is

22

7
RPA
anharmonic
exact

6

En − E0

5
4
3
2
1
0

0

0.02

0.04

κ

0.06

0.08

0.1

Figure 3.1: The excitation energies En − E0 of the first five excited states as a function
of κ in the model (3.1) with Ω = 30. The red dashed-dotted line is the RPA calculation
2

1
(diagonalizing the harmonic Hamiltonian ω α2 + 2 π 2 ). The blue dashed line is obtained
2
2
(40)
1
by diagonalizing ω α2 + 2 π 2 + Λ 4! α4 in the infinite phonon space {|0 ≤ n < +∞ },
2
where ω 2 and Λ(40) are given by Eq. (3.10). The black solid line shows the exact results
from diagonalizing the original fermionic Hamiltonian (3.1) directly. For interpretation of
the references to color in this and all other figures, the reader is referred to the electronic
version of this dissertation.

23

0.04

0.2
1
2
3
4

(a)

0.02

E2 error

E1 error

0.06

0

−0.02
0

0.05

κ

(c)
E4 error

E3 error

0
−0.1
0.05

0.1

0.05

0.1

κ

0.4

0
−0.1
−0.2
0

(b)

−0.2
0

0.1

0.2
0.1

0.1

0.05

κ

0.2
0
−0.2
−0.4
0

0.1

(d)

κ

Figure 3.2: The errors of the first four excitation energies E − Eexact as a function of
κ in a set of calculations. Panels (a), (b), (c), and (d) are for the first, second, third, and
fourth excitation energy, respectively. Four lines on each panel are obtained by diagonalizing
different collective Hamiltonians in different phonon spaces. 1 (green dashed line): H =
ω 2 α2 + 1 π 2 + Λ(40) α4 in {|0 ≤ n < +∞ }. 2 (red dotted line): H = ω 2 α2 + 1 π 2 + Λ(40) α4
2
2
2
2
4!
4!
2 2
(40) 4
in {|0 ≤ n ≤ Ω }. 3 (blue dashed-dotted line): H = ω α + 1 π 2 + Λ 4! α in {|0 ≤
2
2
n ≤ Ω }, but with ω 2 = 1 − κ2 Ω(Ω + 2) replacing that in Eq. (3.10). 4 (black solid line):
2
(40)
(04)
H = ω α2 + 1 π 2 + Λ 4! α4 + Λ 4! π 4 in {|0 ≤ n ≤ Ω }, also with ω 2 = 1 − κ2 Ω(Ω + 2). 1
2
2
and 2 closely overlap and are indistinguishable in the figure.

24

mapped onto the bosonic state |n = 0 . The result is

u4 =

(10) 2
2<F <1 |r12 |
(01) 2
2<F <1 |r12 |

(3.11)

where the summation indices 1 and 2 run over unoccupied and occupied single-particle levels,
respectively (“F ” means Fermi surface). In the Lipkin model Eq. (3.11) reproduces the exact
result Eq. (3.6). Within the finite physical space, the Λ(04) < 0 term does not generate
divergences and should be kept. In order to identify the errors of the “anharmonic” curve in
Fig. 3.1, we plot the errors of the excitation energies for the first four excited states in a set
of calculations in Fig. 3.2. The overlap of curves 1 and 2 means that convergence is reached
for the first few excitation energies in the finite physical space {|0 ≤ n ≤ Ω }. Going from
curve 2 to curve 3, we remove the error owing to the inaccuracy of the harmonic potential
ω 2 in the next-to-leading order in 1/Ω: We replace ω 2 in Eq. (3.10) by ω 2 = 1 − κ2 Ω(Ω + 2),
which is correct not only in the leading order but also in the next-to-leading order of 1/Ω.
Finally in curve 4 the “divergent” term Λ(04) < 0 is included. We see that 4 is a much
better calculation than 1. The little “kink” on curve 4 near κ = 0.05 coincides with the
phase transition of the system, where the spectrum becomes doubly degenerate inside a
well-developed (large enough ω 2 < 0) double-well potential (see Fig. 3.1).
In summary, the microscopically calculated “divergent” terms should always be kept in
the bosonic Hamiltonian when diagonalizing. Within the finite physical space, they mix
nearby states without causing divergences. The low-lying collective states should compose
mainly of the states with the smallest phonon numbers, therefore be insensitive to small
variations (u in Eq. (3.11) and the boundary n ≤ Ω) of the physical phonon space. This is
the case in models of this thesis.

25

Figure 3.3: Single-particle level scheme. ˜ is the time-reversal level of 1. Each level has a
1
quantum number m.
From Fig. 3.2 we see that the next-to-leading order terms in 1/Ω of the RPA frequency
ω 2 could be important. This is the case also in models of the next section. In realistic nuclei,
the critical point ω 2 ≈ 0 could be reached at a relatively small Ω. For example, 100 Pd at
46
the critical point [18] has only eight valence particles (although pairing increase collectivity).
Hence an achieved improvement would be calculating ω 2 in its next-to-leading order of 1/Ω.

3.2

Models with Factorizable Force

In this section we consider models with the factorizable force that resemble the widely
used quadrupole-plus-pairing model for realistic nuclei. The single-particle space is drawn
schematically in Fig. 3.3. There are two groups of degenerate single-particle levels. The
Fermi surface is in between, thus without interaction the lower levels are completely filled
and upper levels are empty. Each single-particle level has a quantum number m that is
a half-integer. Degenerate time-reversal pairs have m of opposite sign, m˜ = −m1 . For
1
˜
fermions, |˜ = −|1 , and we choose the phases such that
1

|m = | − m

,

|−m = −|m

26

(m > 0).

The microscopic fermionic Hamiltonian (2.1) is given in normal-ordering form,

†

H=

e1 a1 a1 +
1

1
4

† †

V1234 N [a1 a2 a3 a4 ],
1234

† †

where N [a1 a2 a3 a4 ] is the normal-ordering form of operators with respect to the HF ground
state. The single-particle energies e1 = ±1/2 for the upper and lower levels, respectively. The
density matrix without interaction is ρ12 = δ12 n1 , where the occupation number n1 = 1(0)
for the lower(upper) single-particle levels. The residual interaction is of multipole-multipole
type (factorizable),
V1234 = −κ(q14 q23 − q13 q24 ),
where the multipole operator Q =

†
12 q12 a1 a2 is Hermitian and time-even. For simplicity

we assume that q is real, thus
q12 = q21 = q˜˜ = q˜˜ .
21
12
Operator q has certain selection rules with respect to the quantum number m, which will be
specified later. We further set diagonal matrix elements of q to zero, q11 = 0; hence in the
mean field, Q(00) = Tr[qρ] = 0.
In the following we consider four models with different structures (different configurations
of single-particle levels and different selection rules of q).
Model 1. We start with the simplest case. Both the upper and lower group have 12
degenerate single-particle levels with quantum numbers m = ± 1 , ± 3 , ..., ± 11 . Operator q
2
2
2
has the selection rule ∆m = 0, i.e., q12 vanishes unless m1 = m2 . The non-vanishing matrix
elements q12 (m1 = m2 ) are set to be 1.
In this model we find by numerical computation an additional “symmetry”, namely,

27

in the Hamiltonian (2.4) there are only three non-vanishing terms : ω 2 , Λ(40) , and Λ(22)
(besides π 2 /2). This is similar to the “quasi-angular-momentum symmetry” in the Lipkin
model, where the only three non-vanishing terms are ω 2 , Λ(40) , and Λ(04) (see Sec. 3.1).
The results for this model are shown in Fig. 3.4. The GDM calculation reproduces the
exact results of the shell model quite well, both for energies and transition rates. In the shell
model we calculated the lowest several states by the Lanczos method. The dashed line in
the upper panel is the beginning of the single-particle continuum (single-particle excitations
with high level density), only those collective states below the continuum were calculated
(due to computation time). In the GDM calculation the resulting bosonic Hamiltonian is
diagonalized in the finite “physical” bosonic space {|0 ≤ n ≤ 12 } (12 is the number of
fermions). The coefficient u in Eq. (3.5) is fixed by Eq. (3.11). In models of this section, u
is a number close to 1. The shown GDM energies and transitions are practically independent
of small variations of u around the above value.
As κ increases, the system goes from vibrational to γ-unstable region. In the vibrational
region, higher excited states are influenced more by the anharmonicities, as expected. At
large κ the spectrum becomes doubly degenerate in a deep double-well potential (large
negative ω 2 ). This is the analog of γ-instability of realistic nuclei in three dimensions.
An important point is that the GDM method works better with increasing collectivity Ω,
the number of effective particle-hole excitations contributing to the collective mode. Another
calculation has been done (not shown here) with 8 particles in 16 single-particle levels. The
GDM results of the current calculation (12 particles in 24 single-particle levels) have very
clear improvement over that of the former. In other words, the error in Fig. 3.4 may be of
order 1/Ω. The largest part of this error may come from the RPA frequency ω 2 . At the

28

3
RPA
exact
GDM

2.5

n

E −E

0

2
1.5
1
0.5
0

0

0.04

RPA
exact
Q01
Q12
Q03
Q23

12
10
〈m|Q|n〉

0.02

8

κ

0.06

0.08

0.1

0.06

0.08

0.1

6
4
2
0

0

0.02

0.04

κ

Figure 3.4: Excitation energies En − E0 and transition matrix elements m|Q|n in model 1
as a function of κ. The black solid lines are exact results of shell-model diagonalization. The
black dashed line is the beginning of the single-particle continuum. The red dashed-dotted
lines are the RPA results. The blue symbols are the GDM results. The stars are energies; the
squares, circles, up-triangles and down-triangles are matrix elements of Q between different
states. m|Q|n that are not shown vanish in both the shell model and the GDM calculations.

29

current stage, the GDM method calculates all Λ(mn) in their leading order of 1/Ω but not
the next. If the correct ω 2 = Λ(20) was smaller by a 1/Ω term than the one determined here
by the RPA equation, all the GDM curves would be shifted to the left (smaller κ), which
would decrease greatly the systematic error (see Fig. 3.4). This systematic error owing to
inaccurate ω 2 was present in all models in this work (see Figs. 3.5-3.7); Also, it is confirmed
in the Lipkin model where everything is known analytically (see Sec. 3.1). Hence an achieved
improvement would be calculating ω 2 in its next-to-leading order of 1/Ω.
Model 2. This model has the same single-particle configuration but the operator q has
now the selection rule ∆m = 0, ±1. Non-vanishing q12 are still set to be 1. Here we do not
find a additional symmetry as in model 1, so the problem exists of what should be the “best”
mapping. In the following we did three sets of GDM calculations. The first calculation is
possibly the simplest, which keeps only Λ(40) (besides ω 2 α2 /2 and π 2 /2) in Hb , fixed by the
constraint from the 4th order in the equation of motion. The second calculation keeps the
lowest two potential (no π dependence) terms Λ(40) and Λ(60) , which are fixed by the two
constraints from up to the 6th order in the equation of motion. The third calculation keeps
all quartic anharmonicities, Λ(40) , Λ(22) , and Λ(04) , fixed by the three constraints from up
to the 8th order in the equation of motion.
We first notice in Fig. 3.5 that in this model the single-particle continuum goes down
with increasing κ, as opposed to going up in model 1. This is because now mixing of singleparticle levels within the upper(lower) group is allowed by the selection rule that ∆m can be
±1. As a result, originally degenerate levels from the upper(lower) group get a finite spread,
which reduces the gap of the single-particle continuum. Only the first excited state is within
the gap and calculated in the shell model.

30

2

En − E0

1.5
1
0.5
0

0

0.005

0.015

0.02
κ

0.025

0.03

0.035

0.04

0.015

0.02
κ

0.025

0.03

0.035

0.04

RPA
exact

20

〈0|Q|1〉

0.01

Λ(40)
Λ(40,60)

15

Λ(40,22,04)
10
5
0

0

0.005

0.01

Figure 3.5: Excitation energies En − E0 and transition matrix elements 0|Q|1 in model 2
as a function of κ. The black solid lines are exact results of shell-model diagonalization. The
black dashed line is the beginning of the single-particle continuum. The red dashed-dotted
lines are the RPA results. The blue symbols (stars, up-triangles and circles) are results of
three different sets of GDM calculations, as labeled in the legend.

31

In the GDM calculations we see that the simplest one degree-of-freedom (Λ(40) ) calculation is reasonably well in most cases except at very large κ. The other two calculations
(Λ(40/60) and Λ(40/22/04) ) give essentially the same results (for the quantities shown in Fig.
3.5), although their common parameter, Λ(40) , is different. This insensitivity of GDM results to the degrees of freedom chosen, is important. As we said in Sec. 2.6, two different
bosonic Hamiltonians could be equivalent if they were related by canonical transformations/renormalizations (2.34) of variables α and π. This insensitivity simply says that the
GDM formalism knows these renormalizations and does them correctly. In model 1 we also
find this insensitivity (not shown). Finally we notice that in regions of ω 2 ∼ 1/Ω, calculations that go to higher orders in the equation of motion may give unphysical results. This
is again because the equation of motion are accurate only in the leading order of 1/Ω. The
fact that this “divergence” appears slightly before the instability point of RPA shown in Fig.
3.5, indicates again that the correct ω 2 may be smaller than the one calculated by RPA.
Models 3 and 4. At last we consider two models with single-particle configurations that
are asymmetric in upper and lower groups, which generates odd anharmonicities that are
necessary for deformation. In model 3, the lower group has 10 single-particle levels with
m = ± 3 , ..., ± 11 , the upper group has 14 single-particle levels with m = ± 1 , ..., ± 13 . In
2
2
2
2
model 4, the lower group has 12 single-particle levels with m = ± 1 , ..., ± 11 , the upper group
2
2
has 10 single-particle levels with m = ± 1 , ..., ± 9 . In both models, operator q still has the
2
2
selection rule of ∆m = 0, ±1, with non-vanishing matrix elements set to be 1. Model 3
has a slightly larger asymmetry than that of model 4, and their signs of the asymmetry are
different.
These two models are more complicated in the sense that now there are more active

32

degrees of freedom (odd anharmonicities). In the GDM method, we do a possibly simplest
calculation. We keep in Hb only Λ(30) , Λ(12) , and Λ(40) (besides ω 2 α2 /2 and π 2 /2). Λ(30)
and Λ(12) are fixed by requiring Q(20) = Q(02) = 0 in the solution (2.40). Then Λ(40)
is fixed by the constraint from the 4th order in the equation of motion. The requirement
Q(20) = Q(02) = 0 is the same as that for models 1 and 2 without upper-lower asymmetry,
by which Λ(30) and Λ(12) vanish.

The results are shown in Figs. 3.6-3.7. The deformation begins around the critical point
of RPA when ω 2 becomes negative. In the vibrational region the potential is stiff and deformation is not easy. As κ increases, the potential becomes flat in bottom and finally of
a double-well shape. Then, even a relatively small odd anharmonicity (here mainly Λ(30) )
can tilt the potential and generate large deformation. We notice firstly that the GDM calculations give the correct sign of deformations: positive/negative for the ground/first-excited
state of model 3, and vice versa for model 4. In the realistic situation of three dimensions,
Λ(30) ((ˆ × α)2 × α)0 ∼ Λ(30) β 3 cos 3γ (ˆ is the quadrupole phonon and β, γ are Bohr
α
ˆ
ˆ
α
shape variables), the sign of Λ(30) “determines” the intrinsic shape of the nucleus (prolate or
oblate). This is especially interesting in the transitional regions where the rotor formula is
not applicable. Secondly, the quantitative agreement of deformation is also good except at
the largest κ. There the deformation “saturates” towards its maximal possible value within
the model space, favored by energy. Meanwhile in the boson mapping, we are too close to
the boundary of the finite physical bosonic space, and the GDM results become inaccurate.
In realistic nuclei this “saturation” may not happen. The number of participating/active
nucleons is usually around 30 in well-deformed nuclei, which is much larger than that around
10 in the current models. Finally, we would like to point out that the first excited state in

33

1

RPA
exact
GDM

En − E0

0.8
0.6
0.4
0.2
0

0

RPA
exact
Q00
Q11
Q01

20
10
〈m|Q|n〉

0.01

0.02
κ

0.03

0.04

0.02
κ

0.03

0.04

0
−10
−20

0

0.01

Figure 3.6: Excitation energies En − E0 and transition matrix elements m|Q|n in model 3
as a function of κ. The black solid lines are exact results of shell-model diagonalization. The
black dashed line is the beginning of the single-particle continuum. The red dashed-dotted
lines are the RPA results. The blue symbols are the GDM results. The stars are energies;
the squares, circles and up-triangles are matrix elements of Q between different states.

34

1
RPA
exact
GDM

En − E0

0.8
0.6
0.4
0.2
0

0

0.005

0.015

0.02
κ

0.025

0.03

0.035

0.04

0.015

0.02
κ

0.025

0.03

0.035

0.04

RPA
exact
Q00
Q11
Q01

20
10
〈m|Q|n〉

0.01

0
−10
−20

0

0.005

0.01

Figure 3.7: Excitation energies En − E0 and transition matrix elements m|Q|n in model 4
as a function of κ. The black solid lines are exact results of shell-model diagonalization. The
black dashed line is the beginning of the single-particle continuum. The red dashed-dotted
lines are the RPA results. The blue symbols are the GDM results. The stars are energies;
the squares, circles and up-triangles are matrix elements of Q between different states.

35

our simple models is not a “rotational” state, rather it corresponds to the next “band head”
in realistic rotational nuclei. The rotational states that are very low in energy come in only
in three dimensions.

To summarize, in this chapter we demonstrate the validity of the GDM procedure for microscopic calculation of the collective/bosonic Hamiltonian. The lowest several states of this
bosonic Hamiltonian quite well reproduce the collective states of the exact shell model, both
for energies and transition rates, in a wide range from vibrational, γ-unstable, to deformed
systems. Specifically, we show that deformation can be described without introducing a
deformed mean field. The traditional procedure of “symmetry breaking and restoration”,
first “statically” breaks rotational symmetry in the ground state, by representing the latter as a Slater determinant of deformed single-particle levels (Nilsson levels); then projects
afterwards to good angular momentum. However, in case of large shape fluctuations (flat
minimum of energy surface) or shape coexistence (two close minima), it may fail. On the
other hand, the GDM procedure always conserves the rotational symmetry. Deformations
are put in “dynamically” at higher orders (for example cubic terms) beyond the mean field.
Thus it is suitable to describe such phenomena as shape fluctuations and coexistence.
In realistic nuclei, the gap of the single-particle continuum is generated by pairing correlations. The GDM formalism based on the Hartree-Fock-Bogoliubov variational method is
straightforward as will be shown in the following chapters. However, another treatment may
be possible. Instead of introducing Bogoliubov quasi-particles and representing the ground
state as their vacuum, the pairing correlations are considered in higher orders beyond the
mean field, by keeping both the particle-hole and particle-particle channels in the factoriza-

36

† †

†

†

†

†

† †

tion a4 a3 a2 a1 ≈ a4 a1 · a3 a2 − a4 a2 · a3 a1 + a4 a3 · a2 a1 . In this way the exact particle number
is always conserved. Work along this line is in progress and results seem promising.

37

Chapter 4

Formalism for Realistic Nuclei

In this chapter we consider the GDM formalism for realistic nuclei. There are three complications. A realistic nucleus has two kinds of fermions; symmetries, e.g. rotational invariance,
need to be respected; and pairing correlations should be considered.

As in the BCS theory we substitute the original system by a grand-canonical ensemble,
in which the chemical potential is fixed by the average particle number of the ground state
in the mean-field order. In this case we need to consider the equation of motion of not
†

only a2 a1 but also a2 a1 . A good treatment of the superfluid ground state, on top of which
collective excitations are formed, is essential.

The collective mode operators αλµ , πλµ have quantum numbers (angular momentum λ
and its projection µ) corresponding to symmetries of the Hamiltonian. In this section we
keep only the quadrupole mode (λ = 2) which is the most important one at low energy. The
case of interacting modes (quadrupole and octupole) is discussed briefly in Appendix D.

38

4.1

Formalism

The microscopic fermionic Hamiltonian for the grand-canonical ensemble is still given by
ˆ
ˆ
Eq. (2.1): we include the −µN term in ˆ (µ is the chemical potential and N is the particle
number operator), and the single-particle index 1, 2 . . . can run over protons and neutrons.
Isospin may not be conserved for some effective interactions. We do not write V in the form
J
V(j j ),(j j ) ; 12 and V1234 carry all the symmetries of H implicitly.
1 2
3 4

Now the reference state |Φ does not have definite particle number,

†

Φ|a1 a2 |Φ ≡ ρ21 ,

Φ|a1 a2 |Φ ≡ κ21 ,

(4.1)

where κ is the pair correlator [19]. And we need two density matrix operators

†

R12 ≡ a2 a1 ,

K12 ≡ a2 a1 ,

(4.2)

and two self-consistent field operators

W {R}12 ≡

V1432 R34 ,

f {R} ≡ + W {R},

(4.3)

34

∆{K}12 ≡

1
2

V1234 K43 .

(4.4)

34

It will be convenient to introduce (RT , f T are transposed matrices)




K
 R

D{R, K} ≡ 

K † I − RT


,



∆{K} 
 f {R}
S{R, K} ≡ 
.
† {K} −f T {R}
∆

(4.5)

The expansion of density matrix operators replacing Eq. (2.12) is now over the multipole

39

mode operators αλµ and πλµ . The operator structures are more complicated owing to
angular-momentum vector coupling. Let us choose the basis for the expansion as {Oχ },
where χ combines all the necessary quantum numbers to specify a term in the basis, including
the power of coordinate a (αa ), the power of momentum b (π b ), the total angular momentum
L, and others. For example, Oχ = {(α × α)l , π}L has χ as a summary of a = 2, b = 1, L, l.
Hence the expansions of R, K replacing Eq. (2.12) are

R=
χ

[rχ Oχ ]0 ,
0

K=
χ

[kχ Oχ ]0 ,
0

(4.6)

where rχ , kχ are matrices in the single-particle space. The symbol [ ]0 means angular0
momentum vector coupling to a scalar. Substituting Eq. (4.6) into Eq. (4.5) we get the
expansion of the latter

D=
χ

[Dχ Oχ ]0 ,
0

S=
χ

[Sχ Oχ ]0 ,
0

(4.7)

where Dχ and Sχ have a dimension that is twice of that of rχ and kχ .

Similarly, in the expansion of the Hamiltonian only the scalar terms (L = 0) exist. We
choose the basis as {Oβ }. In principle the choice of β can be different from χ. The expansion
of H replacing Eq. (2.11) is

H=

Λβ Oβ

(4.8)

β

In addition to Eq. (2.7), we calculate the exact equation of motion [K12 , H]. On the

40

right-hand side of the resulting equations, we substitute the factorization

† †

a4 a3 a2 a1

†

†

†

†

† †

a4 a1 · a3 a2 − a4 a2 · a3 a1 + a4 a3 · a2 a1

(4.9)

replacing Eq. (2.8). In Eq. (4.9) both the particle-particle and particle-hole channels are
kept. The results can be summarized as

[D, H]

[S, D].

(4.10)

This is the equation of motion in the collective band replacing Eq. (2.9).

Let us introduce the group algebraic properties

[Oχ , Oβ ] = i

χβ λ Oλ ,

(4.11)

λ

where χβ λ vanishes unless Lλ = Lχ , aλ = aχ + aβ − 1, bλ = bχ + bβ − 1. And

(Oχ Oη )L =

Lλ =L

(χη) λ Oλ ,

(4.12)

λ

where (χη) λ vanishes unless aλ = aχ + aη , bλ = bχ + bη .

χβ λ and (χη) λ are

complicated basis-dependent (χ and β) expressions consisting of Clebsch-Gordan coefficients.
Substituting Eqs. (4.11) and (4.12) into Eq. (4.10) we get

Λβ χβ λ [Dχ Oλ ]0
0

i
χ
η

λ

(χη) λ [(Lχ Lχ )0(Lη Lη )0]0|[(Lχ Lη )Lλ (Lχ Lη )Lλ ]0 [[Sχ , Dη ]Lλ Oλ ]0 .
0

=
χ

β

λ

41

On the right-hand side, [(Lχ Lχ )0(Lη Lη )0]0|[(Lχ Lη )Lλ (Lχ Lη )Lλ ]0 is a 9j-symbol. Since
terms Oλ in the basis are linearly independent, we get one equation for each Oλ equating
the corresponding coefficients,

i

Λβ χβ λ Dχ
χ

β

(χη) λ [(Lχ Lχ )0(Lη Lη )0]0|[(Lχ Lη )Lλ (Lχ Lη )Lλ ]0 [Sχ , Dη ]Lλ .

=
χ

(4.13)

η

This is the final set of GDM equations replacing Eq. (2.14). The “kinematic” coefficients
for rotational symmetry are singled out. The “dynamics” from the Hamiltonian is hidden in
the self-consistent field Sχ (4.5). Equation (4.13) is valid for the general case of interacting
multiple collective modes (for an example see appendix D). In the next section we list
explicitly its lowest orders in the case of a single quadrupole mode.

4.2

Equations of Motion for Quadrupole Phonon
†

†

The multipole mode operators αλµ , πλµ carry quantum numbers of angular momentum λ,
†

†

its projection µ, and parity (−)λ . The coordinate αλµ is time-even, and the momentum πλµ
is time-odd. Their Hermitian properties are

†

αλµ = (−)λ−µ αλ−µ ,

†

πλµ = (−)λ−µ πλ−µ .

(4.14)

The commutation relation is given by

†

[αλµ , πλ µ ] = i δλλ δµµ .

42

(4.15)

Here we consider only the quadrupole mode λ = 2, and drop the label λ.

The basis {Oχ } is specified in the following expansion of the density matrix operators,

rµ
µ

+
+

1
6

1
2

1
(01) †
πµ +

(10) †
αµ +

R=ρ+

rµ
µ

(02)

L=0,2,4 µ

rLµ (π † × π † )L +
µ

(30)

L=0,2,3,4,6 µ

rLµ {(α† × α† )lL , α† }L +
µ
+
+

1
4

2

1
4

1
6

1
2

(20)

L=0,2,4 µ

rLµ (α† × α† )L
µ
(11)

L=0,1,2,3,4 µ

rLµ {α† , π † }L
µ

(03)

L=0,2,3,4,6 µ

rLµ {(π † × π † )lL , π † }L
µ
(21)

L=0,1,2,3,4,5,6 l=0,2,4 µ

rLlµ {(α† × α† )l , π † }L
µ

(12)

L=0,1,2,3,4,5,6 l=0,2,4 µ

rLlµ {α† , (π † × π † )l }L + . . . (4.16)
µ

Three identical d bosons can couple to L = 0, 2, 3, 4, 6. In the α3 and π 3 terms of Eq.
(4.16) we choose the intermediate quantum number for each L to be lL ; this choice will not
influence the final results. Similarly the expansion for the operator K is
1
(01) †
πµ +

(10) †
αµ +

K =κ+

kµ
µ

kµ
µ

2

(20)

L=0,2,4 µ

kLµ (α† × α† )L + . . .
µ

(4.17)

The basis {Oβ } is specified in the expansion of the Hamiltonian,
1√
ω2 √
5(α × α)0 +
5(π × π)0
H = E0 +
0
0
2
2
(12) √
Λ(30) √
2 , α}0 + Λ
+
5{(α × α)
5{α, (π × π)2 }0
0
0
6
4
Λ(04) √
Λ(40) √
+
5((α × α)0 × (α × α)0 )0 +
5((π × π)0 × (π × π)0 )0
0
0
4
4
(22)
ΛL √
+
5{(α × α)L , (π × π)L }0 + . . . .
0
8
L=0,2,4

43

(4.18)

H is Hermitian, time-even, invariant under rotation and inversion.

4.2.1

Hartree-Fock-Bogoliubov Equation

In the mean-field order, the constant terms of Eq. (4.13) give the HFB equation,

[S (00) , D(00) ] = 0.

(4.19)

Equation (4.19) says that S (00) and D(00) can be diagonalized simultaneously,
 





 

0  
  E 0   n
  = 0,
 
 , 
0 I −n
0 −E

(4.20)

where E and n are diagonal matrices. The chemical potential µ (buried in E) is determined
by N = Tr{ρ} =

1 ρ11 =

†
1 n1 . We introduce the quasiparticle operators bλ , bλ through

the unitary canonical transformation

†

†

(u∗ a1 − v1λ a1 ), bλ =
1λ

†

∗
(u1λ a1 − v1λ a1 ).

(4.21)

If the reference state |Φ is a “quasiparticle vacuum”, |Φ = norm ·

λ bλ |0 , then by Eq.

bλ =
1

1

(4.1) we get:

ρ = vv † ,

κ = −vuT .

44

(4.22)

In this case D(00) is diagonalized by the canonical transformation (4.21):




U =

u†

−v T

−v †

uT


,





 0 0 
U D(00) U † = 
,
0 I

(4.23)

where the matrix n in Eq. (4.20) vanishes. The HFB equation (4.19) requires that S (00)
be diagonalized by U simultaneously. In this work we assume that |Φ is a “quasiparticle
vacuum”. From now on we always work in the quasiparticle basis (multiplying Eq. (4.13)
by U from the left and U † from the right).

4.2.2

Quasi-Particle Random Phase Approximation

In the harmonic order, the linear terms of Eq. (4.13) give the QRPA equations

†

iD(10) = [S (00) , D(01) ] + [S (01) , D(00) ],

(4.24)

†

− iω 2 D(01) = [S (00) , D(10) ] + [S (10) , D(00) ].

(4.25)

πµ :
αµ :

This is a linear homogenous set of equations for D(10) and D(01) , a non-zero solution requires
a zero determinant, from which we solve for the QRPA frequency ω 2 .

4.2.3

Cubic Order

In the cubic order, Eq. (4.13) gives

(α† × α† )L /2, L = 0, 2, 4 :
µ
(11)

−2iω 2 DL

(20)

(20)

− 2iδL,2 Λ(30) D(01) = [S (00) , DL ] + [SL

45

, D(00) ] + 2 [S (10) , D(10) ]L , (4.26)

(π † × π † )L /2, L = 0, 2, 4 :
µ
(11)

2iDL

(02)

(02)

− iδL,2 Λ(12) D(01) = [S (00) , DL ] + [SL

, D(00) ] + 2 [S (01) , D(01) ]L ,

(4.27)

{α† , π † }L /2, L = 0, 1, 2, 3, 4 :
µ
(02)

−iδL,even ω 2 DL
(11)

(11)

= [S (00) , DL ] + [SL

4.2.4

(20)

+ iδL,even DL

+ iδL,2 Λ(12) D(10)

, D(00) ] + [S (10) , D(01) ]L − [D(10) , S (01) ]L .

(4.28)

Quartic Order

In the quartic order, Eq. (4.13) gives

{(α† × α† )lL , α† }L /6, L = 0, 2, 3, 4, 6 :
µ
−

3i 2
ω
2

(21)

l=0,2,4

(30)

(11)

L
DLl · γl,l − 3i(−)L Λ(30) DL
L
(30)

= [S (00) , DL ] + [SL

, D(00) ] +

3
2

(20)

([Sl

L
L=2
· γ2,l − 3iδL2 Λ(40) D(01) · γ0,l
L

L

(20)

, D(10) ]L − [Dl

L
, S (10) ]L ) · γl,l , (4.29)
L

l=0,2,4

{(α† × α† )l , π † }L /4; l = 0, 2, 4; L = 0, 1, 2, 3, 4, 5, 6 :
µ
−2i(−)L ω 2

D
l =0,2,4

2i
4 i (30) L
(12)
(30)
L
· gl,l +
δl,l DL +
D
· gl,l
L
Ll
L
3
3 L

(02)

−2i δl2 Λ(30) DL
(21)

(20)

+ 2 iΛ(12) DL

(21)

(20)

= [S (00) , DLl ] + [SLl , D(00) ] + [Sl
+2

([S

(22)

L
· gl,2 + i δL2 Λl

D(10)

(20)

, S (01) ]L

, D(01) ]L − [Dl

(11)
(11) (10) L
L
, D(10) ]L − [D
,S
] ) · gl,l ,
l
l

l =0,1,2,3,4

46

(4.30)

{α† , (π † × π † )l }L /4; l = 0, 2, 4; L = 0, 1, 2, 3, 4, 5, 6 :
µ
−

2i
4i
(03)
(03)
L
δl,l (−)L ω 2 DL − ω 2 DL (−)L−lL · gl,l + 2i(−)L
L
3 L
3
(11)

−iδl2 Λ(12) DL
(12)

(11)

+ 2i(−)L Λ(12) DL

(12)

(21)
L
· gl,l
Ll

D
l =0,2,4

(22)

L
· gl,2 − i δL2 Λl

D(01)

(02) L
(02)
] − [D(10) , Sl ]L

= [S (00) , DLl ] + [SLl , D(00) ] + [S (10) , Dl
+2

([S

(11)
(11) (01) L
L
, D(01) ]L − [D
,S
] )(−)L−l gl,l ,
l
l

(4.31)

l =0,1,2,3,4

{(π † × π † )lL , π † }L /6, L = 0, 2, 3, 4, 6 :
µ
3i
(−)L
2

l=0,2,4

(03)

(03)

= [S (00) , DL ] + [SL

3i
(12)
(02)
L
L
L=2
DLl · γl,l − Λ(12) DL · γ2,l + 3iδL2 Λ(04) D(10) · γ0,l
L
L
L
2
, D(00) ] +

3
2

(02)

([Sl

(02)

, D(01) ]L − [Dl

L
, S (01) ]L ) · γl,l . (4.32)
L

l=0,2,4

The numerical coefficients γ L and g L are defined by
l,l

l,l

L
L
{(α × α)l , α}L = γl,l · {(α × α)l , α}L , (γl,l = 1),
µ
µ
1
1
L
gl,l · {(α × α)l , π}L .
{{α, π}l , α}L =
µ
µ
8
4

(4.33)
(4.34)

l=0,2,4

Values of γ L and g L are given in Appendix A.
l,l

l,l

As shown in chapter 2, in the GDM method there is one constraint from each even order
in the equations of motion. Here in the quartic order the constraint is found as following.
(22)

Setting L = 2, keeping only (D/S)(30/21/12/03) and Λ(40) , Λl

47

, Λ(04) terms, fl × Eq.
L

(4.29) + 1 ω 2
2

l=0,2,4 fl × Eq. (4.31) gives

L=2:

1
2

− iω 2 ·

fl · DL
L

(30)

fl · SL
L

1
+ ω2
2

i
+3iΛ(40) D(01) · f0 + ω 2
2
1
2

(03)

L

l=0,2,4
(30)

= [ S (00) ,
+[

(21)

fl · DLl + ω 2 fl · DL
1
+ ω2
2

(12)

fl · DLl

]

l=0,2,4
(12)

, D(00) ]

fl · SLl
l=0,2,4
(22)

fl · Λl

D(01) + . . .

(4.35)

l=0,2,4

2
l=0,2,4 fl × Eq. (4.30) +ω flL × Eq. (4.32) gives

−

i
2

1
2

L

1
2

= [ S (00) ,
+[

(30)

fl · DL

i·

L=2:

1
+ ω2
2

(22)

fl · Λl

l=0,2,4

(21)

(03)

fl · DLl + ω 2 fl · DL
L

l=0,2,4

(03)

(21)

l=0,2,4

(12)

fl · DLl

fl · SLl + ω 2 fl · SL
L

]

, D(00) ]

D(10) − 3iω 2 Λ(04) D(10) · f0 + . . .

(4.36)

l=0,2,4

where fl is defined in Eq. (A.7). A solvability condition exists because the variable parts
(those including D(30/21/12/03) ) of Eqs. (4.35) and (4.36) have the same structure as the
QRPA equations (4.24) and (4.25).

A computer code is needed for general GDM calculations using realistic forces. This
would be a generalization of the existing GDM code for calculations in chapter 3, by including
angular-momentum vector coupling. In the next chapter we consider the special case of the
quadrupole-plus-pairing Hamiltonian, where approximate analytical solutions are possible.

48

Chapter 5
Quadrupole-plus-Pairing Model
In this chapter the GDM formalism is applied to the specific Hamiltonian of the quadrupoleplus-pairing model. As was understood long ago [20, 21], this model combines the most important nuclear collective phenomena in particle-particle (pairing) and particle-hole (quadrupole
mode) channels. Approximate analytical solutions exist, if we keep only the coherent part in
the interaction, as was usually done for this model. These analytical GDM results are given
in Sec. 5.1, and compared with the exact results of the shell model diagonalization in Sec.
5.2. Finally, in Sec. 5.3 we do realistic calculations for tin isotopes.
Let us first introduce the model. The multipole operator is

†

†

Qλµ {R} = T r{qλµ R} =

†

†

qλµ12 a1 a2 ,

12
†
qλµ = fλ (r) · iλ Yλµ (θ, φ),

(5.1)
(5.2)

where fλ (r) is real. The definition of Eq. (5.2) differs from the “usual” one in two aspects:
†
†
a factor iλ is included, and qλµ ∼ Yλµ instead of qλµ , thus qλµ creates projection µ. The

49

Hermitian properties are

†

qλµ = (−)λ−µ qλµ ,

†

Qλµ = (−)λ−µ Qλµ .

(5.3)

Here we consider only the quadrupole mode, and drop the label λ = 2. The pairing operators
P and P † are defined by

P† =

1
2

† †
1

a1 a˜ ,
1

P =

1
2

1

a˜ a1 ,
1

(5.4)

˜
where ˜ is the time-reversed single-particle level of 1. For fermions, |˜ = −|1 .
1
1

The quadrupole-plus-pairing Hamiltonian is

H=
κ
+
4

1

G
†
( 1 − µ)a1 a1 −
4
†

1234 µ

12

†

† †
1

a1 a˜ a˜ a2
2
† †

(−qµ14 qµ23 + qµ13 qµ24 )a1 a2 a3 a4 .

Approximately, this Hamiltonian can be written as H ≈

†
1
†
1 ( 1 −µ)a1 a1 −GP P − 2 κ

(5.5)

†
µ Qµ Qµ .

The difference is in a one-body term originated from the Q · Q part. H is Hermitian and
time-even, which implies real G, κ, and 1 = ˜ . In a realistic nucleus there are protons and
1
neutrons; formally we can still use Eq. (5.5) if the quadrupole force strengths are the same
for proton-proton, neutron-neutron, and proton-neutron interactions (κp = κn = κpn = κ),
while remembering the pairing is treated for protons and neutrons separately (Gp = Gn ).
We will assume this is the case.

50

5.1

The GDM Method

In this section the GDM equations of motion of chapter 4 are solved in the quadrupole-pluspairing model. As was usually done, only the “coherent” part of the interaction is kept:
we neglect in the W {R} field (4.3) the Fock/exchange terms and the contribution from the
pairing force, and in the ∆{K} field (4.4) the contribution from the quadrupole-quadrupole
force. Then the results are analytical:

†

W {R}12 ≈ −κ

qµ12 Qµ ,

(5.6)

∆{K}12 ≈ δ1˜ GP.
2

(5.7)

µ

5.1.1

BCS Theory

In the quadrupole-plus-pairing model the HFB equation becomes the BCS equation. The
canonical transformation (4.21) mixes only time-reversal pair of orbitals,

u12 = δ12 u1 ,

v12 = −δ1˜ v1 ,
2

(u1 )2 + (v1 )2 = 1,

(5.8)

where u1 = u˜ , v1 = v˜ are real numbers. The density matrices (4.22) become
1
1

ρ12 = δ12 (v1 )2 ,

κ12 = δ1˜ u1 v1 .
2

(5.9)

The constant terms of Eq. (5.6) give

(00)

f {ρ}12 = δ12 [ 1 − µ] + w12 ≈ δ12 [ 1 − µ − G(v1 )2 ] ≡ δ12 e1 ,

51

(5.10)

where e1 ≡ 1 − µ − G(v1 )2 . And the constant terms of Eq. (5.7) give

∆{κ}12 ≈ δ1˜ GP (00) ≡ −δ1˜ ∆,
2
2

where the pairing energy ∆ ≡ −GP (00) = G
2

(5.11)

3 u3 v3 . Then the HFB equation (4.19) leads

to the BCS set of equations:

∆·

G
4

1
E1

e
1+ 1 ,
E1

(5.12)
(5.13)

E1 =
1
2

= 0,

e1 = 1 − µ − G(v1 )2 ,

(u1 )2 =

1−

(e1 )2 + (∆)2 ,

(5.14)

(v1 )2 =

1
2

e
1− 1 ,
E1

(5.15)

(v1 )2 .

(5.16)

1

N=
1

E1 is the quasiparticle energy. The chemical potential µ is fixed by Eq. (5.16). The gap
equation (5.12) has a non-trivial solution ∆ > 0 only if G is greater than its critical value
Gc [21]. For later convenience we introduce:

(2)

ξµ12 ≡

(u1 v2 + u2 v1 )
qµ12 , ηµ12 ≡ (u1 u2 − v1 v2 )qµ12 .
(E1 + E2 )2
(2)

(5.17)

The superscript “(2)” in ξµ12 stands for the denominator 1/(E1 + E2 )2 , other “powers” are
defined similarly.

52

5.1.2

Quasi-Particle Random Phase Approximation

Here we give only the final results. Detailed derivations can be found in Ref. [12]. The
QRPA secular equation determining frequency ω 2 is

(−)j2 −j1

1=κ
τ1 j1 τ2 j2

(2j1 + 1)(2j2 + 1)
E1 + E2
†
†
ξ 12 ξ 21 .
2 − ω2
5
(E1 + E2 )

(5.18)

†

ξ 12 ≡ τ1 j1 ξ † τ2 j2 is the reduced matrix element, the convention for which is given in
Appendix B. τ1 combines all other quantum numbers specifying a single-particle level, except
j1 . The normalization Q(10) is given by

1 = (κQ(10) )2

(−)j2 −j1
τ1 j1 τ2 j2

(2j1 + 1)(2j2 + 1)
(E1 + E2 )
†
†
ξ 12 ξ 21 . (5.19)
2 − ω 2 ]2
5
[(E1 + E2 )

Starting from the next (cubic) order, we restrict ourselves to regions near the critical
point of QRPA (ω 2 ≈ 0). Following the procedure in Appendix C, we give results of
Λ(30) and Λ(40) in their leading order of ω 2 (constant terms independent of ω 2 ). As in the
Landau phase transition theory, the critical point is such a position that the leading potential
term ω 2 α2 /2 happens to vanish. The next non-vanishing term, usually Λ(40) α4 /4, is thus
dominant restoring the stability of the system (see Refs. [15] and [12]).

53

5.1.3

Cubic Anharmonicity

In this model the cubic anharmonicity (C.7) becomes

.
Λ(30) = 3(κQ(10) )3

(2j1 + 1)(2j2 + 1)(2j3 + 1)
5
τ1 j1 τ2 j2 τ3 j3




 2 2 2 
(1)† (1)† †
· ξ 12 ξ 23 η 31 ,
·


 j j j 
1
2
3

(5.20)

.
where = means the equation is correct in constant terms but not in ω 2 terms or higher.

(40)

We give the expression of P (20) which will appear in Λc

(2j1 + 1)

[(u1 )2 − (v1 )2 ]2 .
]=
4E1

(−)j1 +j2

(2j1 + 1)(2j2 + 1)
5

P (20) · [ 1 − G
τ1 j1

−(κQ(10) )2
τ1 j1 τ2 j2

·

of next section:

(1)† (1)†

2u1 v1 · ξ 12 ξ 21 −

(u1 )2 − (v1 )2 † (1)†
· η 12 ξ 21 .
E1

(5.21)

P (20) is divergent when G is greater than but close to Gc . In this region of the pairing phase
transition, ∆ is small, and P (20) ∼ 1/∆. The GDM + BCS method is not valid in this
region: in the mean-field order the BCS solution already fails, as is well known.

54

5.1.4

Quartic Anharmonicity

The quartic anharmonicity Eq. (C.8) becomes

.
(5)†
f0 · Λ(40) = −2 f2 · (κQ(10) )2 (Λ(30) )2 Tr ξµ ξµ
− f2 · Λ(30) Λ(12)
(1)†

†

−f2 · Λ(30) (κQ(10) )3 Tr {ξ (1) , ξ (3) }L=2 ηµ + Tr {η, ξ (3) }L=2 ξµ
µ
µ
(3)†

(2)†

+ Tr {η, ξ (2) }L=2 ξµ
µ

+2 Tr {η, ξ (1) }L=2 ξµ
µ

(2)†

−f0 · G P (20) (κQ(10) )2

(u1 v1 + u2 v2 ) ξµ12 ξµ21
12

+f0 · G P (20) (κQ(10) )2 ·

{

(u1 )2 − (v1 )2

12
(10) )4
+(κQ

2E1

+

(u2 )2 − (v2 )2 (1) †
}ξµ12 ηµ21
2E2
(1)†

fl · Tr {ξ, (ξ (1) × ξ (1) )l }L=2 ξµ
µ

l=0,2,4
(1)†
†
− Tr {ξ (1) , {η, ξ (1) }l,(1) }L=2 ηµ
−Tr {η, {η, ξ (1) }l,(1) }L=2 ξµ
µ
µ

l,(1)

where {η, ξ (1) }12

,

(5.22)

= {η, ξ (1) }l /(E1 + E2 ). There is an undetermined parameter Λ(12) in
12

Eq. (5.22). Values of numerical factors fl are given in Appendix A. Equation (5.22) in
terms of reduced matrix elements can be found in Ref. [12].

5.2

Comparison with Exact Results

We compare the GDM results of Sec. 5.1 in the following model with the exact one of
NuShellX [22]. There are 10 fermions of one kind and four single-particle levels with energies:

single − particle levels
(MeV)

1p 1
2

7
3
5
0f 2 1p 2 0f 2

−0.1

0.0

55

1.0

1.1

4

excitation energy

3

exact E
2+
exact E
4+
GDM E
2+

2

1

0
0.0

0.1

0.2

0.3

G

Figure 5.1: Excitation energies (from Table 5.1) in the quadrupole plus pairing model as a
function of the pairing strength G, at the critical point ω = 0 (κ = κc ). The black squares
and red circles show the exact excitation energy of the first 2+ and 4+ state, respectively,
“NuShellX E2+ ” and “NuShellX E4+ ” in Table 5.1. The blue triangles give “GDM E2+ ”
from Table 5.1.
We take the radial wavefunctions to be harmonic oscillator ones. In Eq. (5.2) we take f (r)
†
†
ˆ ˆ
to be r2 so qµ = −ˆ2 Y2µ (θ, φ). For convenience we make qµ dimensionless by combining
r

its original dimension with κ (see the end of Appendix B). The model space is similar
1
to the realistic pf -shell, but the 1p 2 and 1p 3 levels are inverted to increase collectivity:
2

in the current case the q matrix elements (q 1 5 and q 3 7 ) are large between the
1p 2 ,0f 2
1p 2 ,0f 2
single-particle levels above and below the Fermi surface.
We restrict ourselves to the critical region (QRPA ω 2 ≈ 0). As said in Sec. 5.1.2, here
the stability of the system is restored by the quartic potential term Λ(40) α4 /4. We did a set

56

Table 5.1: Results of the quadrupole plus pairing model for different values of the pairing
strength G. All quantities are in units of MeV. ∆ is the solution of Eq. (5.12). The chemical
potential µ is the solution of Eq. (5.16). κc is the critical κ such that ω 2 in Eq. (5.18)
becomes zero. Λ(30) is given by Eq. (5.20). Λ(40) is given by Eq. (5.22) setting Λ(12) = 0.
“GDM E2+ ” is the excitation energy of the first 2+ state from diagonalizing Eq. (4.18) for
(22)

ω 2 = Λ(12) = Λ(04) = ΛL = 0. “NuShellX E2+ ” is the exact excitation energy of the first
2+ state from diagonalizing Eq. (5.5), in which G and κ are given by G and κc in the table.
Similarly, “NuShellX E4+ ” is the exact excitation energy of the first 4+ state.
G

0

0.03

0.06

0.09

0.11

0.12

∆

0.0

0.0

0.0

0.0

0.0

0.066

µ

0.5

0.5

0.5

0.5

0.5

0.454

κc

0.102

0.105

0.107

0.110

0.112

0.113

Λ(30)

-0.160 -0.173

-0.188

-0.203

-0.213

-0.219

Λ(40)

0.483

0.526

0.572

0.621

0.655

0.616

GDM E2+

0.882

0.908

0.933

0.959

0.976

0.955

NuShellX E2+

0.855

0.892

0.944

1.023

1.106

1.158

NuShellX E4+

0.778

0.827

0.927

1.110

1.284

1.383

G

0.15

0.18

0.21

0.25

0.30

∆

0.453

0.672

0.862

1.096

1.374

µ

0.444

0.429

0.415

0.395

0.370

κc

0.113

0.122

0.135

0.154

0.179

Λ(30)

-0.234 -0.270

-0.310

-0.378

-0.474

Λ(40)

1.185

1.918

2.901

4.683

7.830

GDM E2+

1.194

1.405

1.614

1.894

2.249

NuShellX E2+

1.353

1.552

1.764

2.059

2.438

NuShellX E4+

1.705

2.076

2.465

2.987

3.631

57

of calculations with increasing pairing strength G. At each value of G, the strength κ of the
Q · Q force is taken to be at the critical value κc such that the QRPA frequency ω 2 = 0. The
results are summarized in Table 5.1. For clarity, we draw the last three lines of Table 5.1
as Fig. 5.1. The coefficient Λ(40) in Table 5.1 is calculated by Eq. (5.22) setting Λ(12) = 0
(dropping the −f2 · Λ(30) Λ(12) term). A non-zero term Λ(12) in its reasonable range does
not influence Λ(40) much, since in the current model Λ(30) is small due to the approximate
symmetry with respect to the Fermi surface (see Table 5.1). Then “GDM E2+ ” is calculated
(22)

by diagonalizing Eq. (4.18), setting Λ(12) = Λ(04) = ΛL

= 0 (ω 2 = 0 since κ takes its

critical value).
The critical value of the pairing strength Gc is around 0.11 ∼ 0.12 MeV. When G < Gc ,
the BCS solution ∆ = 0, and µ can be anywhere between

7
0f 2

= 0 and

1p 3
2

= 1.0 MeV. We

checked that in this case our results (5.18-5.22) do not depend on the choice of µ. In Table
5.1 we fix µ at 0.5 MeV. In the region where G is greater than but close to Gc , our method
is invalid as discussed under Eq. (5.21). This is illustrated in Fig. 5.1 by the ‘kink’ on the
“GDM E2+ ” curve near G ∼ 0.12.
In Fig. 5.1 “exact E2+ ” and “exact E4+ ” are the exact results by NuShellX. At G = 0
the first excited state is 4+ instead of 2+ . In this case the 4+ state is a single-particle
7
excitation from 0f 2 to 0f 5 ; the 2+ state is a collective state with approximately half holes
2

in 1p 1 and 0f 7 levels, half particles in 1p 3 and 0f 5 . As G increases, the collective 2+ state
2
2
2
2
becomes the first excited state. When G is large enough, ∆ dominates over the original
single-particle spacing , and the results become stable. As an example, at G = 0.30 MeV,
the quasiparticle continuum starts at ∼ 3.5 MeV; from the second excited state 0+ at 3.506
MeV to 4.153 MeV there are 15 states with J P = 0+ , 2+ , 4+ , 6+ . The first excited state 2+

58

at 2.438 MeV should be identified as a collective state, stabilized at around sixty percents
within the gap.
It is seen in Fig. 5.1 that “GDM E2+ ” agrees well with the exact result “exact E2+ ” in
general. On the G < Gc side, our E2+ does increase with G although not rapidly enough.
On the G > Gc side, when ∆ is not too small, the agreement is very good.

59

Table 5.2: Neutron single-particle energies for 100 Sn and 132 Sn in units of MeV.
s.p. level
100 Sn
132 Sn

5.3

0g 9
2

0g 7
2

1d 5
2

1d 3
2

2s 1
2

0h 11
2

9
0h 2

-16.7893 -10.2893 -10.6089 -8.6944 -8.7167 -8.8152 -0.8152
-14.74

-9.74

-8.97

-7.31

-7.62

-7.38

-0.38

Tin isotopes

As a realistic application of the GDM formalism, we consider tin isotopes. They are semimagic nuclei; the proton number 50 is a magic number. In our calculations the protons are
7
kept inactive. For the neutrons, we keep 7 active single-particle levels: 0g 9 , 0g 2 , 1d 5 , 1d 3 ,
2
2
2

2s 1 , 0h 11 , and 0h 9 . Their energies for 100 Sn and 132 Sn are listed in Table 5.2. These energies
2
2
2
are taken from two shell-model calculations in Refs. [23] and [24]. For tin isotopes with mass
number from 102 to 130, the single-particle energies are given by linear extrapolation between
the above two ends 100 Sn and 132 Sn. We assume that their wavefunctions are of the harmonic
†

oscillator type, and the matrix elements of the quadrupole operator qµ = −r2 Y2µ are given
in appendix B. The results have the factor b2 , where the length parameter b =

/(mΩ0 )

characterize the size of the harmonic oscillator wavefunction, Ω0 is the frequency of the
harmonic mean field. In our calculation we take Ω0 = 45A−1/3 − 25A−2/3 MeV, suggested
in Ref. [25], fitting globally the radii of the nuclei. The variations of Ω0 are small along
the tin isotope chain, from 8.53 MeV for 100 Sn to 7.87 MeV for 132 Sn.
The strength G of the pairing interaction in the Hamiltonian (5.5) is fixed to be 0.175
MeV, reproducing the even-odd mass staggering. In Fig. 5.2 we plot the experimental evenodd mass staggering and the BCS pairing gap ∆. Near the neutron magic numbers (102 Sn,
130 Sn), the BCS theory incorrectly gives a too small solution, as is well known. Thus in the

following we only consider isotopes from 104 Sn to 128 Sn.

60

Sn isotope
1.5

∆ (MeV)

1

0.5

0
100

experiment
BCS

105

110

115
mass

120

125

130

Figure 5.2: Even-odd mass staggering in the tin isotopes. The black squares show the
experimental data MA − (MA−1 + MA+1 )/2. The red circles show the BCS pairing gap ∆
as the solution of Eq. (5.12).

61

The strength κ of the quadrupole interaction in the Hamiltonian (5.5) is fixed to be
κb4 = 0.13 MeV, where b is the above mentioned harmonic-oscillator length parameter. The
above value of κ reproduces in the QRPA calculation the experimental E(2+ ) at the two
ends (104 Sn and 128 Sn) of the tin isotope chain (see Fig. 5.3). In the middle, as the number
of valence particle increases, the QRPA frequency ω 2 becomes small and finally negative,
with increasing collectivity.
In the GDM calculation, we calculate Λ(30) and Λ(40) in their leading constant term in
the expansion over ω 2 , through Eqs. (5.20) and (5.22). The resulting Hamiltonian,
1√
ω2 √
5(α × α)0 +
5(π × π)0
0
0
2
2
(30)
(40) √
Λc √
2 , α}0 + Λc
5{(α × α)
5((α × α)0 × (α × α)0 )0 ,
+
0
0
6
4
H=

(5.23)

is diagonalized in the “physical bosonic space” {|0 ≤ n ≤ nmax }, defined in the same way
as that in chapter 3. Specifically, |n is the state with n phonons,

†
µ Aµ Aµ |n = n|n ,

†

where Aµ is defined by generalizations of Eqs. (3.5) and (3.11). The maximal number of
phonons, nmax , is given by the smaller of the number of pairs of valence particles or holes
in the 50 − 82 shell. For example, nmax = 5 for 110 Sn and nmax = 6 for 120 Sn.
In the Hamiltonian (5.23), Λ(30) and Λ(40) are calculated only in their leading order of ω 2 .
However, it should be a good approximation because here the small parameter ω 2 /(2∆)2 is
around 0.3. As explained in appendix C, anharmonicities Λ(mn) should be smooth functions
of ω 2 , in such a small range ∼ 0.3 their variations are usually small. This is the case for
the models in chapter 3. The Hamiltonian (5.23) is expected to work better with increasing
collectivity. Here in Fig. 5.3 we show the GDM results from 108 Sn to 124 Sn, where nmax ≥ 4.
Near magic numbers where nmax ≤ 3, it may be not meaningful to speak about quartic

62

Sn isotope

E(2+) (MeV)

2
1
experiment
QRPA
GDM

0
−1
−2
100

105

110

115
mass

120

125

130

125

130

0.3

0.2
0.15

B(E2;0 →2 )

+

0.25

+

2 2

(e b )

Sn isotope

experiment
QRPA
GDM

0.1
0.05
0
100

105

110

115
mass

120

Figure 5.3: Excitation energy E(2+ ) and transition rates B(E2; 0+ → 2+ ) for the first
excited state 2+ in the tin isotopes. The black squares show the experimental data taken
from Ref. [26]. The blue triangles are the GDM results. The red circles show the QRPA
results. On the upper panel, negative values for the QRPA E(2+ ) actually mean imaginary
frequency.

63

anharmonicity in the Hamiltonian (5.23).
In Fig. 5.3 we see that the GDM calculation with the Hamiltonian (5.23) reproduces the
experimental data quite well in the region 108 Sn ∼ 124 Sn, where the QRPA solution collapses
(negative ω 2 ). The quartic potential term Λ(40) α4 restores the stability of the system. In
the lower panel we calculate the transition rates (defined in appendix B) of the quadrupole
operator Qµ ≈ Q(10) αµ , with Q(10) given by Eq. (5.19). An effective charge eeff = 0.7 is
used in both the QRPA and the GDM calculations.
The key point here is that the residual interaction, κQ† Q and GP † P , should be smooth.
Consequently in the middle of the shell the QRPA solution collapses as a result of increasing
collectivity. A recent QRPA calculation with Skyrme interactions shows quite a similar
pattern [27]. In this region higher order anharmonicities, here mainly Λ(40) , restore the
stability of the system.

64

Chapter 6
Conclusion
This work develops the generalized density matrix (GDM) formalism for microscopic calculation of the collective/bosonic Hamiltonian. The standard method in nuclear structure,
shell model, requires impractical computation time for majority of medium and heavy nuclei. Thus, we have to look for physically justified approximations by identifying the most
important degrees of freedom out of the huge shell-model space. The collective modes, like
vibration or rotation of the nucleus as a whole, are definitely the correct choice for describing
the low-lying spectrum. These collective modes have the nature of bosons; phenomenological models using an effective bosonic Hamiltonian are often quite successful in explaining
the experimental data. However, the complete microscopic theory of deriving the collective/bosonic Hamiltonian from the underlying shell-model Hamiltonian is still missing after
several decades.
In this work we propose a procedure based on the generalized density matrix. The
procedure is rather simple, clean, and consistent. In compact form, there are only two
equations, (2.14) and (2.33). They fix the bosonic Hamiltonian completely, which greatly

65

generalizes the well-known random phase approximation that calculates only the harmonic
terms. Many nuclei, including those at the γ-unstable and rotational limit, have a small
or negative harmonic potential, thus the anharmonic terms are indispensable restoring the
stability of the system. Specifically, we show that deformation can be described without
breaking the rotational symmetry in the mean field (cranking); relevant correlations are
put in “dynamically” at higher orders. Thus the GDM approach is suitable to describe such
phenomena as shape fluctuations and coexistence, also in the transitional region the coupling
between rotational and vibrational motions.
In this thesis, the test of the proposed GDM method in its most general form has been
limited to systems without rotational symmetry (chapter 3). The agreement with the exact
results, both for energies and transition rates, is very good. The application to realistic
nuclei has only been done for the “analytical” quadrupole-plus-pairing model near the critical
region of the quasi-particle random phase approximation (chapter 5). The next step would
be generalizing the existing GDM code by including angular-momentum vector coupling that
are necessary for general realistic calculations.
As a systematic way to construct the effective Hamiltonian, the applicability of the GDM
formalism should be more broad. We list here some future directions. Firstly, it may be
possible to develop a theory for pairing without introducing Bogoliubov quasi-particles, which
improves the BCS or Hartree-Fock-Bogoliubov treatment in the sense that the exact particle
number is always conserved. The current available results seem promising. Secondly, the
GDM formalism could be used to construct the effective Hamiltonian for odd-mass nuclei,
where the coupling between the collective modes and the unpaired nucleon is very important
(see appendix E). Thirdly, a more ambitious project would be doing better GDM calculations

66

for the even-even nuclei by decoupling two or more nucleons from the collective modes. In
this way the “single-particle” states could be studied and collective states should be more
accurate. From the computational/mathematical point of view, the GDM formalism should
be identical to the shell model if all the nucleons are decoupled from the collective modes.
Practically, we could seek to do the best calculation allowed by computer time. The odd-mass
and odd-odd nuclei could be treated along the same line.
Of course, the vitality of the GDM method depends on the results. But we hope that
the current work could serve as a beginning of its promising future.

67

APPENDICES

68

Appendix A

L
L
Values of gl,l and γl,l

The definition of g L is given by Eq. (4.34),
l,l

1
{{α, π}l , α}L =
µ
8

l=0,2,4

1
L
gl,l · {(α × α)l , π}L ,
µ
4

which implies

L
gl,l = (−)l−l




(2l + 1)(2l + 1)

2 2

l






 2 L l 


The definition of γ L is given by Eq. (4.33),
l,l

L
{(α × α)l , α}L = γl,l · {(α × α)l , α}L .
µ
µ

69

L
= gl ,l .

(A.1)

Analytical expressions of γ L can be obtained in the following way. Assume l and l are
l,l

even. We have the identity

L=0,2,3,4,6
l =0,2,4

[ {(α† × α† )l , α† }L , (π × π)0 ] |
µ
0

2
L
( δl,l + 2 · gl,l ) · {(α† × α† )l , π † }L .
= i√
µ
5 l=0,2,4

(A.2)

Replacing l in Eq. (A.2) by l we obtain

L=0,2,3,4,6
l =0,2,4

[ {(α† × α† )l , α† }L , (π × π)0 ] |
µ
0

2
L
= i√
( δl,l + 2 · gl,l ) · {(α† × α† )l , π † }L .
µ
5 l=0,2,4

(A.3)

Let us take the ratio of Eq. (A.2)/Eq. (A.3). The left-hand side is γ L

l ,l

by Eq. (4.33).

Since in the right-hand side of Eq. (A.2) and Eq. (A.3), {(α† × α† )l , π † }L with different l
µ
are linearly independent, we have

γlL,l

=

δl,l + 2 · g L

l,l
,
δl,l + 2 · g L
l,l

(l = 0, 2, 4).

(A.4)

The ratio on the right-hand side is independent of l. Since the matrix δl,l + 2 · g L is
l,l
symmetric (with respect to l, l ), Eq. (A.4) implies

L
δl,l + 2 · gl,l = flL · flL ,

70

(A.5)

where flL =

L
1 + 2 · gl,l . Then from Eq. (A.4) we obtain

fL
γlL,l = lL .
f

(A.6)

l

We will use only L = 2:

L=2
fl=0 =

7
,
5

2
L=2
fl=2 = √ ,
7

6
L=2
fl=4 = √ .
35

In the main text the superscript L=2 on flL=2 is dropped for simplicity.

71

(A.7)

Appendix B
Conventions
Our convention for the Wigner-Eckart theorem is:

j m

λ
n1 j1 m1 |Tµ |n2 j2 m2 = Cj 1m 1,λµ · n1 j1 T λ n2 j2 ,
2 2

(B.1)

j m

where Cj 1m 1,λµ is the Clebsch-Gordan coefficients.
2 2
†

The reduced matrix element of qλµ (5.2) is
†


 l +λ−l
 i

(−)j +λ−j




nlj qλ n l j =
j1
(2λ+1)(2j +1)
·C 2
4π(2j+1)
j 1 ,λ0
2

dr r2 fλ (r)Rnlj (r)Rn l j (r), l + λ − l is even,

0,

l + λ − l is odd,

where the single-particle wave functions |nljm are defined as

ψnljm = Rnlj (r) ·
ml ms

jm
Clm ,sm il Ylm (θ, φ) χms ,
s
l
l

72

in which spin s = 1 , Rnlj (r) is a real function, and a factor il is included.
2
In this work we have used the matrix elements of the realistic quadrupole operator in
the harmonic oscillator single-particle basis. In this case λ = 2, fλ (r) = r2 , Rnlj (r) is
independent of j, and n = 2nr + l is the major-shell quantum number. The non-vanishing
matrix elements of nlj q † n l j

have n − n = −2, 0, 2, and l − l = −2, 0, 2. For these

combinations the symmetric radial integral becomes
















dr r4 Rnl (r)Rn l (r) = b2 ·















where b =

n + 3,
2

n = n, l = l,

− (n + l + 3)(n − l),

n = n, l = l + 2,

−1
2

(n + l + 3)(n − l + 2),

n = n + 2, l = l,

1
2

(n + l + 3)(n + l + 5),

n = n + 2, l = l + 2,

1
2

(n − l + 2)(n − l + 4),

n = n + 2, l = l − 2,

mΩ0 is the length parameter, Ω0 is the harmonic oscillator frequency. As

†
mentioned at the beginning of Sec. 5.2, the factor b2 will be combined with κ to make qµ

dimensionless.
The transition rates B(E2) is defined as
2Jf + 1
| f Q i |2
B(E2, Ji → Jf ) =
2Ji + 1
In Fig. 5.3 the unit of B(E2, 0+ → 2+ ) is given in e2 b2 , where 1b = 1barn = 100fm2 .

73

(B.2)

Appendix C
Systems Near Critical Point
Anharmonicities become important when the harmonic potential ω 2 α2 /2 becomes small or
negative. This is the case in realistic medium and heavy nuclei away from magic numbers
[28]. In the critical region ω 2 ≈ 0, we are able to determine the cubic potential term Λ(30)
and the quartic potential term Λ(40) as shown below.
We make an assumption in the spirit of Landau phase transition theory: in Eq. (2.4),
the leading potential term ω 2 α2 /2 vanishes at the critical point, while other higher order
terms Λ(mn) remain finite. Taylor expanding Λ(mn) over ω 2 ,

(mn)

Λ(mn) = Λc

(mn)

the leading constant term Λc

(mn) 2
(mn) 4
ω + Λ2
ω + ...,

+ Λ1

(C.1)

is finite. The stability of the system is restored by higher
(mn)

order anharmonicities, e.g. Λ(40) α4 /4; thus the generalized density matrices r12

are also

(mn)

finite. Again we call the finite leading constant term in a Taylor expansion rc12 . For
.
convenience, we use = instead of = if an equation is correct in constant terms but not in ω 2
terms or higher.

74

C.1

Harmonic Order

Keeping only the constant terms of Eq. (2.23) we have

(10) .
w34 =
12

n
(10)
V3214 12 w12 .
e12

(C.2)

Defining a square matrix
n
D(34),(12) ≡ δ(12),(34) − V3214 12 , (e3 = e4 , e1 = e2 ),
e
12

(C.3)

.
(10)
the e3 = e4 part of Eq. (C.2) is written as Dw(10) = 0. Because the quantities wc12
.
do not vanish, we have Det[D] = 0. Thus DT , the transpose matrix of D, has a zero
eigenvalue. More accurately, DT has an eigenvalue of order ω 2 ; because Det[D], the product
of all eigenvalues of DT , is of order ω 2 . Assume that the eigenvector corresponding to this
eigenvalue is η34 :
.
DT η = 0 η = 0

C.2

⇒

.
η T D = 0.

(C.4)

Cubic Anharmonicity

Keeping only the constant terms of Eq. (2.14) with (rs) = (20), we have

(01) .
(20)
(20)
−2iΛ(30) r12 = e12 r12 − n12 w12 + 2[w(10) , r(10) ]12 .

75

(C.5)

(20)

Calculating w34

from Eq. (C.5), the e3 = e4 part is written as
.
D · w(20) = C,

(C.6)

where D is defined in Eq. (C.3), C contains Λ(30) and lower-order quantities. Multiplying
Eq. (C.6) from the left by η T and using Eq. (C.4) we obtain
(01)

Λ(30) ·
e1 =e2 e3 =e4

i

η34 V3214
e1 =e2 e3 =e4

r
η34 V3214 12
e12

.
=

[w(10) , r(10) ]12 1
−
η34 V31 14 [p, r(10) ]11 ,
e12
2 e =e
1 1 e3 =e4
(10)

(30)

where p is given by n12 p12 = −ir12 . Eq. (C.7) gives Λc

(20)

+ O(ω 2 ). Then w12

(C.7)

is solved

from Eq. (C.6) with an overall factor (temporarily call δ (20) ) still undetermined.

(11)

Neglecting the −iω 2 r(02) term, w12

(11)

(then r12 ) are solved from Eq. (2.14) with (rs) =

(11)
(11), as a function of Λ(12) and δ (20) . We emphasize that in solving w12 the coefficient
(11)

matrix is actually not singular after combining w21

(11)

(11)

= (w12 )∗ with w12 , because w(11)

have symmetries different from w(10) in Eq. (C.2).

.
From Eq. (2.14) with (rs) = (02), we obtain an equation D · w(02) = ... Multiplying it
.
from the left by η T we fix δ (20) as a function of Λ(12) . Then from the equation D · w(02) = ...
(02)

we solve for w12

as a function of Λ(12) , with an overall factor δ (02) still undetermined.

In summary, there remain two undetermined parameters in this order: Λ(12) and an
(02)
overall factor δ (02) in w12 . The latter is Q(02) in the factorizable force model (Sec. 3.2).

76

C.3

Quartic Anharmonicity

Similarly, we obtain from Eq. (2.14) with (rs) = (30) and (rs) = (21):
(01)

Λ(40) ·

r
. i
(20)
η34 V31 14 r
η34 V3214 12 = Λ(12)
11
e12
3
e1 =e1 e3 =e4
e1 =e2 e3 =e4
(11)

−Λ(30)

r
i
(02)
η34 V3214 12 − Λ(30)
η34 V31 14 r
11
e12
3
e1 =e1 e3 =e4
e1 =e2 e3 =e4
+
−
−

i
2

η34 V3214
e1 =e2 e3 =e4

[w(20) , r(10) ]12 + [w(10) , r(20) ]12
e12

1
η34 V31 14 ([w(20) , r(01) ]11 + [w(01) , r(20) ]11 )
6 e =e
1 1 e3 =e4

1
η34 V31 14 ([w(11) , r(10) ]11 + [w(10) , r(11) ]11 ).
3 e =e
1 1 e3 =e4

(40)

Eq. (C.8) gives Λc

(C.8)

+ O(ω 2 ). There is one unknown parameter Λ(12) ; quantities (r/w)(20)

and (r/w)(11) depend implicitly on Λ(12) . Those terms with the cubic anharmonicity Λ(12)
would be small in the case of spherical nuclei, as shown in chapter 5.

In summary, this appendix fixes the cubic potential Λ(30) (C.7) and the quartic potential
Λ(40) (C.8) near the critical point ω 2 ≈ 0, by considering the leading terms of the GDM
equation of motion.

77

Appendix D
Mode Coupling
In many soft nuclei there exists a low-lying octupole (3− ) mode. It can interact strongly
with the quadrupole (2+ ) mode, and both of them should be kept in the collective subspace.
For convenience we still use αµ , πµ for the quadrupole mode; and use αµ , πµ for the octupole
ˆ ˆ
mode. The collective bosonic Hamiltonian replacing Eq. (2.4) is
1√
ω2 √
H=
5(α × α)0 +
5(π × π)0
0
0
2
2
ω2 √
ˆ
1√
7(ˆ × α)0 +
α ˆ 0
7(ˆ × π )0
π ˆ 0
+
2
2
Λ(10|20) √
+
7(α × (ˆ × α)2 )0 + . . .
α ˆ 0
2

(D.1)

Λ(10|20) is the most important mode-coupling term in the case of soft vibrations with large
amplitudes.
Following the procedure of chapter 2 and appendix C, we are able to determine the
leading constant term of Λ(10|20) in a Taylor expansion over both ω 2 and ω 2 [see Eq. (C.1)].
ˆ
Below we give the result in the quadrupole plus pairing model. The microscopic Hamiltonian

78

is:

H=
1
+
4

1234 µ

G
†
( 1 − µ)a1 a1 −
4

12
1
†
†
†
†
† †
(−κqµ14 qµ23 + κqµ13 qµ24 − κqµ14 qµ23 + κqµ13 qµ24 )a1 a2 a3 a4 .
ˆˆ ˆ
ˆˆ ˆ

Approximately, this Hamiltonian can be written as H ≈
1ˆ
2κ

† †
1

a1 a˜ a˜ a2
2

†
1
†
1 ( 1 −µ)a1 a1 −GP P − 2 κ

(D.2)

†
µ Qµ Qµ −

ˆ† ˆ
ˆ
µ Qµ Qµ , the difference is in a one-body term originating from the Q · Q part. κ is the

strength of the octupole force. The mean field is determined by the HFB equation. In the
harmonic order the two modes do not mix, the octupole mode satisfies the same QRPA
equation (5.18) and normalization condition (5.19) as the quadrupole mode, with necessary
changes. In the next order we have the main result:

.
Λ(10|20) = 2κQ(10) (ˆ Q(10) )2 ·
κˆ

·



 3

n1 j1 n2 j2 n3 j3

2



3 


 j j j 

1
2
3

·

(2j1 + 1)(2j2 + 1)(2j3 + 1)
7

(1)† ˆ(1)† †
†
ˆ(1)† ˆ(1)†
ˆ
2 ξ 13 ξ 32 η 21 + η 13 ξ 32 ξ 21 .

(D.3)

The octupole operator q connects single-particle levels with opposite parity, thus the intruder
ˆ
level becomes important. This may destroy in Eq. (D.3) symmetry with respect to the
Fermi surface. Three-body forces will contribute to the Λ(10|20) term quite differently. The
coexistence and interaction of soft quadrupole and octupole modes could be important to
the search of mechanisms for many-body enhancement of the nuclear Schiff moment and the
atomic electric dipole moment [29].

79

Appendix E
Odd-Mass Nuclei
In the odd-mass nuclei the coupling between the unpaired nucleon and the collective modes
is important. An effective Hamiltonian can be derived straightforwardly within the GDM
formalism. The “collective” subspace in this case is spanned by the direct product of the
unpaired nucleon and the phonons

†

†

b1 |Cn ≡ b1 ⊗ (A† )n |Φ

(1 ∈ Ω, 0 ≤ n ≤ nmax ),

(E.1)

†

where b1 is the quasi-particle creation operator (4.21), and |Φ is the quasi-particle vacuum.
Within this subspace the effective Hamiltonian is calculated as

†

Ci |He(12) |Cj ≡ Ci |b1 |H|b2 |Cj
(00)

(10)

(20) α

(01)

Ci |Λ(12) + Λ(12) α + Λ(12) π + Λ(12)

2

2

(02) π

+ Λ(12)

2

2

+ ...|Cj ,

(E.2)

where we substitute the nucleonic Hamiltonian (2.1), then use Eq. (2.8) and its generalization
to three-body density matrix. In the “collective” subspace (E.1) we diagonalize the resulting

80

effective Hamiltonian He (E.2).
The above treatment could also be used in calculating the states of “single-particle”
nature in even-even nuclei, as mentioned in CONCLUSION.

81

BIBLIOGRAPHY

82

BIBLIOGRAPHY
[1] A. Bohr and B. Mottelson, Nuclear Structure (Benjamin, New York, 1975), Vol. 2.
[2] L. Prochniak and S.G. Rohozinski, J. Phys. G: Nucl. Part. Phys. 36, 123101 (2009).
/Topical review./
[3] A. Arima and F. Iachello, Ann. Rev. Nucl. Part. Sci. 31, 75 (1981); F. Iachello and A.
Arima, the Interacting Boson Model (Cambridge Univ. Press, 1987).
[4] A. Kerman and A. Klein, Phys. Rev. 132, 1326 (1963).
[5] S.T. Belyaev and V.G. Zelevinsky, Yad. Fiz. 11, 741 (1970) [Sov. J. Nucl. Phys. 11, 416
(1970)]; Yad. Fiz. 16, 1195 (1972) [Sov. J. Nucl. Phys. 16, 657 (1973)]; Yad. Fiz. 17,
525 (1973) [Sov. J. Nucl. Phys. 17, 269 (1973)].
[6] V.G. Zelevinsky, Prog. Theor. Phys. Suppl. 74-75, 251 (1983).
[7] M.I. Shtokman, Yad. Fiz. 22, 479 (1975) [Sov. J. Nucl. Phys. 22, 247 (1976)].
[8] V.G. Zelevinsky, Nucl. Phys. A344, 109 (1980).
[9] V.G. Zelevinsky, Nucl. Phys. A337, 40 (1980).
[10] P.N. Isaev, Yad. Fiz. 32, 978 (1980) [Sov. J. Nucl. Phys. 32, 5056 (1980)]; Yad. Fiz. 34,
717 (1981) [Sov. J. Nucl. Phys. 34, 399 (1981)].
[11] A. Klein and A.S. Umar, Phys. Rev. C 35, 1672 (1987).
[12] L.Y. Jia, Phys. Rev. C 84, 024318 (2011).
[13] L. Y. Jia and V. G. Zelevinsky, Phys. Rev. C 84, 064311 (2011).

83

[14] L. Y. Jia and V. G. Zelevinsky, Phys. Rev. C 86, 014315 (2012).
[15] V. G. Zelevinsky, Int. J. Mod. Phys. E2, 273 (1993).
[16] H. J. Lipkin, N. Meshkov and A. J. Glick, Nucl. Phys. 62, 188 (1965).
[17] A. Klein and E.R. Marshalek, Rev. Mod. Phys. 63, 375 (1991).
[18] O. K. Vorov and V. G. Zelevinsky, Nucl. Phys. A439, 207 (1985).
[19] S. T. Belyaev, Nucl. Phys. 64, 17 (1965).
[20] A. Bohr, B. R. Mottelson and D. Pines, Phys. Rev. 110, 936 (1958).
[21] S. T. Belyaev, Mat. Fys. Medd. Kgl. Dan. Vidensk. Selsk. 31, No. 11 (1959).
[22] NuShellX@MSU, B. A. Brown and W. D. M. Rae,
http://www.nscl.msu.edu/ brown/resources/resources.html
[23] B. A. Brown, N. J. Stone, J. R. Stone, I. S. Towner, and M. Hjorth-Jensen, Phys. Rev.
C 71, 044317 (2005).
[24] T. Trivedi, P. C. Srivastava, D. Negi, and I. Mehrotra, Int. J. Mod. Phys. E21, 1250049
(2012).
[25] J. Blomqvist and A. Molinari, Nucl. Phys. A106, 545 (1968).
[26] National Nuclear Data Center, http://www.nndc.bnl.gov/.
[27] B. G. Carlsson, J. Toivanen, and A. Pastore, Phys. Rev. C 86, 014307 (2012).
[28] L. S. Kisslinger and R. A. Sorensen, Rev. Mod. Phys. 35, 853 (1963).
[29] N. Auerbach and V. Zelevinsky, J. Phys. G: Nucl. Part. Phys. 35, 093101 (2008).

84