Canonical-basis Time-Dependent
Hartree-Fock-Bogoliubov Theory
andLinear-Response Calculation
for Light to Heavy Nuclei

Shuichiro EBATA
February 2011

Canonical-basis Time-Dependent
Hartree-Fock-Bogoliubov Theory
andLinear-Response Calculation
for Light to Heavy Nuclei

Shuichiro EBATA( Doctoral Program in Physics )Submitted to the Graduate School ofPure and Applied Sciencesin Partial Fulfillment of the Requirements
for the Degree of Doctor of Philosophy in
Scienceat theUniversity of Tsukuba

Abstract
We have developed a new theory to study excited states and dynamics of atomic nuclei, which
we call the Canonical-basis time-dependent Hartree-Fock-Bogoliubov (Cb-TDHFB). The CbTDHFB is a simpliﬁed theory of the full TDHFB theory, treating pairing energy functional witha BCS-like approximation which is assumed to be diagonal in the canonical basis. Implementing
the theory with a real-space and real-time method, we can describe nuclear excitations anddynamics taking account of deformation eﬀects and pairing correlations, even for heavy nuclei
with a reasonable computational costs. In this thesis, we ﬁrst present a derivation of the CbTDHFB theory and we apply the method to linear-response calculations for even-even light toheavy nuclei and demonstrate its capability and accuracy by comparing our results with recent
calculations of the quasi-particle random-phase approximation (QRPA) which is equivalent tothe small amplitude limit of TDHFB, with Skyrme functionals.
We derived the Cb-TDHFB equations starting with the full TDHFB equations, introducing
a simpliﬁed form for the pairing energy functional. In the derivation, we start with the ordinary
TDHFB equation written in terms of the time-dependent generalized density matrix R(t) whichis composed of the time-dependent normal density matrix ρ(t) and the pairing tensor κ(t). We
then rewrite the equation in the canonical basis, which diagonalize ρ(t). The derived TDHFB
equations in the canonical basis has a very simple form. However, they are not useful practicallyfor a general case since it costs much to obtain the canonical basis themselves at each time.
Solving the time-dependent equation for the canonical basis is as diﬃcult as solving the full
TDHFB equation. This diﬃculty disappears if we employ a simpliﬁed pairing energy density
functional, which is equivalent to that employed in the BCS approximation for ground statecalculations.
For linear response calculations, we present results of the Cb-TDHFB calculations for excitedstates of some spherical and deformed nuclei to show performance of the method. We show
strength functions of isovector dipole and isoscalar quadrupole excitations for some spherical anddeformed nuclei. First, we show the comparison results of Cb-TDHFB with the QRPA in order
to conﬁrm the justiﬁcation of canonical-basis formulation and the reduction of computational
cost with this new framework. We show isoscalar quadrupole strength functions of neutron-richnucleus 34 Mg for conﬁrmation of our new method, and also discuss full self-consistent calculationsof 24 Mg as an example of light nucleus and 208 Pb and 154 Sm as an example of heavy nuclei, inorder to be clear, especially, the eﬀects of spin-orbit and Coulomb residual interaction. Next, inorder to show typical strength functions for spherical and deformed nuclei, we discuss isovector
dipole and isoscalar quadrupole mode for 34 Si as an example of spherical nucleus, for 34 Mg as an
example of prolate deformed nucleus and for 24 Ne as an example of oblate deformed nucleus. For
a conﬁrmation of computational cost of the calculation using Cb-TDHFB, we isovector dipole
strength functions of deformed heavy nuclei 172 Yb and 236−240 U. We show then the comparison
our results with experimental data for photo-nuclear reaction cross sections. And we also show
systematic calculations low-energy E1 strength which is often called Pygmy dipole resonanceand which is well-known as a important strength to understand nucleosynthesis. In this thesis,
we discuss the appearance of low-energy E1 strength functions of C, O, Ne, Mg, Si, S, Ar and
Ca isotopes which are relatively light nuclei.

Contents
1 Introduction
1.1 Basic properties and theoretical models of atomic nuclei . . . . . . . . . . . . . .
1.2 New phenomena in unstable nuclei . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Purpose of the present study and outline of this paper . . . . . . . . . . . . . . .

1
1
4
6

2 Formulation
11
2.1 Time-dependent mean ﬁeld theory . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.1.1 Time-dependent Hartree-Fock theory (TDHF) . . . . . . . . . . . . . . . 12
2.1.2 Time-dependent Hartree-Fock-Bogoliubov theory (TDHFB) . . . . . . . . 12
2.2 Canonical-basis formulation of time-dependent mean ﬁeld theory . . . . . . . . . 14
2.2.1 TDHF equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2.2 Canonical-basis TDHFB (Cb-TDHFB) equations . . . . . . . . . . . . . . 16
2.3 Properties of the Cb-TDHFB equations . . . . . . . . . . . . . . . . . . . . . . . 18
2.3.1 Gauge invariance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.3.2 Conservation laws . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.3.3 Stationary solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.3.4 Small amplitude limit and the Nambu-Goldstone modes . . . . . . . . . . 20
2.4 Cb-TDHFB equations with a simple pairing energy functional and gauge condition 22
2.4.1 Pairing energy functional . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.4.2 Properties of Cb-TDHFB equations with Eg . . . . . . . . . . . . . . . . 23
3 Application of Cb-TDHFB to linear-response calculations25
3.1 Determination of the pairing strength g . . . . . . . . . . . . . . . . . . . . . . . 25
3.2 Linear-response calculation in real-time . . . . . . . . . . . . . . . . . . . . . . . 26
3.3 Details of numerical calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.3.1 Three-dimensional coordinate space representation . . . . . . . . . . . . . 28
3.3.2 Calculation for ground state . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.3.3 Energy density functional . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.3.4 Calculation of time propagation . . . . . . . . . . . . . . . . . . . . . . . . 31
4 Results and Discussion
4.1 Comparison with QRPA calculations . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Full self-consistent calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
i

33
34
35

4.3

4.4
4.5
4.6

Typical strength functions for spherical, prolate and oblate nucleus . . . . . . . .
4.3.1 Example of spherical nucleus (34 Si) . . . . . . . . . . . . . . . . . . . . . .
4.3.2 Example of prolate nucleus (34 Mg) . . . . . . . . . . . . . . . . . . . . . .
4.3.3 Example of oblate nucleus (24 Ne) . . . . . . . . . . . . . . . . . . . . . . .
Computational cost (172 Yb, 236,238,240 U) . . . . . . . . . . . . . . . . . . . . . . .
Comparison with experiments (16 O, 24−26 Mg, 40 Ca, 90 Zr, 144−154 Sm, 208 Pb) . . .
Pygmy resonance (C,O,Ne,Mg,Si,S,Ar,Ca isotopes) . . . . . . . . . . . . . . . . .

39
39
41
43
47
48
52

5 Summary and Future work

59

AppendicesA
Data Table and E1 strength distributions . . . . . . . . . . . . . . . . . . . . . .
A.1
Ground state of Carbon (Z = 6) isotopes . . . . . . . . . . . . . . . . . .
A.2
Ground state of Oxygen (Z = 8) isotopes . . . . . . . . . . . . . . . . . .
A.3
Ground state of Neon (Z = 10) isotopes . . . . . . . . . . . . . . . . . . .
A.4
Ground state of Magnesium (Z = 12) isotopes . . . . . . . . . . . . . . .
A.5
Ground state of Silicon (Z = 14) isotopes . . . . . . . . . . . . . . . . . .
A.6
Ground state of Sulfur (Z = 16) isotopes . . . . . . . . . . . . . . . . . .
A.7
Ground state of Argon (Z = 18) isotopes . . . . . . . . . . . . . . . . . .
A.8
Ground state of Calcium (Z = 20) isotopes . . . . . . . . . . . . . . . . .
A.9
E1 strength distributions of Carbon (Z = 6) isotopes (N = 2 − 16) . . . .
A.10 E1 strength distributions of Oxygen (Z = 8) isotopes (N = 6 − 20) . . . .
A.11 E1 strength distributions of Neon (Z = 10) isotopes (N = 6 − 22) . . . .
A.12 E1 strength distributions of Magnesium (Z = 12) isotopes (N = 6 − 28) .
A.13 E1 strength distributions of Silicon (Z = 14) isotopes (N = 10 − 32) . . .
A.14 E1 strength distributions of Sulfur (Z = 16) isotopes (N = 10 − 32) . . .
A.15 E1 strength distributions of Argon (Z = 18) isotopes (N = 16 − 36) . . .
A.16 E1 strength distributions of Calcium (Z = 20) isotopes (N = 14 − 44) . .
BImaginary-Time method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
B.1
Gradient method and Imaginary-time method . . . . . . . . . . . . . . . .
B.2
Constraint imaginary-time method . . . . . . . . . . . . . . . . . . . . . .
CBCS equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
DIntrinsic frame and Deformation parameters . . . . . . . . . . . . . . . . . . . . .
D.1
Intrinsic frame . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
D.2
Quadrupole deformation parameter (β, γ) . . . . . . . . . . . . . . . . . .
D.3
β- and γ-vibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
D.4
Deﬁnition of quadrupole deformation parameter (β, γ) . . . . . . . . . . .

67
67
67
68
68
68
69
69
70
70
71
72
73
74
75
76
77
78
81
81
82
87
90
90
90
92
92

Bibliography

95

iiChapter 1

Introduction
The main purpose of the present paper is a construction of a new theoretical method which
provides a practical and universal description of nuclear reaction dynamics and elementary
modes of excitation. In order to achieve this, the theory should be able to describe nuclear
many-body dynamics with pairing correlations (superﬂuidity) which are very important forheavy nuclei. In this chapter, we give a brief introduction of basic properties of stable nuclei
and of theoretical approaches developed for description of nuclear dynamics. Then, we present
new features in unstable nuclei, to which the present method should be applicable.

1.1

Basic properties and theoretical models of atomic nuclei

Atomic nucleus is a self-binding ﬁnite quantum many-body system of size of several femtometers.
It is composed of two kinds of fermions (protons and neutrons) with spin 1/2 with (roughly)
the identical masses. The general term for protons and neutrons is “nucleons”. The number ofnucleons in the biggest nucleus in nature is a few hundreds at most. However, it is still impossible
to describe properties of heavy nuclei, starting from the bare nucleon-nucleon interaction. Thisis mainly because the Fermions require the anti-symmetrization for wave functions and thenucleon-nucleon interaction has a repulsive core in a short distance. Even if we use the latest
supercomputer, the current limit of the ab-initio calculation is nuclei with the mass number (A)less than 12 [5]. There are about 3,000 nuclei observed in experiments (Fig. 1.1), and yet more
to be discovered. Therefore, the development of eﬀective models is essential to understand avariety of nuclear phenomena.
Nucleus has “magic numbers”, which are analogous to those in atoms. This suggests thatthe mean free path of a nucleon is much larger than the mean distance between nucleons.
Thus, the Fermi gas model could be a useful starting point to discuss nuclear properties. The
magic numbers were explained with the j −j coupling shell model, proposed by M.Mayer and
J.H.D.Jensen [6], with a strong spin-orbit potential. Therefore, a crude approximation of thenucleus is given by independent nucleons freely moving in an attractive average potential.
The shell model, originally proposed by Mayer and Jensen, assumed a spherical potential.
This model qualitatively explains many characteristic properties of nuclei near the magic numbers in simple pictures, including single-particle (hole) states (stripping and pickup reactions),
1

Chapter 1. Introduction

Figure 1.1: Nuclear Chart of observed isotopes. Total number of nuclei is 3175 in this ﬁgure.
Each square means one nucleus and black ones correspond to stable nuclei. The numbers 2,
8, 20, . . . are called “magic number”. This ﬁgure is taken from National Nuclear Data Center
Home Page (http://www.nndc.bnl.gov/chart/)
magnetic moments, and low-energy spectra in odd-A nuclei. However, for nuclei whose proton
and/or neutron numbers are far from the magic numbers, a simple interpretation according to
the spherical shell model fails to explain nuclear properties. For instance, they were found tohave the electrical quadrupole moments much larger than those predicted by the shell model [7].
It cannot be explained by the independent nucleons in the spherical mean-ﬁeld potential. This
fact suggests that the nucleons collectively move to enhance the quadrupole moment, which may
eventually lead to a static deformation of the nucleus.
The deformation of the nucleus introduces an anisotropy that makes possible to specify
an orientation of the whole system. Then, it generates the rotational degrees of freedom andexhibits rotational spectra. Vast amount of experimental data can be understood by a modelthat the nucleons are moving independently in a deformed potential (Nilsson model) [8]. Figure
1.2 shows the adiabatic changes of single-particle levels in a deformed Nilsson potential withrespect to quadrupole deformation parameter β [9].
The nuclear deformation results from the quantum shell eﬀects, which is a consequence of
couplings between the single-particle motion and the deformed potential. This leads to an idea ofthe collective (uniﬁed) model by A.Bohr and B.R.Mottelson [10, 11]. Nucleons are independently
moving in an average nuclear potential, while the shape of the potential slowly changes (vibrates).
The shape change of the potential aﬀects the single-particle motion, then, the change in singleparticle motion leads to modiﬁcation of the potential. This feedback mechanism is caused by2

1.1. Basic properties and theoretical models of atomic nuclei

Figure 1.2: Neutron single-particle levels in the deformed Woods-Saxon potential as functionsof the quadrupole deformation parameter β. This ﬁgure is taken from Ref.[9].
the coupling between them. Therefore, the Hamiltonian of the collective model contains thefollowing three parts; the collective part to describe the vibrations of the potential shape, thenucleonic part for the single-particle motion in the potential, and the coupling terms between
them.
The time-dependent mean-ﬁeld theory can be regarded as a microscopic realization of thecollective model, using only nucleons’ degrees of freedom1 . The time-dependent mean-ﬁeld
potential is constructed from the microscopic interaction among nucleons, while the nucleons
are excited by motion of the mean-ﬁeld potential. The self-consistency between the mean ﬁeld
and one-body density represents the coupling terms. There have been extensive studies using
the time-dependent Hartree-Fock (TDHF) method, especially in 1970’s and 1980’s [12]. At later
chapters, we recapitulate the TDHF theory. Many of those TDHF calculations were performedin the coordinate-space representation [13, 14, 15, 16, 17] to describe nuclear fusion and heavyion collision dynamics, recently they have been also carried out [18, 19, 20, 21, 22]. In the small
amplitude limit, the TDHF is able to describe nuclear excitation spectra and transition densities,
equivalent to the random-phase approximation (RPA) method [3]. Real-time calculations of thesmall-amplitude TDHF have been carried out recently [23, 24].
The RPA and TDHF are good tools to describe excited state and collective motion in nuclei,
1

Since the bare nucleon-nucleon interaction has a strong repulsive core, a direct application of the meanﬁeld approximation breaks down. Thus, historically, eﬀective interactions, that renormalize high-momentum
components of the two-body interaction, are used for the mean-ﬁeld calculations [12]3Chapter 1. Introduction
however, they can not deal with an important correlation in nuclei, the pairing correlation.
For instance, the nuclear binding energies are known to have a systematic variation dependingon the evenness or oddness of the proton (Z) and neutron (N ) numbers (Fig. 1.3). This is
known as “odd-even eﬀect” even present in the old empirical mass formula of by H.Bethe andR.Weizsäcker [25]. This odd-even eﬀect may be described in terms of the pairing correlation. For
the description of the nuclear pairing, A.Bohr, B.R.Mottelson and D.Pines successfully applied
the BCS theory to nuclei [26] in the following year after the BCS theory was published [27]. The
gap energy can be seen in low-lying excitation spectra (Fig. 1.4). The pairing condensate leads
to the nuclear superﬂuidity and the moment of inertia is expected to become smaller than therigid-body value. This is also conﬁrmed in experiment (Fig. 1.5).
The pairing is know to be very important for nuclear collective dynamics as well. Figure
1.6 presents lifetimes of spontaneous ﬁssioning nuclei. The lifetimes of even-even nuclei aresystematically larger than those of odd-even nuclei, and those of odd-odd nuclei are even shorter[28, 29] (Fig. 1.6). The nuclear softness (deformability) signiﬁcantly depends on the evenness oroddness of Z and N . The pairing correlation is supposed to be responsible for this. From theseobservations, we believe that the inclusion of the pairing is essential for description of dynamical
properties of heavy nuclei.
N-Z=21
N-Z=23

Sn(N,Z) [MeV]

Sn(N,Z)=B(N,Z) - B(N-1,Z)
10

5

70

75

80
Neutron Number

85

90

Figure 1.4: Low-lying excited states in the even-even
Figure 1.3: Neutron separation energies Sn (N, Z) [1].

tin isotopes below 2.8 MeV. This ﬁgure is taken from
Ref.[31].

The data are taken from Ref.[30]1.2

New phenomena in unstable nuclei

Structure and reaction of unstable nuclei are of great interest as current topics in nuclear physics.
Recent developments of radioactive isotope facilities in the world, including RI Beam Factory(RIBF) at RIKEN, enable generation and detection of diverse rare isotopes. We are now able toexperimentally investigate properties of these unknown nuclei. It demands a theoretical modelthat is capable of a non-empirical description of a variety of nuclear properties. This is one of
the motivations for the new approach we are going to present in this article. In this section, we
recapitulate interesting properties of these unstable nuclei, especially focusing on neutron skinand halo structure [33, 34], and new magic numbers [35, 36].
The density distribution of neutrons is almost identical to that of protons in stable nuclei.
41.2. New phenomena in unstable nuclei

Figure 1.5: Moments of inertia of deformed nuclei
as function of mass A in the so-called rare-earth region(150 < A < 190). Solid lines indicate experimental data,crossed line shows the rigid momentum of inertia anddashed and dot-dashed lines show the calculation withpairing correlation. This ﬁgure is taken from Ref.[32]Figure 1.6: Spontaneous ﬁssion half-lives as a functionof Z 2 /A. From left to right, the lines indicates even-even,
odd-even and odd-odd nuclei. This ﬁgure is taken from
Ref.[28].

There is little diﬀerence in Fermi levels of neutrons and protons. However, the Fermi level forneutrons in neutron-rich nuclei becomes much higher than that of protons. This may lead to asigniﬁcant diﬀerence in the density distribution of protons and neutrons, to produce phenomena
called neutron skin and neutron halo. The neutron skin structure basically means that theneutrons and protons have diﬀerent radii, thus to create a nuclear surface region in which only
the neutrons exist. The neutron halo is a consequence of the weak binding of neutrons, having avery extended density distribution of neutrons because the wave function damps very weakly ina classically forbidden region, ψ(r) ∼ e−κr . The most famous example of the neutron halo is 11 Li[33](Fig. 1.7). The neutrons of the halo nucleus is very dilute in a surface area. Theoretically,
properties of the pairing in the dilute neutron matter is predicted to be very diﬀerent from thosein nuclear matter at the normal density. The pairing correlation at the normal density producesa Cooper pair of nucleons, that is spatially very extended, similar to the electron Cooper pairin superconductors. However, this property is drastically modiﬁed by changing the neutron
density. The Cooper pair at low density is becoming very compact, which is often called “dineutron correlation” [37]. Since two neutrons do not create a bound state in the vacuum, thereal Bose-Einstein Condensation (BEC) is not realized in the neutron matter. However, theneutron density in the surface region of halo nuclei approximately corresponds to the BCS-BEC
cross over region [38].
Since the structure and the correlation in the ground state are signiﬁcantly modiﬁed in
the neutron-rich nuclei, we may expect a new type of elementary excitations. The pygmy
dipole resonance (PDR) is such an example expected to appear in unstable nuclei. The PDR
is a dipole oscillation whose frequency much lower than the giant dipole resonance (GDR).
It has been predicted as a collective soft-dipole excitation which is the oscillation of the coreagainst the weakly-bound neutrons [39, 40] (see Fig. 1.8), Despite of recent experimental andtheoretical studies, it has not yet been conﬁrmed the soft-dipole picture in neutron-rich nuclei.
5Chapter 1. Introduction

Figure 1.7: Observed matter root mean square (rms)radii for light nuclei. Lines connect isotopes only for
guides. 11 Li has a radius much lager than the other nuclei. This ﬁgure is taken from Ref.[33]Figure 1.8: Schematic picture of the soft-E1 resonance.
This ﬁgure is taken from Ref.[34].

The collectivity of the PDR is one of current issues in physics of unstable nuclei [41]. It shouldbe also emphasized that the PDR plays an important role for element synthesis [42, 43].
The nuclear deformation is mainly determined by the shell eﬀect. Therefore, we expectthat nuclei near the magic numbers are spherical, which has been conﬁrmed for stable nuclei.
However, experimental data of unstable nuclei reveal a nuclear mass region, often called “islandof inversion” [44], in which the nuclei are deformed even though they correspond to the magic
numbers. An example with the neutron magic number 20 is shown in Fig. 1.9. It is suggested
that, in this region, the magic number may change from N = 20 to N = 16 [35] (Fig. 1.10).
Similar phenomena are theoretically predicted in heavier neutron-rich regions as well (N = 40
with 17 ≤ Z ≤ 25) [45, 46, 47]. The shell structure and magic numbers are essential ingredientsfor understanding of nuclear structure. Since the shell eﬀect aﬀects the total binding energy,
they are also important for element synthesis reactions, such as the determination of the pathof rapid neutron capture process (r-process). A non-empirical model, which is able to describe
and determine an arbitrary nuclear shape, is desirable, because we may encounter unknown
phenomena. The discovery of the island of inversion gives us a lesson that it may be dangerous
to make an assumption according to the common sense in physics of stable nuclei.

1.3

Purpose of the present study and outline of this paper

The purpose of this study is to construct a feasible method which can universally describe
properties of nuclei in the entire nuclear chart, including stable and unstable nuclei. The methodshould be able to describe properties of stable nuclei and also applicable to new phenomena inunstable nuclei. In order to predict properties of unknown nuclei, it is essential for the model
to possess a predictive power without adjusting parameters from nucleus to nucleus. With6

1.3. Purpose of the present study and outline of this paper

Figure 1.9: Partial periodic table highlighting the “island of inversion” centered at 32 Na. The “inversion”
means that the ordering of some single-particle levels areinverted in this region. The magic numbers Z = 8 and
N = 20 are emphasized with bold lines. This ﬁgure istaken from Ref.[44]Figure 1.10: Neutron number (N ) dependence of observed neutron separation energy Sn for (a) nuclei with
odd N and even Z, and for (b) nuclei with odd N and
odd Z. The arrows indicate to magic numbers. Thisﬁgure is taken from Ref.[35].

this respect, the self-consistent mean-ﬁeld (density-functional) models are superior to otherapproaches. They are also computationally feasible. Using the parallel supercomputers, asystematic calculation of the ground-state properties of all nuclei can be achieved in several
hours [48].
In order to describe reaction dynamics and nuclear response properties, the time-dependentmean-ﬁeld approach would be a good candidate. The TDHF theory, previously mentioned, canbe extended to the one applicable to superﬂuid nuclei. The inclusion of the pairing correlation
leads to the “time-dependent Hartree-Fock-Bogoliubov” (TDHFB) theory. However, it requires
enormous computational resources because the number of orbitals to treat is, in principle, inﬁnite. Only recently, a few attempts of the TDHFB have been done for realistic applications,but either with a spherical restriction [49] or with a small model space [50].
The small-amplitude limit of the TDHFB is identical to the quasi-particle RPA (QRPA). For
spherical nuclei, there have been extensive studies with the QRPA to investigate excited statesin stable and unstable nuclei [51, 39]. The QRPA methods with the coupling with continuum
state are also in progress [52, 53, 54]. Very recently, the QRPA calculations have become possible
for nuclei with axially deformed ground states [55, 56, 57]. However, the fully self-consistent
QRPA calculation including the spin-orbit and residual Coulomb interactions requires a massive
computational cost. For example, in Ref.[57], the calculation was performed in a massively
parallel computer using ten thousands of CPU’s.
In this thesis, we propose an alternative yet feasible approach, “Canonical-basis TDHFB”(Cb-TDHFB) in the three-dimensional coordinate-space representation. The Cb-TDHFB isformulated in the canonical-basis representation and an approximation for the pair potential
leads to a simple theory that requires the computational cost much smaller than the full TDHFB
calculation. We show a basic idea and the derivation of the Cb-TDHFB equations. Then, as theﬁrst numerical application, we apply the method to calculation of the linear response, to show
7

Chapter 1. Introduction
the strength distributions calculated with the real-time approach. We compare the result of the
Cb-TDHFB with those of recent QRPA calculations, to demonstrate the validity and accuracyof the Cb-TDHFB method. The dynamical eﬀects of the pairing correlation will be addressed.
This thesis is organized as follows. In Chapter 2, we derive the basic equations of Cb-TDHFB
method. According to a pairing energy functional, it may require a careful treatment on thetime-dependent phases of the canonical states. In Chapter 3, we show details of numericalcalculation. In Chapter 4, we show results for the strength functions of isovector dipole andisoscalar quadrupole modes. The subjects of nuclei are in light ones, carbon (C), oxygen (O),neon (Ne), magnesium (Mg), silicon (Si), sulfur (S), argon (Ar), calcium (Ca) isotopes, and inheavy ones are samarium (Sm) and uranium (U) isotopes and zirconium of 90 (90 Zr), ytterbiumof 172 (172 Yb) and lead of 208 (208 Pb). First, we compare our results with other calculations
in order to conﬁrm the accuracy of Cb-TDHFB method as time-dependent scheme. And weshow the basic properties of deformed nucleus and the eﬀects of pairing correlation on isovector
dipole and isoscalar quadrupole modes of each deformed nucleus. And we discuss the full selfconsistency and computational cost of Cb-TDHFB with calculation results, we compare CbTDHFB results with experimental data from light to heavy nuclei. Then we show the properties
and eﬀects of pairing on pygmy dipole resonances for nuclei up to A = 64. Finally, the summary
and future work of this study are given in Chapter 5.

8Chapter 2

Formulation
2.1

Time-dependent mean field theory

In this section, we show the derivations of time-dependent Hartree-Fock (TDHF) and timedependent Hartree-Fock-Bogoliubov (TDHFB) equations in order to lead to the canonical-basis
TDHFB. We need some preparations.
We start to derive TDHF(B) equations from time-dependent many-body Schrödinger equation. We set a time-dependent many-body wave function |Φ(t)⟩. When this |Φ(t)⟩ obey the
many-body Schrödinger equation, the time evolution of |Φ(t)⟩ is given by|Φ(t)⟩ = e−iHt/~ |Φ(0)⟩.

(2.1)

We assume Hamiltonian H is the sum of a kinetic energy and a two-body interaction. H canbe written with a Fermion creation and annihilation operator by
H=

∑
α,β

tαβ c†α cβ +1 ∑V̄αβµν c†α c†β cν cµ ,
4

(2.2)

α,β,µ,ν

where V̄αβµν = Vαβµν − Vαβνµ and Fermion operators obey {ci , c†j } = δij . We set the timedependent normal density matrix ρ(t) and pairing tensor κ(t) as
ρij (t) ≡ ⟨Φ(t)|c†j ci |Φ(t)⟩,(2.3)

κij (t) ≡ ⟨Φ(t)|cj ci |Φ(t)⟩.

(2.4)

When we regard that the time-evolution of |Φ(t)⟩ obey Eq.(2.1), to obtain an equation of motion,
we can calculate their time derivative:i~∂ρij (t) = ⟨Φ(t)|[c†j ci , H]|Φ(t)⟩,
∂t

(2.5)

i~

∂
κij (t) = ⟨Φ(t)|[cj ci , H]|Φ(t)⟩.
∂t

(2.6)

11

Chapter 2. Formulation
In order to expand the commutation relation [c†j ci , H], we use followings[c†j ci , c†α cβ ] = c†j cβ δiα − c†α ci δβj ,[c†j ci , c†α c†β cν cµ ] = c†j c†β cν cµ δiα − c†j c†α cν cµ δiβ − c†α c†β ci cµ δjν + c†α c†β ci cν δjµ .
And in order to expand the commutation relation [cj ci , H], we use followings[cj ci , c†α cβ ] = cj cβ δiα − ci cβ δjα ,[cj ci , c†α c†β cν cµ ] = c†α cj cν cµ δiβ − c†α ci cν cµ δjβ + cj c†β cν cµ δiα − ci c†β cν cµ δjα .

2.1.1

Time-dependent Hartree-Fock theory (TDHF)

In the TDHF theory, the many-body wave function |Φ(t)⟩ is replaced with a single Slater determinant asN∏|Φ(t)⟩ ≡c†i |0⟩,(2.7)
i=1

where N is the number of particles in the system. This wave function is called Hartree-Fock(HF) state, and this replacement is Hartree-Fock approximation. HF state can not deal with κ,
so TDHF equation is written by only the time-derivative of ρ. Eq.(2.5) can be expanded thenin HF state as∑∑∂tiβ ρβj −tαj ρiαi~ ρij =
∂t
α
β
)
(
1 ∑
+
V̄iβµν ρνβ ρµj − ρµβ ρνj
2
β,µ,ν
)
(
1 ∑
−
V̄αβµj ρiβ ρµα − ρiα ρµβ
2
α,β,µ
∑
=
{(tiµ + Γiµ )ρµj − ρiµ (tµj + Γµj )} = [h, ρ]ij ,
(2.8)
µ

∑

where we introduce Γαµ ≡ βν V̄αβµν ρνβ called HF potential, h ≡ t + Γ so-called single-particle
Hamiltonian and the Eq.(2.8) is called TDHF equation.

2.1.2

Time-dependent Hartree-Fock-Bogoliubov theory (TDHFB)

We can derive the TDHFB equations with similar procedure for TDHF. First, we need to chosea quasi-particle vacuum |Φ(t)⟩ as many-body wave function. The quasi-particle operators βk† , βk
are usually deﬁned as
∑
∑
∗
∗ †
βk† ≡
Uαk c†α + Vαk ck , βk =Uαkcα + Vαk
ck ,
(2.9)
α

|Φ(t)⟩ ≡∏αβk† |0⟩.

(2.10)

k12

2.1. Time-dependent mean ﬁeld theory
The vacuum of this quasi-particle state |Φ(t)⟩ is called Hartree-Fock-Bogoliubov (HFB) state.
The replacement of many-body wave function to HFB state is HFB approximation. HFB state
can deal with pairing tensor κ. TDHFB equations are derived from time-derivatives of ρ and κ.
Eq.(2.5),(2.6) can be expanded in HFB state as
i~

∂
ρij
∂t

=

∑

tiβ ρβj −tαj ρiα

α

β

+

∑

(
)
1 ∑
V̄iβµν ρνβ ρµj − ρµβ ρνj + κ∗jβ κµν
2
β,µ,ν

−

(
)
1 ∑
V̄αβµj ρiβ ρµα − ρiα ρµβ + κ∗αβ κµi ,2(2.11)

α,β,µ

i~

∂
κij
∂t

=

∑

tiβ κβj − tjβ κβi +

β

+
−

1∑
V̄ijµν κµν
2 µ,ν

)
(
1 ∑
V̄αiµν ρjα κµν − ρνα κµj + ρµα κνj
2 α,µ,ν
)(1 ∑V̄αjµν ρiα κµν − ρνα κµj + ρµα κνi .
2 α,µ,ν

(2.12)∑∑
When we introduce Γαµ = βν V̄αβµν ρνβ and ∆αβ ≡ 12 µν V̄αβµν κµν called pair potential, so
Eq.(2.11) and (2.12) can be rewritten as
i~

∂
ρij
∂t

∂
i~ κij∂t=∑{}(tiµ + Γiµ )ρµj − ρiµ (tµj + Γµj ) + ∆iµ κ∗jµ − κµi ∆∗µj ,

µ

=

∑{

}
(tiµ + Γiµ )κµj − (tjµ + Γjµ )κµi − ∆iµ ρ∗µj − ρiµ ∆µj + ∆ij .

µWe introduce h = t + Γ to them, then we can get usual TDHFB equations as∂
ρ = [h, ρ] + κ∆∗ − ∆κ∗ ,∂t∂i~ κ = hκ + κh∗ + ∆(1 − ρ∗ ) − ρ∆.
∂ti~(2.13)(2.14)

Properties of TDHFBWe show that TDHFB conserves expectation values of total number and total energy. Total
number N and energy E can be written with using HFB many-body wave function |Φ⟩ as∑∑N ≡ ⟨Φ|N |Φ⟩ =⟨Φ|c†i ci |Φ⟩ =ρii = Tr{ρ},(2.15)iE ≡ ⟨Φ|H|Φ⟩ =

∑
α,β

i

{
}
1 ∗
1
1 ∗
1
hαβ ρβα − Γαβ ρβα − ∆αβ κβα = Tr hρ − Γρ − ∆ κ ,
2
2
2
2
13

(2.16)

Chapter 2. Formulation
where the notation Tr stands for taking the trace of the matrices and we assume the Hamiltonian
H given as (2.2). We expand above values with the small time-variations of ρ and κ. Thesevariations can be expressed TDHFB equations (2.13) and (2.14) asρ′ = ρ − iη([h, ρ] + κ∆∗ − ∆κ∗ ),κ′ = κ − iη(hκ + κh∗ + ∆(1 − ρ∗ ) − ρ∆),where η is equivalent to value of dt/~. Then, we can express time-variations δN as
δN = Tr{ρ′ − ρ} = −iηTr{[h, ρ] + κ∆∗ − ∆κ∗ } = 0.

(2.17)

It is given by following relations,[h, ρ]αβ =∑hαµ ρµβ − ραµ hµβ → 0 (with δαβ ),µTr{κ∆∗ − ∆κ∗ } =∑δαβ (κ∆∗ − ∆κ∗ )αβ

α,β

=

1 ∑
καµ V̄µαmn κ∗mn − V̄αµmn κmn κ∗µα = 0.
2 α,µ,m,n

And also we can express time-variation δE as
δE = − iηTr{h[h, ρ] + hκ∆∗ − h∆κ∗ }iη
−Tr {−∆∗ hκ − ∆∗ κh∗ − ∆∗ ∆ + ∆∗ ∆ρ∗ + ∆∗ ρ∆}2iη. −Tr { h∗ κ∗ ∆ + κ∗ h∆ + ∆∗ ∆ − ∆∗ ρ∆ − ρ∗ ∆∗ ∆} = 0,2where we use the relations h∗ = hT , κT = −κ and ∆T = −∆ and the time-variation from Γρ inEq.(2.16) is zero because Tr{[h, ρ]} = 0.

2.2

Canonical-basis formulation of time-dependent mean field
theoryIn this section, we derive the basic equations of Cb-TDHFB method. Using the time-dependent
variational principle, the similar equations were derived by Blocki and Flocard [58]. However,it was not clear that what kind of approximation was introduced and how they are diﬀerent
from the full TDHFB. We present here a suﬃcient condition to reduce the TDHFB equations
to those in a simple canonical form.
We start from the density-matrix equation of the TDHFB and ﬁnd equations for the canonicalbasis states and their occupation- and pair-probability factors. In order to clarify our heuristic
strategy, let us start from a simpler case without the pairing correlation.
14

2.2. Canonical-basis formulation of time-dependent mean ﬁeld theory

2.2.1

TDHF equation

The TDHF equation in the density-matrix formalism is written as Eq.(2.8) or [3]i∂ρ(t) = [h(t), ρ(t)] ,∂twhere ρ(t) and h(t) are the one-body density operator and the single-particle (Hartree-Fock)
Hamiltonian, respectively. We now express the one-body density using the time-dependent
canonical single-particle basis, {|ϕk (t)⟩}, which are assumed to be orthonormal (⟨ϕk (t)|ϕl (t)⟩ =δkl ).
ρ(t) =N∑|ϕk (t)⟩⟨ϕk (t)|,(2.18)

k=1

where N is the total particle number. Substituting this into Eq.(2.8), we have
N {
N
} ∑
∑
{h(t)|ϕk (t)⟩⟨ϕk (t)| − |ϕk (t)⟩⟨ϕk (t)|h(t)} .
i|ϕ̇k (t)⟩⟨ϕk (t)| + i|ϕk (t)⟩⟨ϕ̇k (t)| =k=1(2.19)

k=1

the inner product with |ϕk (t)⟩ leads to
(
P̂

)
∂
i − h(t) |ϕk (t)⟩ = 0∂t

k = 1, · · · , N,

(2.20)∑
with P̂ = 1 − N
k=1 |ϕk (t)⟩⟨ϕk (t)|. Here, we used the conservation of the orthonormal property
for the canonical states, d/dt⟨ϕk (t)|ϕl (t)⟩ = 0. This leads to the most general canonical-basis
TDHF (Cb-TDHF) equations∑∂|ϕl (t)⟩ηlk (t),i |ϕk (t)⟩ = h(t)|ϕk (t)⟩ −∂tN

k = 1, · · · , N,

(2.21)

l=1

where the matrix ηlk (t) is arbitrary but should be hermitian to conserve the orthonormal property. It is easy to see that the time evolution of the density does not depend on the choice ofηlk . This is related to the gauge invariance with respect to the unitary transformations among|ϕk (t)⟩ (k = 1, · · · , N ). The most common choice is ηlk = 0, which leads to the TDHF equation
shown in most textbooks.
15

Chapter 2. Formulation

2.2.2

Canonical-basis TDHFB (Cb-TDHFB) equations

We now derive Cb-TDHFB equations starting from TDHFB equations Eq.(2.13) and (2.14) 1 .
∂ρ(t) = [h(t), ρ(t)] + κ(t)∆∗ (t) − ∆(t)κ∗ (t),∂t∂i κ(t) = h(t)κ(t) + κ(t)h∗ (t) + ∆(t)(1 − ρ∗ (t)) − ρ(t)∆(t).
∂ti

At each instant of time, we may diagonalize the density operator ρ̂ in the orthonormal canonicalbasis, {ϕk (t), ϕk̄ (t)} with the occupation probabilities ρk . For the canonical states, we use thealphabetic indexes such as k for half of the total space indicated by k > 0. For each state withk > 0, there exists a “paired” state k̄ < 0 which is orthogonal to all the states with k > 0.
The set of states {ϕk , ϕk̄ } generate the whole single-particle space2 . We use the Greek lettersµ, ν, · · · for indexes of an adopted representation (complete set) for the single-particle states.
∑
The creation operator of particles at the state |ϕk (t)⟩ is expressed as ĉ†k (t) = µ ⟨µ|ϕk (t)⟩ĉ†µ ,and the TDHFB state is expressed in the canonical (BCS) form as}∏{|Ψ(t)⟩ =uk (t) + vk (t)c†k (t)c†k̄ (t) |0⟩.
(2.22)
k>0

For later purposes, it is convenient to introduce the following notations for two-particle states:⟨µν|ϕk (t)ϕk̄ (t)⟩ ≡ ⟨µ|ϕk (t)⟩⟨ν|ϕk̄ (t)⟩,(2.23)⟨⟨µν|ϕk (t)ϕk̄ (t)⟩⟩ ≡ ⟨µν|ϕk (t)ϕk̄ (t)⟩ − ⟨µν|ϕk̄ (t)ϕk (t)⟩,(2.24)

and for projection operator on a canonical pair of states (k, k̄),π̂k (t) ≡ |ϕk (t)⟩⟨ϕk (t)| + |ϕk̄ (t)⟩⟨ϕk̄ (t)|.
Then, it is easy to show the following properties (k, l > 0):∑⟨µν|ϕk (t)ϕk̄ (t)⟩⟨ϕl (t)ϕl̄ (t)|µν⟩ = δkl ,∑(2.25)(2.26)µν⟨⟨µν|ϕk (t)ϕk̄ (t)⟩⟩⟨⟨ϕl (t)ϕl̄ (t)|µν⟩⟩ = 2δkl ,(2.27)µν∑⟨⟨µσ|ϕk (t)ϕk̄ (t)⟩⟩⟨⟨ϕl (t)ϕl̄ (t)|νσ⟩⟩ = δkl ⟨µ|π̂(t)k |ν⟩,σ∑⟨⟨µσ|ϕk (t)ϕk̄ (t)⟩⟩⟨ν|π̂l |σ⟩ = δkl ⟨⟨µν|ϕk (t)ϕk̄ (t)⟩⟩.

(2.28)(2.29)σ1

In some case, the TDHFB equation can be written in terms of the generalized density matrix [4] as
i
(

where
R≡

ρ
−κ∗

∂
R = [H, R] ,
∂t

)
κ
∗ ,
1−ρ

(
H≡

2

h
−∆∗

)
∆
∗ .
−h

In the case without pairing (∆ = 0), the canonical pair becomes arbitrary as far as they have the same
occupation probability ρk that is either 1 or 0.

16

2.2. Canonical-basis formulation of time-dependent mean ﬁeld theory
Using these notations, the density and the pairing-tensor matrixes are given by∑ρµν (t) =ρk (t)⟨µ|π̂k (t)|ν⟩,(2.30)k>0κµν (t) =∑κk (t)⟨⟨µν|ϕk (t)ϕk̄ (t)⟩⟩,(2.31)

k>0

where ρk (t) and κk (t) are occupation and pair probabilities, respectively. In terms of the BCS
factors of (u, v) [3], they are given as ρk (t) = |vk (t)|2 and κk (t) = u∗k (t)vk (t). It should be notedthat the canonical pair of states, |ϕk (t)⟩ and |ϕk̄ (t)⟩, are assumed to be orthonormal but notnecessarily related with each other by the time reversal, |ϕk̄ ⟩ ̸= T |ϕk ⟩.
Thanks to the orthonormal property, we can invert Eqs.(2.30) and (2.31) for ρk and κk ,
ρk (t) =∑⟨ϕk (t)|µ⟩ρµν (t)⟨ν|ϕk (t)⟩ =⟨ϕk̄ (t)|µ⟩ρµν (t)⟨ν|ϕk̄ (t)⟩,(2.32)

1∑⟨⟨ϕk (t)ϕk̄ (t)|µν⟩⟩κµν (t).
⟨ϕk (t)ϕk̄ (t)|µν⟩κµν (t) =2 µνµν(2.33)µνκk (t) =∑µν∑With help of Eq.(2.30), the derivative of ρk (t) with respect to time t leads to
i

∑
dρµν
d
d
ρk (t) =⟨ϕk (t)|µ⟩i⟨ν|ϕk (t)⟩ + iρk (t) ⟨ϕk (t)|ϕk (t)⟩
dt
dt
dt
µν
∑
dρµν
=
⟨ϕk (t)|µ⟩i⟨ν|ϕk (t)⟩
dt
µν
∑{
}
=
κk (t)∆∗µν (t)⟨νµ|ϕk (t)ϕk̄ (t)⟩ + κ∗k (t)∆µν (t)⟨ϕk (t)ϕk̄ (t)|µν⟩ .

(2.34)µν

We used the assumption of norm conservation for the second equation, and used the TDHFB
equation (2.13) in the last equation. Since the pair potential ∆µν (t) is anti-symmetric, it is
written in a simple form asdρk (t) = κk (t)∆∗k (t) − κ∗k (t)∆k (t),dt∑
1∑∆k (t) ≡ −∆µν (t)⟨ϕk (t)ϕk̄ (t)|µν⟩ = −∆µν (t)⟨⟨ϕk (t)ϕk̄ (t)|µν⟩⟩.
2 µν
µν
i

In case that the pair potential is computed from a two-body interaction v as ∆µν (t) =the gap parameters, ∆k (t), are identical to those of the BCS approximation [3].
∆k (t) = −∑κl (t)(vkk̄,ll̄ − vkk̄,l̄l ) ≡ −l>0∑κl (t)v̄kk̄,ll̄ .

(2.35)(2.36)∑

αβ vµν,αβ καβ (t),(2.37)l>0

It should be noted here that the two-body matrix elements vkk̄,ll̄ (and the anti-symmetric v̄kk̄,ll̄ )are time-dependent because the canonical basis, (k, k̄) and (l, ¯l), are time-dependent.
17Chapter 2. Formulation
In the same way, we evaluate the time derivative of κk (t) as
(
)
∑
dκµν
d
dϕk
dϕk̄
i κk (t) =⟨ϕk (t)ϕk̄ (t)|µν⟩i+ iκk (t) ⟨|ϕk (t)⟩ + ⟨|ϕ (t)⟩ .
dt
dt
dt
dt k̄
µν

(2.38)

Then, using the TDHFB equation (2.14), we obtainidκk (t) = (ηk (t) + ηk̄ (t)) κk (t) + ∆k (t) (2ρk (t) − 1) ,dt

(2.39)kwhere ηk (t) ≡ ⟨ϕk (t)|h(t)|ϕk (t)⟩ + i⟨ ∂ϕ∂t |ϕk (t)⟩.
The time-dependent equations for ρk (t) and κk (t) are now given in rather simple forms as
Eqs.(2.35) and (2.39). So far, their derivation is solely based on the TDHFB equations, utilizing
the fact that ρµν (t) and κµν (t) can be expressed by the orthonormal canonical basis, |ϕk (t)⟩and |ϕk̄ (t)⟩, and their occupation and pair probabilities, ρk (t) and κk (t). However, in general,
the time evolution of the canonical basis is not given by a simple equation. Therefore, we now
introduce an assumption (approximation) that the pair potential is written as∑∆µν (t) = −∆k (t)⟨⟨µν|ϕk (t)ϕk̄ (t)⟩⟩.
(2.40)

k>0This satisﬁes Eq.(2.36), but in general, Eq.(2.36) can not be inverted because the two-particle
states |ϕk ϕk̄ ⟩ do not span the whole space. In other words, we only take into account the pairpotential of the “diagonal” parts in the canonical basis, ∆kl̄ = −∆k δkl . In the stationary limit(|ϕk̄ ⟩ = T |ϕk ⟩), this is equivalent to the ordinary BCS approximation (see Sec. 2.3.3).
With the approximation of Eq.(2.40), it is easy to see that the TDHFB equations, (2.13)
and (2.14), are consistent with the following equations:i∂|ϕk (t)⟩ = (h(t) − ηk (t))|ϕk (t)⟩,∂ti∂|ϕ (t)⟩ = (h(t) − ηk̄ (t))|ϕk̄ (t)⟩.
∂t k̄

(2.41)

In summary, the Cb-TDHFB equations consists of Eqs. (2.41), (2.35), and (2.39). To derivethese equations from the TDHFB equations, we have assumed the diagonal property of the pairpotential, Eq.(2.40).

2.3

Properties of the Cb-TDHFB equations

2.3.1

Gauge invariance

The ηk (t) and ηk̄ (t), in Eqs.(2.39) and (2.41), must be real to conserve the orthonormal property,
however, they are arbitrary. This is related to the phase degrees of freedom of the canonicalstates. The Cb-TDHFB equations, (2.35), (2.39) and (2.41), are invariant with respect to thefollowing gauge transformations with arbitrary real functions, θk (t) and θk̄ (t).
|ϕk ⟩ → eiθk (t) |ϕk ⟩−i(θk (t)+θk̄ (t))κk → e

κk

and
and
18

|ϕk̄ ⟩ → eiθk̄ (t) |ϕk̄ ⟩−i(θk (t)+θk̄ (t))∆k → e(2.42)∆k(2.43)

2.3. Properties of the Cb-TDHFB equations
simultaneously withηk (t) → ηk (t) +dθkdt

and ηk̄ (t) → ηk̄ (t) +

dθk̄
.
dt

The phase relations of Eq.(2.43) are obtained from Eqs.(2.33) and (2.36).

2.3.2

Conservation laws

Orthonormality of canonical states
Apparently, Eq.(2.41) conserves the orthonormal property of canonical states, as far as ηk arereal.
∂
i ⟨ϕk (t)|ϕl (t)⟩ = ⟨ϕk (t)|{(h(t) − ηl (t)) − (h† (t) − ηk (t))}|ϕl (t)⟩ = 0.
(2.44)∂t
Here, we assume ⟨ϕk (t)|ϕl (t)⟩ = δkl at time t.
Average particle number
The average particle number also conserves because
i

∑
d
d ∑
ρk (t) = 2(κk (t)∆∗k (t) − κ∗k (t)∆k (t)) = 0,
N (t) = 2i
dt
dt
k>0

(2.45)

k>0

where we used the expression of the pairing energy, Eq.(2.71), for the last equation.
Average total energy
Time variation of the energy functional E[ρ, κ] can be divided into two: dE/dt = dE/dt|ρ +dE/dt|κ . The variation of energy associated with the normal-density ﬂuctuation isi∑ dρk∑ ∂E dρµνdE=i=i(ϵk (t) + ϵk̄ (t)),
dt ρ∂ρµν dtdtµν(2.46)

k>0

where ϵk (t) = ⟨ϕk (t)|h(t)|ϕk (t)⟩. This equation has an intuitive physical interpretation. Theenergy carried by a canonical state |ϕk ⟩ is ϵk (t)×ρk . If the occupation probability is ﬁxed during
the time evolution, the right-hand side of Eq.(2.46) vanishes. This corresponds to the cases suchas the TDHF and its extension with ﬁxed BCS occupation probabilities. In the TDHFB, theenergy variation in Eq.(2.46) transfers from/to the pairing energy. In fact, time variation due
to the pairing tensors-produce ﬂuctuation produces() ∑1 ∑ ∂E dκµνdE
∂E dκ∗µν=ii+ ∗=(κ∗k ∆k − κk ∆∗k )(ϵk (t) + ϵk̄ (t)),(2.47)dt κ
2 µν ∂κµν dt∂κµν dtk>0

where Eq.(2.40) is used. Because of Eq.(2.35), two contributions of Eqs.(2.46) and (2.47) always
cancel and the total energy is conserved. This is natural because the Cb-TDHFB equations
satisfy the TDHFB equations, (2.13) and (2.14), for which the conservation of the total energy
in TDHFB is well known [4].
19Chapter 2. Formulation

2.3.3

Stationary solution

When we assume that all the canonical states are eigenstates of the time-independent singleparticle Hamiltonian h0 .
|ϕk̄ (t)⟩ = |ϕ0k̄ ⟩eiθk̄ (t) ,|ϕk (t)⟩ = |ϕ0k ⟩eiθk (t) ,h0 |ϕ0k ⟩ = ϵ0k |ϕ0k ⟩,h0 |ϕ0k̄ ⟩ = ϵ0k |ϕ0k̄ ⟩,(2.48)(2.49)

where |ϕk̄ ⟩ = T |ϕk ⟩ have the same eigenvalues ϵ0k as |ϕk ⟩. Here, dθk /dt = i⟨∂ϕk /∂t|ϕk ⟩ anddθk̄ /dt = i⟨∂ϕk̄ /∂t|ϕk̄ ⟩ are arbitrary real functions of t. κk (t) and ∆k (t) should have a common
time-dependent phase associated with the chemical potential λ as e−2iλt . In addition to this,according to their deﬁnition, Eqs.(2.33) and (2.36), they have the following additional phases
connected with the phases of the canonical states.
κk (t) = κ0k exp{−i(2λt + θk (t) + θk̄ (t))},(2.50)∆0k exp{−i(2λt + θk (t) + θk̄ (t))},(2.51)∆k (t) =The stationary case of Eq.(2.35), dρ0k /dt = 0, indicates that κ0k and ∆0k have the same arguments
to make κk (t)∆∗k (t) real. If all the pairing matrix elements are real, we can choose both κ0k and∆0k are real. Then, Eq.(2.39) reduces to2(ϵ0k − λ)κ0k + ∆0k (2ρ0k − 1) = 0.
This is consistent with the ordinary BCS result3 .


0 −λ
ϵ
1
k
1 − √
,
ρ0k =
2
0
0
2
2
(ϵ − λ) + (∆ )
k

κ0k =

2.3.4

(2.53)

k

∆0k

1
√
.
2 (ϵ0 − λ)2 + (∆0 )2k(2.52)(2.54)k

Small amplitude limit and the Nambu-Goldstone modes

It is known that the small-amplitude approximation for the TDHFB around the HFB groundstate is identical to the QRPA. In the QRPA, when the ground state (HFB state) breaks continuous symmetry of the Hamiltonian, the Nambu-Goldstone modes appear as the zero-energy
modes. In this section, we show that this is also true for the small-amplitude limit of the
Cb-TDHFB.
The ground state is given by |ϕ0k ⟩, |ϕ0k̄ ⟩, κ0k , and ρ0k which satisfy Eqs.(2.49) and (2.52).
}∫t{Extracting trivial phase factors, ξk (t) ≡ 0 ηk (t′ ) − ϵ0k dt′ , we express the time-dependent
quantities as follows:{}{}|ϕk (t)⟩ = |ϕ0k ⟩ + |δϕk (t)⟩ eiξk (t) , |ϕk̄ (t)⟩ = |ϕ0k̄ ⟩ + |δϕk̄ (t)⟩ eiξk̄ (t) , (2.55){}{}
κk (t) = κ0k + δκk (t) e−i{ξk (t)+ξk̄ (t)+2λt} , ∆k (t) = ∆0k + δ∆k (t) e−i{ξk (t)+ξk̄ (t)+2λt} , (2.56)
ρk (t) = ρ0k + δρk (t), h(t) = h0 + δh(t), (2.57)
3

BCS equation is given in Appendix.C

20

2.3. Properties of the Cb-TDHFB equations
Substituting these into Eqs.(2.41), (2.35), and (2.39), they lead to the following equations in thelinear order with respect to the ﬂuctuation.
∂|δϕk (t)⟩ = (h0 − ϵ0k )|δϕk (t)⟩ + δh(t)|ϕ0k ⟩, (k ↔ k̄),(2.58)∂t∂0∗i δρk (t) = ∆0∗(2.59)k δκk (t) + κk δ∆k (t) − c.c.,∂t∂i δκk (t) = 2(ϵ0k − λ)δκk (t) + (2ρ0k − 1)δ∆k (t) + 2∆0k δρk (t).
(2.60)∂t
When these ﬂuctuating parts have speciﬁc oscillating frequency ω, they correspond to the normal
modes. The zero-energy modes correspond to stationary normal-mode solutions with ω = 0.
iTranslation and rotation
When the HFB ground state spontaneously violates the translational (rotational) symmetry,⃗ which transform the ground state into a new state but keep thethere are generators, P⃗ (J),energy invariant. Here, let us denote one of those hermitian generators, S. The transformation
with respect to the generator S with real parameter α leads to|ϕ0k ⟩ → |ϕ0k (α)⟩ = eiαS |ϕ0k ⟩ (k ↔ k̄),h0 → h0 (α) = e

iαS

−iαS

h0 e

,

(2.61)(2.62)

with ρk (α) = ρ0k , κk (α) = κ0k , ϵk (α) = ϵ0k , and ∆k (α) = ∆0k . These transformed quantities
should also satisfy Eq.(2.49).
(h0 (α) − ϵ0k )|ϕ0k (α)⟩ = 0.
(2.63)In the linear order with respect to the parameter α, we haveiα(h0 − ϵ0k )S|ϕ0k ⟩ + iα[S, h0 ]|ϕ0k ⟩ = 0.

(2.64)

Equation (2.64) means that |δϕSk ⟩ ≡ iαS|ϕ0k ⟩ and δhS ≡ iα[S, h0 ] correspond to a normal-mode
solution with ω = 0 for Eq.(2.58). δρSk = 0, δκSk = 0, and δ∆Sk = 0 also satisfy Eqs.(2.59)
and (2.60). Therefore, the Nambu-Goldstone modes related to the spontaneous breaking of
the translational and rotational symmetries become zero-energy modes in the small-amplitude
Cb-TDHFB equations.
Pairing rotation
When the ground state is in the superﬂuid phase, we have κ0k ̸= 0 at least for a certain k. The
ground state can be transformed into a new state by operation of eiθN where N is the number
operator. This transformation changes the phase of κk and ∆k but keeps the other quantities
invariant.
δκN
k
δρN
k

= e2iθ κ0k − κ0k ≈ 2iθκ0k ,= δhN = 0,

2iθ 0
0
0
δ∆N
k = e ∆k − ∆k ≈ 2iθ∆k ,N|δϕNk ⟩ = |δϕk̄ ⟩ = 0.

(2.65)(2.66)

Using Eq.(2.52), it is easy to see that these quantities correspond to an ω = 0 solution ofEqs.(2.58), (2.59), and (2.60). Thus, the pairing rotational modes appear as the zero energy
modes as well.
21

Chapter 2. Formulation
Particle-particle (hole-hole) RPA
The Cb-TDHFB equation (2.58) seems to be independent from the rest of equations, (2.59)
and (2.60), at ﬁrst sight. However, this is not true in general, because δ∆k (t) depend on|δϕk (t)⟩ and δh(t) depends on δρk (t). In contrast, when the ground state is in the normal phase(κ0k = ∆0k = 0), δ∆k (t) becomes independent from |δϕl (t)⟩, and we have δρk (t) = 0. This meansthat the particle-hole (p-h) channel is exactly decoupled from the particle-particle (p-p) andhole-hole (h-h) channels. It is well-known that TDHF equations in the small-amplitude limit,(2.58), reduce to the RPA equation in the ph-channel [3, 4, 59]. Thus, we here discuss propertiesof the p-p and h-h channels.
The p-p and h-h dynamics are described by the following equations.
i∂δκk (t) = 2ϵ0k δκk (t) ± δ∆k (t),∂t(2.67)

where the sign + (−) is for hole (particle) orbitals, and we omit the chemical potential λ. For thep-p channel (ω = EN +2 − EN ), a normal mode with frequency ω is described by δκp = Xp e−iωtfor particle orbitals (|p| > N/2) and δκh = −Yh e−iωt for hole orbitals (|h| ≤ N/2). For the h-hchannel (ω = EN −2 − EN ), it is described by δκh = Xh e−iωt for hole orbitals (|h| ≤ N/2) andδκp = −Yp e−iωt for particle orbitals (|k| > N/2). Equation (2.67) can be rewritten in a matrix
form as
( 0
)( )
(
)( )
2ϵp δpp′ + v̄pp̄p′ p̄′−v̄pp̄h′ h̄′
Zp′
1 0
Zp
,
(2.68)=ω−v̄hh̄p′ p̄′−2ϵ0h δhh′ + v̄hh̄h′ h̄′
Zh′
0 −1
Zh
where Zp = Xp (Zp = Yp ) and Zh = Yh (Zh = Xh ) for the p-p (h-h) channel. This is equivalent
to the p-p and h-h RPA in the BCS approximation [3].

2.4

Cb-TDHFB equations with a simple pairing energy functional and gauge condition

2.4.1

Pairing energy functional

Normally, the pairing energy functional is bi-linear with respect to κµν and κ∗µν . For instance,
when it is calculated from the two-body interaction, it is given byEκ (t) =∑vµν,ρσ κ∗µν (t)κρσ (t).

(2.69)µ,ν,ρ,σThus, the pairing energy can be also written as1∑ ∗1∑
κµν (t)∆∗µν (t) =κ (t)∆µν (t)2 µν
2 µν µν∑∑= −κk (t)∆∗k (t) = −κ∗k (t)∆k (t).

Eκ (t) =

k>0

k>0

22

(2.70)(2.71)

2.4. Cb-TDHFB equations with a simple pairing energy functional and gauge condition
For numerical calculations in the present study, we adopt a schematic pairing functional ina form of∑∑∑Eg (t) = −Gkl κ∗k (t)κl (t), = −κ∗k (t)∆k (t), ∆k (t) =Gkl κl (t),(2.72)
k,l>0

k>0

l>0

where Gkl is a hermitian matrix. This pairing functional produces a pair potential which isdiagonal in the canonical basis. This is consistent with the approximation of Eq. (2.40). However,the functional violates the gauge invariance (2.43), because∑∑
Gkl e−iθl +θl̄ κl (t) ̸= e−iθk +θk̄
Gkl κl (t).
(2.73)
l>0

k>0

The violation comes from the fact that the ∆k (t) in this schematic deﬁnition no longer holdthe correct phase relation to canonical states (k, k̄), according to the deﬁnition of Eq.(2.36).
∂ϕk̄
k
Therefore, we require the gauge condition of ⟨ ∂ϕ
∂t |ϕk ⟩ = ⟨ ∂t |ϕk̄ ⟩ = 0, so as to minimize thephase change of canonical states. This means that we choose the gauge parameters ηk (t) asηk (t) = ϵk (t) = ⟨ϕk (t)|h(t)|ϕk (t)⟩,2.4.2ηk̄ (t) = ϵk̄ (t) = ⟨ϕk̄ (t)|h(t)|ϕk̄ (t)⟩.

(2.74)

Properties of Cb-TDHFB equations with Eg

The Cb-TDHFB equations with the simple pairing functional (2.72) keep the following desired
properties, if we adopt the special gauge condition (2.74).
1. Conservation law(a) Conservation of orthonormal property of the canonical states(b) Conservation of average particle number(c) Conservation of average total energy
2. The stationary solution corresponds to the HF+BCS solution.
3. Small-amplitude limit(a) The Nambu-Goldstone modes are zero-energy normal-mode solutions.
(b) If the ground state is in the normal phase, the equations are identical to the particlehole, particle-particle, and hole-hole RPA with the BCS approximation.
Among these properties, 1(a) and 1(b) do not depend on the choice of the gauge, however, theother properties are guaranteed only with the special choice of gauge (2.74).

23

Chapter 3

Application of Cb-TDHFB tolinear-response calculations
3.1

Determination of the pairing strength g

In numerical calculations, we start from the HF+BCS calculation for the ground state. The
pairing energy is calculated for the constant monopole pairing interaction with a smooth truncation for the model space. We follow the prescription given by Tajima et al. [60], which isequivalent to the following choice of Gkl of Eq.(2.72).
Gkl = gf (ϵ0k )f (ϵ0l ),(3.1)

with a constant real parameter g. The cut-oﬀ function f (ε), which depends on the ground-state
single-particle energies, is in the following form(f (ε) =[
1 + exp

ε − ϵc
0.5 MeV] )−1/2

θ(ec − ε),(3.2)

with the cut-oﬀ energies
ϵc = λ̃ + 5.0 MeV,

ec = ϵc + 2.3 MeV,

(3.3)

where λ̃ is the average of the highest occupied level and the lowest unoccupied level in the HF
state. Here, the cut-oﬀ parameter ec is necessary to prevent occupation of spatially unlocalized
single-particle states, known as the problem of unphysical gas near the drip line. For neutrons,
if ec becomes positive, we replace it by zero.
To determine the pairing strength constant g for each nuclei, we again follow the prescriptionof Ref.[60] which is practically identical to the one in Ref.[61]. The pairing strength constantg in this thesis is determined in order to treat a wide range of nuclei. For this purpose, thedetermination of pairing strength g is introduced the continuous spectrum approximation. So,
25

Chapter 3. Application of Cb-TDHFB to linear-response calculations
we solve the following particle-number and gap equation expressed with level density D̄(ε),∫ ∞(ε − λ̄τ )2 D̄τ (ε)
Nτ =
dε √
,
(3.4)
¯2
(ε − λ̄τ )2 + f 2 (ε)∆−∞∫ ∞f 2 (ε)D̄τ (ε)¯ = gτ ∆
¯
∆
dε √
,
(3.5)
¯2
2
(ε − λ̄τ )2 + f 2 (ε)∆−∞
where D̄τ (ε) is single-particle level density obtained in the Thomas-Fermi approximation andindex τ means neutron and proton system. D̄τ (ε) forms
1D̄τ (ε) = 2
2π

(

∫
dr

2m∗τ (r)~2)3/2(ε − Vτ (r))1/2 Θ(ε − Vτ ),(3.6)

where m∗τ is eﬀective mass and Vτ means the central potential of ph-part. Eq.(3.4) determines the
Fermi level λ̄τ for the density of state D̄τ , while Eq.(3.5) decide the strength gτ with the empirical¯ = 12A−1/2 MeV in each nucleon system. For light nuclei(A < 50), andformula for pairing gap ∆when the calculated value exceeds 0.6 MeV, we replace gτ with 0.6 MeV because the resulting gτ
becomes apparently too strong. This value also obey the Ref.[60]. The pairing force strengths
Gkl are calculated for the ground state and kept constant during the time evolution. We deﬁne
the state-independent pairing gap as follows:∑κk (t)f (ϵ0k ).
(3.7)∆(t) ≡ g
k>0

The gap parameter for each canonical pair of states, k and k̄, can be written as ∆k (t) =∆(t)f (ϵ0k ).

3.2

Linear-response calculation in real-time

The canonical states in the ground state deﬁne the initial state for time evolution. In order tostudy the linear response, we use an weak instantaneous external ﬁeld of Vext (⃗r, t) = −ηF (⃗r)δ(t)to start the time evolution. Here, F (⃗r) is a one-body ﬁeld, such as E1 operator with recoilcharges,{(N e/A)ri
for protons
Fi (⃗r) =,(3.8)−(Ze/A)ri for neutrons
where i = (x, y, z). We also study the isoscalar quadrupole (Q2K ) response withF (⃗r) = Q2K = √12(1 + δK0 ){ 2}
r Y2K (r̂) + r2 Y2−K (r̂) ,K = 0 and 2,(3.9)

where Q20 and Q22 modes are sometimes called β- and γ-vibration in low-energy part. Then,at time t = 0+, the canonical states are given byϕk (⃗r, σ; t = 0+) = eiηF (⃗r) ϕ0k (⃗r, σ),ϕk̄ (⃗r, σ; t = 0+) = eiηF (⃗r) ϕ0k̄ (⃗r, σ),26(3.10)

3.2. Linear-response calculation in real-timeand the BCS factors byρk (t = 0+) = ρ0k ,κk (t = 0+) = κ0k .

(3.11)The parameter η controls the strength of the external ﬁeld. In this paper, since we calculate thelinear response, it should be small enough to validate the linearity.
The strength function with respect to the operator F is calculated with the followings. We
set |Φ0 ⟩, |Φn ⟩ as a many-body ground state and excited state which are assumed eigenstate of a
Hamiltonian H. If we add an external ﬁeld which is weak enough to regard as perturbation to|Φ0 ⟩ at t = 0+, we can expand |Φ0 ⟩ with |Φn ⟩ as like,∑|Φ(0+)⟩ ≡ eiηF |Φ0 ⟩ ≃ |Φ0 ⟩ + iη|Φn ⟩⟨Φn |F |Φ0 ⟩.
nAnd, we can consider a time-evolution of many-body state |Φ(t)⟩ from |Φ(0+)⟩ with H,∑
e−iEn t |Φn ⟩⟨Φn |F |Φ0 ⟩,|Φ(t)⟩ = e−iHt |Φ(0+ )⟩ = e−iE0 t |Φ0 ⟩ + iηnwhere E0 , En are eigenvalue of H. Next, we can calculate the time-dependent expectation
valuef (t) with F and |Φ(t)⟩ as,()∑
f (t) = ⟨Φ(t)|F |Φ(t)⟩ = f (0) + iη|⟨Φn |F |Φ0 ⟩|2 e−iẼn t − eiẼn t ,nwhere Ẽn = En − E0 . Then, we get f (E) in the energy representation through the Fouriertransformation with a smoothing parameter Γ/2,∫ ∞dt e(iE−Γ/2)t (f (t) − f (0))
f (E) =
0
(
)
∑
i
i
= iη
|⟨Φn |F |Φ0 ⟩|2+,Ẽn − E − iΓ/2 Ẽn + E + iΓ/2nwhere Γ/2 needs for this integration in general. In a limit of Γ = 0, the relation between f (E)and the strength function of F can be describe with using Cauchy principal value P,
}
{
∑
(
)
P
P
2
+
+ iπ δ(Ẽn − E) − δ(Ẽn + E)
(3.12)
lim f (E) = −η|⟨Φn |F |Φ0 ⟩|Γ→0Ẽn − E Ẽn + E
nFinally, we can get a formula of the strength function with respect to the operator F [24].
S(E; F ) ≡∑|⟨Φn |F |Φ0 ⟩|2 δ(E − Ẽn ) = −n1
Imf (E),πηẼn > 0 ,(3.13)

In practical, we have introduced a smoothing parameter Γ which is set to 1 MeV throughout
the calculations in Sec.4.1, 4.3, 4.6. The formula can be obtained from the time-dependent
perturbation theory in the ﬁrst order with respect to η [24]. Note that the strength function
S(E; F ) is independent from magnitude of the parameter η as far as the linear approximation is27

Chapter 3. Application of Cb-TDHFB to linear-response calculations

60

η0
0.5η0
0.25η00.6
0.5

S(E;Q20) [ fm4/MeV ]S(E;Q20) [Arbitrary unit]

0.7

0.4
0.3
0.2
0.1
0

η0
0.5η0
0.25η0

50
40
30
20
10
0

0

5

10
15
20
Energy [ MeV ]

25

30

0

5

10
15
20
Energy [ MeV ]25

30

Figure 3.1: The isoscalar quadrupole strength

Figure 3.2: Same as Fig. 3.1 but they are normal-

distributions for prolate deformed nucleus 20 Newith using Cb-TDHFB in BKN interaction. Solidline shows the result in η0 , dashed line shows theresult in a half of η0 and dotted line shows the result in a quarter of η0 . The smoothing parameterof Γ = 1 MeV is used.

ized with respect to the amplitude η0 . The smoothing parameter of Γ = 1 MeV is used.

valid. Figure 3.1 and 3.2 show examples of the strength function when the linear approximation
is valid. Figure 3.1 shows the results of Q20 vibration of 20 Ne in BKN interaction1 given in
Ref.[13, 15]. Lines in ﬁgure 3.1 correspond with results of several η, solid line shows the resultin η0 , dashed line shows the result in a half of η0 and dotted line shows the result in a quarterof η0 . They have not been normalized with the amplitude η0 . Figure 3.2 shows the normalized
results. We can see that they give same results with valid amplitude. In the present study, we
adopt the value of η = 10−4 fm−1 for the E1 operator, and η = 10−3 fm−2 for the quadrupole
operator.

3.3

Details of numerical calculations

3.3.1

Three-dimensional coordinate space representation

We use the Cartesian coordinate-space representation for the canonical states, ϕk (⃗r, σ; t) =⟨⃗r, σ|ϕk (t)⟩ with σ = ±1/2. In this present study, the three-dimensional (3D) coordinate spaceis discretized in square mesh of ∆x = ∆y = ∆z = 0.8 fm in a sphere with radius of 12 fm
when the mass number of nucleus is less than 70, in the case of over 70 mass number, the 3D
coordinate space has the mesh of ∆x = ∆y = ∆z = 1.0 fm in a sphere with radius of 15 fm.
Thus, each canonical state is represented by ϕk (i, j, k, σ; t) with three discrete indexes for the
3D space.
1
The eﬀective interaction of Bonche, Koonin and Negele (BKN) had been originally proposed for analyzing
the dynamics of spin and isospin saturated nuclei. This eﬀective interaction is compose of δ-type two-body andthree-body terms, but it does not have spin-orbit interaction. So, BKN interaction is sometimes called simpliﬁed
Skyrme interaction.

28

3.3. Details of numerical calculations

3.3.2

Calculation for ground state

First, we need to obtain a static solution to construct an initial state for the time-dependentcalculation. The numerical procedure is as follows:
1. Solve Eq.(2.49) for occupied canonical states (|k| ≤ N/2) with ρk = 1, and construct the
HF Hamiltonian h0 [ρ], using the imaginary-time method [62] (See Appendix B for moredetail).
2. Calculate unoccupied canonical states ϕ0k (⃗r, σ) (|k| > N/2) below the energy cut-oﬀ ec ,using the imaginary-time method with h0 .
3. Solve the BCS equations [3] to obtain ρk and κk (See Appendix C B for more detail).
4. Update h0 [ρ] with new ρk , then solve Eq.(2.49) again with the imaginary-time method, tocalculate canonical states with ϵ0k < ec .
5. Back to 3. and repeat until convergence.
To solve Eq.(2.49), the imaginary-time-evolution operator for a small time interval ∆t is repeatedly operated on each single-particle wave function. After each evolution, the single-particle
wave functions are orthonormalized with the Gram-Schmidt∫ method from low- to high-energy
states.
We
addthe
constraints
forthecenterofmass,⃗rρ(⃗r) = 0, and the principal axis,∫
ri rj ρ(⃗r) = 0 (i ̸= j) for deformed nuclei (See Appendix B for more detail).

3.3.3

Energy density functional

In the present calculations, we adopt a Skyrme energy functional, ESky [ρ], with the parameterset of SkM* [63] which can reproduce experimental data of ground-state radii and multipole
moments. The functional contains both time-even and time-odd densities same as Ref.[64]. The
form of this energy density functional is,∫
ESky [ρ] =H(r) =

d⃗r H(⃗r),~2
τ + B1 ρ2 + B2 (ρ2n + ρ2p ) + B3 (ρτ − j 2 ) + B4 (ρn τn − j 2n + ρp τp − j 2p )2m
+B5 ρ∇2 ρ + B6 (ρn ∇2 ρn + ρp ∇2 ρp ) + B7 ρ2+α + B8 ρα (ρ2n + ρ2p )+B9 (ρ∇·J + j ·∇×ρ + ρn ∇·J n + j n ·∇×ρn + ρp ∇·J p + j p ·∇×ρp )+B10 ρ2 + B11 (ρ2n + ρ2p ) + B12 ρα + B13 ρα (ρ2n + ρ2p ),(3.14)

where the densities which have a isospin index n or p deﬁned by each single-particle wave
functions with isospin, and the density without isospin index means that of total (ρ ≡ ρn + ρp ).
29

Chapter 3. Application of Cb-TDHFB to linear-response calculations
The densities are deﬁned as,ρτ (⃗r) =∑|vk,τ |2 |ϕk,τ (⃗r, σ)|2 ,k,σττ (⃗r) =∑|vk,τ |2 |∇ϕk,τ (⃗r, σ)|2 ,

k,σ

1 ∑
|vk,τ |2 (ϕ∗k,τ (⃗r, σ)∇ϕk,τ (⃗r, σ) − ϕk,τ (⃗r, σ)∇ϕ∗k,τ (⃗r, σ)),2ik,σ∑∇·J τ (⃗r) = −i|vk,τ |2 (∇ϕ∗k,τ (⃗r, σ)×∇ϕk,τ (⃗r, σ ′ ))·⟨σ|σ|σ ′ ⟩,j τ (⃗r) =k,σ,σ ′ρτ (⃗r) =∑|vk,τ |2 ϕ∗k,τ (⃗r, σ)ϕk,τ (⃗r, σ ′ )⟨σ|σ|σ ′ ⟩,(3.15)

k,σ,σ ′

where σ means Pauli spin matrices. Using a generalized Skyrme eﬀective interaction v Skyrme iscomposed of density- and momentum-dependent two-body interaction with parameters t0 , x0 , t1 , x1 ,t2 , x2 , t3 , x3 , α and W :1vSkyrme (⃗r1 , ⃗r2 ) = t0 (1 + x0 Pσ )δ(⃗r1 − ⃗r2 ) + t1 (1 + x1 Pσ )(k′ 2 + k2 )δ(⃗r1 − ⃗r2 )2
1′+t2 (1 + x2 Pσ )k ·δ(⃗r1 − ⃗r2 )k + t3 (1 + x3 Pσ )ρα δ(⃗r1 − ⃗r2 )
6+iW (σ 1 + σ 2 )·k′ ×δ(⃗r1 − ⃗r2 )k,(3.16)

here Pσ is the spin exchange operator and k is the operator (1/2i)(∇1 − ∇2 ) acting on ket state
and k′ is the operator (−1/2i)(∇1 − ∇2 ) acting on bra state, we can express the coeﬃcients Bi
as:
1
1
1
t0 (1 + x0 /2), B2 = − t0 (1/2 + x0 ), B3 = (t1 (1 + x1 /2) + t2 (1 + x2 /2)),
2
2
4
1
1
= − (t1 (1/2 + x1 ) − t2 (1/2 + x2 )), B5 = − (3t1 (1 + x1 /2) − t2 (1 + x2 /2)),
4
16
1
1
1
=
(3t1 (1/2 + x1 ) + t2 (1/2 + x2 )), B7 = t3 (1 + x3 /2), B8 = − t3 (1/2 + x3 ),
16
12
12
1
1
1
1
1
= − W, B10 = t0 x0 , B11 = − t0 , B12 = t3 x3 , B13 = − t3 .
(3.17)
2
4
4
24
24

B1 =
B4
B6
B9

The parameters of SkM∗ are t0 = −2645.0, t1 = 410.0, t2 = −135.0, t3 = 15595.0, x0 = 0.09,~2x1 = x2 = x3 = 0.0, α = 1/6, W = 130.0. And the 2m= 20.73398 with using the nucleon
mass = 938.91897 MeV, ~c = 197.32705 MeV fm and e2 /~c = 1/137.03604, these parameters
are obeyed Ref.[65]. The pairing energy functional is added to this, to give the total energyfunctional, E[ρ, κ] = ESky [ρ] + Eg [κ].
30

3.3. Details of numerical calculations

3.3.4

Calculation of time propagation

To solve the Cb-TDHFB equations in real time, we use the simple Euler’s algorithm.
iϕk (t + 2dt) = iϕk (t) + {h(t + dt) − ϵk (t + dt)}ϕk (t + dt) × 2dt,iρk (t + 2dt) =iρk (t) + {κk (t + dt)∆∗k (t + dt) − c.c.} × 2dt,

(3.18)(3.19)iκk (t + 2dt) = iκk (t) + [κk (t + dt){ϵk (t + dt) + ϵk̄ (t + dt) − 2λ}+∆k (t + dt){2ρk (t + dt) − 1}] × 2dt.

(3.20)

Here, we insert the chemical potential in Eq.(2.39) which cancels a global time-dependent phaseat the ground state, e−2iλt , for κk and ∆k . To construct the states at the ﬁrst step of t = dt,
we use the fourth-order Taylor expansion of the time-evolution operator for the canonical states[24] and use the Euler’s method for ρk (dt) and κk (dt). The time step dt is 0.0005 MeV−1 . The
time evolution is calculated up to T = 10 MeV−1 .

31Chapter 4

Results and Discussion
In this chapter, we show many linear-response calculations with Cb-TDHFB. We apply the
Cb-TDHFB method to calculation of the strength functions for Carbon (C), Oxygen (O), Neon
(Ne), Magnesium (Mg), Silicon (Si), Sulfur (S), Argon (Ar) and Calcium (Ca) isotopes andseveral heavy nuclei.
At ﬁrst, Section 4.1 shows comparison with QRPA calculations in order to conﬁrm that CbTDHFB is correct or not, as a time-dependent scheme including pairing correlation. Section 4.2
shows the eﬀect of spin-orbit and Coulomb residual interaction for light deformed nucleus (24 Mg)
and heavy nucleus (208 Pb and 154 Sm). Section 4.3 shows typical strength functions for spherical,
prolate and oblate nucleus in isovector dipole and isoscalar quadrupole excitations, and eﬀects
of pairing correlations in each excitation states. Section 4.4 shows the computational cost ofsmall-amplitude Cb-TDHFB calculation comparing with a recent deformed QRPA calculationfor 172 Yb. We show also the results of 236−240 U which are heaviest isotopes in our results.
Section 4.5 shows comparisons our results with experimental data. There are results of lightnuclei (16 O, 24,26 Mg and 40 Ca) and heavy nuclei (90 Zr, 208 Pb and Sm isotopes) for photo-nuclear
reaction cross sections. We mention an one of properties of E1 strength distributions with
Skyrme interaction. The ground-state data of nuclei presented in above sections are printed inTable 4.1. Table 4.1 shows ground-state properties; total energy, numbers of canonical orbitals
for neutron and proton, deformations, chemical potentials, and gap energies deﬁned by Eq.
(3.7). These nuclei show a variety of shapes (spherical, prolate and oblate), with and without
superﬂuidity. For nuclei in the superﬂuid phase with ∆ ̸= 0, the numbers of canonical orbitals,
Mτ are larger than neutron or proton number, but not signiﬁcantly diﬀerent. In the case with
∆ = 0, of course, we only calculate the occupied orbitals (Mn = N and Mp = Z). Therefore,
the numerical task of the Cb-TDHFB is in the same order as that of the TDHF. Note that,in the real-time calculation, the numerical cost is proportional to Mn + Mp . Section 4.6 shows
examples for systematic calculation with Cb-TDHFB. The number of nuclei presented in thissection is large (∼ 80), so their ground state data are printed in Appendix A we discuss theappearance of low-energy E1 strength which is often called “pygmy resonance”, and the eﬀects
of pairing correlation for this strength.
33

Chapter 4. Results and Discussion
125

Q20
Q22

100

S(E;Q2K) [fm4/MeV]

75

(a)

(b)

(c)

(d)

50
25
0
100
75
50
25
0
0

5

10 15 20 25 30 35

05

10 15 20 25 30 35 40

E [ MeV ]Figure 4.1: Calculated isoscalar quadrupole strength distribution for 34 Mg: (a) Cb-TDHFB
with time-independent spin-orbit and Coulomb potentials, (b) Cb-TDHFB, (c) Deformed
QRPA without the residual spin-orbit and Coulomb interactions[55]. and (d) Deformed QRPA
calculation[56]. The smoothing parameter of Γ = 1 MeV is used.

4.1

Comparison with QRPA calculations

We expect that the strength functions calculated in the present real-time approaches reproduce
those in the QRPA. This is strictly true if we solve the full TDHFB equations, however, since wesolve the Cb-TDHFB equations with the schematic pairing functional of Eq. (2.72) with (3.1), let
us ﬁrst show the comparison between our result and the HFB+QRPA calculations. The QRPA
calculations have been done with a computer program for axially deformed nuclei developed inRef. [55, 56], which diagonalizes the QRPA matrix of large dimensions in the quasi-particle
basis. This is based on the HFB ground state calculated in the two-dimensional coordinatespace representation with the Skyrme functional SkM* but with the density-dependent contact
interaction for the pairing channel. The space is truncated by the quasi-particle energy cut-oﬀ
of Ecut = 60 MeV and also by the cut-oﬀ for the magnetic quantum number of the quasiparticle angular momentum, Ωc = 19/2. In this QRPA calculation, the residual spin-orbit andCoulomb interactions are neglected. Thus, in order to make a comparison more meaningful, we
also neglect the time-dependence of these potentials in the Hamiltonian h(t), during the time
evolution.
In panels (a) and (c) of Fig. 4.1, we show the isoscalar quadrupole strength distributions(K = 0 and 2) for 34 Mg, calculated with (a) Cb-TDHFB and (c) QRPA calculation of K.Yoshidabased on [55]. The ground state has an axially symmetric prolate shape with ﬁnite pairing gaps
34

4.2. Full self-consistent calculation
for neutrons (Table 4.1). The HFB calculation with the contact pairing interaction for the panel(c) produces a deformation and average neutron pairing gap of β = 0.37 and ∆n = 1.7 MeV,
respectively. Note that a renormalization factor, which was used in Ref. [55], is set to be unity
in the present QRPA calculation. The peak energies in these calculations are approximately
identical, however, the height of the lowest peak is noticeably diﬀerent. We suppose that this isdue to the diﬀerence in the pairing energy functionals.
In panels (b) and (d) of Fig. 4.1, we show another comparison between the Cb-TDHFB
calculation and the QRPA calculation of Losa et al.[56] using the transformed harmonic oscillator
basis. Since this QRPA calculation includes all the residual interactions, it is compared with
the Cb-TDHFB calculation with the fully self-consistent time dependence. It turns out thatthe residual spin-orbit and Coulomb interactions slightly shift the giant quadrupole resonancehigher in energy, while they shift the lowest peak lower in energy. Actually, these shifts aremainly attributed to the residual spin-orbit interaction and the eﬀect of the residual Coulomb isvery small. The results in panels (b) and (d) well agree with each other, except for the height of
the lowest peak. We mention these eﬀects of spin-orbit and Coulomb residual interaction in next
section. Again, this may be due to the diﬀerent treatment of the pairing, because Ref. [55, 56]
also uses the contact pairing interaction. These comparisons indicate that the small-amplitude
Cb-TDHFB calculation well reproduces a fully self-consistent QRPA calculations. We would
like to mention that, for the isovector dipole excitations, the agreement is even better than theisoscalar quadrupole cases.

4.2

Full self-consistent calculation

Here, we mention eﬀects of spin-orbit and Coulomb residual interaction in the isovector dipole
and isoscalar quadrupole excitations. Sometimes, there are the RPA or QRPA calculations which
neglect or do not deal with directly these residual interaction. It is not often reasonable for numerical performance to include directly the eﬀects of these residual interaction with considering
the physical meanings caused by these residual interaction. Several (Q)RPA calculations include
the eﬀects of these residual interaction with a reasonable approximation, but in some case, theseapproximation is not good. We report about the “some case” which should need these residualinteraction, namely need fully self-consistent calculation.
First, let us show in Fig. 4.2 that the comparison results of the present Cb-TDHFB calculation with those of the RPA calculation. The fully self-consistent RPA calculation has been
done with the ﬁnite amplitude method (FAM) developed in Refs. [59, 66]. The same Skyrme
functional (SkM*) and the same model space have been used in these calculations. Since theground state of the 24 Mg nucleus is in the normal phase (∆ = 0), these two results should
be identical. This can be conﬁrmed in Fig. 4.2, which clearly demonstrates the accuracy of
our real-time method. From here, we mention the concrete eﬀects of spin-orbit and Coulomb
residual interaction. First, we show the eﬀects for isovector dipole mode of 24 Mg. In Fig. 4.3,
upper panel shows the result of fully self-consistent calculation (solid line labeled “Full”) andthe result without both spin-orbit and Coulomb residual interaction (dashed line labeled “Full(SO,C)”). Under panel show the comparison the result of Full-(SO,C) with the result without
only Coulomb residual interaction (dotted line labeled “Full-C”). We can see the eﬀects of35

Chapter 4. Results and Discussion

S(E;E1) [e2 fm2/MeV]

5

Cb-TDHFB
FAM

4

3

2

1

0
10

15

20

25

30

E [ MeV ]Figure 4.2: E1 strength distribution for 24 Mg calculated with the Cb-TDHFB (solid line) andwith the FAM[59, 66] (symbols). The smoothing parameter of Γ = 1 MeV is used.
these residual interaction. Upper panel shows that spin-orbit and Coulomb residual interaction
shift the peak position of GDR to higher energy only with 80 keV, and under panel shows thatthe eﬀects of Coulomb residual interaction is almost nothing. This eﬀects of Coulomb residualinteraction can be expected that 24 Mg is a light nucleus, so the eﬀects is small. At later, we
mention the eﬀects of Coulomb part for heavy nucleus.
Next, we mention the eﬀects of residual interactions for isoscalar quadrupole mode for 24 Mg.
Fig. 4.4 shows the results for Q20 -vibration mode of 24 Mg, and Fig. 4.5 shows the results forQ22 -vibration mode of 24 Mg. In the both Fig. 4.4 and Fig. 4.5, Solid line labeled “Full” means
the result of fully self-consistent calculation, in (b)-panel, dashed line labeled “Full-SO” means
the result without only spin-orbit residual interaction, in (c)-panel, dashed line labeled “Full-C”
means the result without only Coulomb residual interaction and in (d)-panel, dashed line labeled“Full-(SO,C)” means the result without both spin-orbit and Coulomb residual interaction. In theboth Q20 - and Q22 modes, we can see the signiﬁcantly eﬀects of spin-orbit residual interaction
and the eﬀects of Coulomb part is very small. The eﬀects of spin-orbit residual interaction
in isoscalar quadrupole mode shift the low-energy strength to lower energy region and giant
resonance to higher energy region at same time. The eﬀects of spin-orbit residual interaction
for isoscalar quadrupole mode are much stronger than those for isovector dipole mode in 24 Mg.
From these results (isovector dipole mode and isoscalar quadrupole mode), we can consider thatthe spin-orbit residual interaction has a dependence on an excitation mode.
Next, we show the eﬀects of residual interactions for heavy nucleus case. For example,
we mention the eﬀects for isovector dipole mode of the heavy spherical nucleus 208 Pb and theheavy prolate deformed nucleus 154 Sm. First, we mention the eﬀects of spin-orbit residualinteraction for isovector dipole mode of heavy nucleus. Fig. 4.6 shows photo-nuclear reaction
cross sections of 208 Pb computed from isovector dipole mode results. This value can be comparedwith experimental values and at later we show the comparison in Section 4.5. The ground state of
208 Pb is spherical with β = 0.0 and HF state with ∆ = ∆ = 0. Solid line labeled “Full” meansnp
36

4.2. Full self-consistent calculation

5

Full
Full-(SO,C)

4

S(E;E1) [e2fm4/MeV]

3
2
1
0

Full-(SO,C)
Full-(SO)

4
3
2
1
0
5

10

15

20

25

30

35

E [ MeV ]Figure 4.3: Comparison the fully self-consistent calculation with the calculation neglected residual interactions for isovector dipole mode of 24 Mg. Solid line means the fully self-consistentcalculation, dashed line means the calculation without spin-orbit and Coulomb residual interaction and dotted line means the calculation without only Coulomb residual interaction.
fully self-consistent calculation and dashed line labeled “Full-SO” means the result without spinorbit residual interaction. From Fig. 4.6, spin-orbit residual interaction moves the peak position
of GDR to higher region with 400 keV. This eﬀects is bigger than that of 24 Mg.
Next, we mention the eﬀects of Coulomb residual interaction for isovector dipole mode ofheavy deformed nucleus. Fig. 4.7 shows photo-nuclear reaction cross sections of 154 Sm computedfrom isovector dipole mode results. The ground state of 154 Sm is prolate deformed with deformation parameter γ = 0◦ and β = 0.32, and HF+BCS state with ∆n = 0.87, ∆p = 0.99MeV.
In the Fig. 4.7, we compare fully self-consistent calculation (solid line labeled “Full”) with theresults neglected Coulomb residual interaction (dashed line labeled “Full-C”). We can see theeﬀects of Coulomb residual interaction in the both K = 0, 1 vibration modes. Coulomb residualinteraction shifts the peak position of GDR to higher energy. The eﬀects of Coulomb part forisovector dipole modes in heavy nuclear system is much bigger than those in light nuclear system
as like 24 Mg.
Above discussions can be summarized two things. One is that spin-orbit residual interaction
37

Chapter 4. Results and Discussion

100
Full

Full

(b)

Full-SO

Full

Full

Full-C(d)

Full-(SO,C)

(a)

80

S(E;Q20) [fm4/MeV]

60
40
20
0

(c)

80
60
40
20
0
0

5

10

15

20

25

0

5

10

15

20

25

30

E [ MeV ]Figure 4.4: Fully self-consistent calculation and the calculations without spin-orbit or Coulomb
or both residual interaction for isoscalar quadrupole mode (K = 0 channel) of 24 Mg. (a)-panel
shows only the result of fully self-consistent calculation (solid line labeled “Full”). (b)-panel
shows the comparison “Full” with the result without only spin-orbit residual interaction (dashedline labeled “Full-SO”). (c)-panel shows the comparison “Full” with the result without only
Coulomb residual interaction (dashed line labeled “Full-C”). (d)-panel shows the comparison“Full” with the result without both spin-orbit and Coulomb residual interaction (dotted line
labeled “Full-(SO,C)”).

has dependence on excitation modes and the dependence is not small. The other is that theeﬀects of Coulomb residual interaction are negligible in light nuclear system, but in heavy
nuclear system the eﬀects become much bigger. In this section, we discuss the eﬀects of residualinteraction for giant resonance mainly. The discussion about low-energy part has remained yet,and about the eﬀects for isoscalar quadrupole mode of heavy nucleus has also remained yet.
In future work for this part, we want to suggest the auxiliary measures which can be able toinclude the eﬀects of residual interaction to (Q)RPA calculation easily. Lastly of this section,
we mention the “some case” which needs fully self-consistent calculation. At least, we need fully
self-consistent calculation for isoscalar quadrupole mode of deformed nuclei, because there are
the eﬀects of spin-orbit residual interaction for both low-energy and giant resonance parts.
38

4.3. Typical strength functions for spherical, prolate and oblate nucleus

100
Full

Full

(b)

Full-SO

Full

Full

Full-C(d)

Full-(SO,C)

(a)

80

S(E;Q22) [fm4/MeV]

60
40
20
0

(c)

80
60
40
20
0
0

5

10

15

20

25

0

5

10

15

20

25

30

E [ MeV ]Figure 4.5: Fully self-consistent calculation and the calculations without spin-orbit or Coulomb
or both residual interaction for isoscalar quadrupole mode (K = 2 channel) of 24 Mg. The
meanings of the lines are same as Fig. 4.4

4.3

Typical strength functions for spherical, prolate and oblate
nucleusIn this chapter, we show many and various strength functions. But these strength functions
can be categorized with some characters. We categorize these strength functions on their shape;
spherical, prolate and oblate shapes. Following sections, we mention the typical strength functions of spherical, prolate and oblate nucleus for isovector dipole and isoscalar quadrupole excitations, and eﬀects of pairing correlation in each excitations. We consider low- and high-energy
parts of strength functions. The low-energy (< 10 MeV) strength are important parts to investigate especially the structure of single-particle levels and eﬀects of pairing correlation. The
high-energy (∼ 20 MeV) strength are important to investigate collective motion of a nucleus.
This part is often called giant resonance (GR) part, for example, isovector dipole case is called
giant dipole resonance (GDR). GR are produced by nuclear-surface vibration basically.

4.3.1

Example of spherical nucleus (34 Si)

Here, we show the results of 34 Si which is spherical nucleus. For spherical nuclei case, the angularmomentum is a good quantum number. This means that we cannot distinguish the rotational
axis. Namely, spherical system has isotropic strength functions in dipole external ﬁelds for all
axes.
39

Chapter 4. Results and Discussion
400
Total-Full

Full
Full-(SO)800

Cross Section [mb]Cross Section [mb]

208

Pb

600

400

K=0-Full
K=1-Full154

Total Full-C

Sm

300

K=0 Full-C
K=1 Full-C

200

100

200
0
0

5
5

10

15

20

Figure 4.6: Isovector dipole strength distributions produced by the fully self-consistent calculation and the calculation without spin-orbit residualinteraction, for spherical nucleus 208 Pb. Solid line
shows fully self-consistent result, dashed line shows
the result without spin-orbit residual interaction.

10

15

20

25

Energy [ MeV ]E [ MeV ]Figure 4.7: Isovector dipole strength distributions produced by the fully self-consistent calculation and the calculation without Coulomb residual interaction, for prolate deformed nucleus 154 Sm.
Solid line shows fully self-consistent result, dashedline shows the result without Coulomb residual interaction.

As a concrete example, we show the result of 34 Si which has spherical ground state inboth cases with and without pairing correlation. Fig. 4.8 shows the isovector dipole strength
distributions S(E; E1) for 34 Si with while the pairing correlation. The solid line shows theresult with pairing correlation, dashed line shows the result without pairing correlation. Here,S(E; E1) is deﬁned as∑∑ ∑S(E; E1) =Si (E; E1) =|⟨n|Fi |0⟩|2 δ(E − Ẽn ),(4.1)
i=x,y,zi=x,y,z n

where the one-body operator Fi is given by Eq. (3.8). The K = 0 strength is Sz (E; E1) andK = 1 corresponds to Sx (E; E1) + Sy (E; E1).
The dipole strength is concentrated in a single peak at giant resonance in the both results TDHF
and Cb-TDHFB. There is little eﬀects of pairing correlation for GDR.
Next, we show the strength functions of isoscalar quadrupole mode for spherical nuclei. Fig.
4.9 shows results of TDHF with the ground state written by HF approximation, and Fig. 4.10
shows results of Cb-TDHFB with the ground state written by HF+BCS approximation. The
diﬀerence of total energy between HF+BCS and HF ground state is about 200 keV, the HF+BCS
ground state is more stable than HF one of 34 Si. In the each ﬁgures, solid line means the resultof Q20 -vibration and dashed line means the result of Q22 -vibration. These ﬁgures show same
strength distributions for Q20 - and Q22 -vibration in both Cb-TDHFB and TDHF. For spherical
nuclei case, they do not have the distinguishable axes, so the external ﬁeld of Q20 is equivalent
to Q22 in this system. Spherical nuclei should have same strength distributions for Q20 and
Q22 modes. In the results of Cb-TDHFB, the strength distributions for Q20 and Q22 modes
in low-energy have some diﬀerence. This reason is that the initial state of Cb-TDHFB has
some expectation value of Q22 . The order of “some” value is 10−5 fm2 . For the linear-response
results with real-time scheme, it is very important the accuracy of a ground state. In spherical40

4.3. Typical strength functions for spherical, prolate and oblate nucleus

15

Cb-TDHFB-IVD
TDHF-IVD

S(E;E1) [e2 fm2/MeV]

34

Si-IVD

10

5

0
0

5

10

15
20
25
Energy [ MeV ]

30

35

40

Figure 4.8: Calculated isovector dipole strength distributions for 34 Si. Solid line shows a resultof Cb-TDHFB calculation with pairing correlation, ∆n = 1.63 MeV. Dashed line shows a resultof TDHF calculation without pairing correlation. The smoothing parameter of Γ = 1 MeV isused on each calculations.
nuclei, the eﬀects of pairing correlation are small for the giant isoscalar quadrupole resonance.
A low-energy strength shifts to lower energy with 200 keV, but it is not so big change.

4.3.2

Example of prolate nucleus (34 Mg)

Here, we show the results of 34 Mg as an example of prolate deformed nuclei. 34 Mg has the HF
and HF+BCS ground state of prolate shape. The deformation parameter of both ground states
are γ = 0◦ , and HF state has βHF = 0.41 and HF+BCS state has βBCS = 0.37. Generally,
pairing correlation inﬂuences the shape of nuclei to become close to spherical. This 34 Mg casehas also this tendency. The diﬀerence of total energy between HF+BCS and HF ground stateis about 400 keV. HF+BCS ground state is more stable than HF one. First, we show theresults of isovector dipole mode for 34 Mg. Fig. 4.11 shows the results of TDHF with HF groundstate and Fig. 4.12 shows the results of Cb-TDHFB with HF+BCS ground state. Basically,
the isovector giant dipole resonance of axial symmetric deformed nuclei has separated peaks.
The vibrations along major and minor axes generate two giant resonances. We usually regard
that prolate deformed (axial-deformed) nuclei have two axes which a major axis is rotational
symmetric one and a minor axes is vertical to the major axis. Axial symmetric deformed nuclei
have K-quantum number in stead of angular momentum. K-quantum number is a projectedangular momentum to rotational symmetric axis. In the both Fig. 4.11 and Fig. 4.12, solid lines
mean the sum of strength distributions for major and minor axes, dashed lines labeled K = 0
mean the distribution produced by the vibration along major axis and dot lines labeled K = 1
mean that along minor axis. The energy peak of giant resonance of K = 0 is lower energy thanthat of K = 1. When we can regard these giant resonance produced by surface vibration, thecircumference of ellipsoid for major axis direction is longer than that for minor axis direction,
the vibrational cycle on the circumference for major axis is longer than that for minor axis. From
41

Chapter 4. Results and Discussion

120
100

S(E;Q2m) [e2 fm4/MeV]

2

Q20
Q22

100

34

Si - TDHF

80

4

S(E;Q2m) [e fm /MeV]

120

Q20
Q22

60
40
20

34

Si - Cb-TDHFB

80
60
40
20

0

0
0

5

10

15
20
Energy [ MeV ]

25

30

0

5

10

15
20
Energy [ MeV ]25

Figure 4.9: The isoscalar quadrupole strength

Figure 4.10: The isoscalar quadrupole strength

distributions for spherical nucleus 34 Si with using
TDHF. Solid line shows result of Q20 vibration,dashed line shows result of Q22 vibration. The
smoothing parameter of Γ = 1 MeV is used.

distributions for spherical nucleus 34 Si with using
Cb-TDHFB. The meanings of the lines are same as
Fig. 4.9. The smoothing parameter of Γ = 1 MeV
is used.

30

this reason, the energy peak of giant resonance of K = 0 is lower energy than that of K = 1.
This relation is one of character for prolate deformed nuclei. The eﬀects of paring correlation forisovector dipole mode of prolate deformed nuclei are also small, but we can see that the peaksof two giant resonances become slightly close, and the low-energy strengths extend slightly to
lower energy region.
Next, we show results of isoscalar quadrupole excitation for prolate deformed nucleus 34 Mg.
Fig. 4.13 shows the results of TDHF with HF ground state and Fig. 4.14 shows the result of
Cb-TDHFB with HF+BCS ground state. In the case of axial symmetric quadrupole deformed
nuclei, there is the relation for between the modes produced by external ﬁeld Q20 and Q22 which
are sometimes called Q20 - and Q22 -vibration. In the both Fig. 4.13 and Fig. 4.14, solid lines
mean result of Q20 -vibration and dashed lines mean result of Q22 -vibration. In prolate nuclei
case, the giant resonance of isoscalar quadrupole vibration are also two part as like the caseof isovector dipole mode, and the peak position of giant resonance produced by Q20 -vibration
is lower than that produced by Q22 -vibration. The eﬀects of pairing correlation for isoscalar
quadrupole mode of deformed nuclei appears in low-energy region. Form comparison Fig. 4.13
with Fig. 4.14, we can see that solid line has the new mode in lower than 3 MeV obviously.
There is no signiﬁcantly diﬀerence in the dashed line. This new mode can be expected pairing
vibration. Axial symmetric deformed nuclei have K-quantum number that we already havementioned, so in this axial symmetric deformed system, external ﬁeld Q20 can induce K = 0
states of nuclear system. The pair of nucleons with total angular momentum zero can be inducedby the external ﬁeld Q20 , but it is not induced by the external ﬁeld Q22 . The pairing vibration
mode appear in the energy less than a breaking energy of nucleon pair. This breaking energyof pair is characterized by gap energy as 2∆. Actually, the breaking energy of pair can beevaluated with gap energy in 34 Mg to 2∆n ≃ 3 MeV. In adding, the external ﬁelds of Q20 and
Q22 for spherical system only induce L = 2 state, so the pair vibration mode with K = 0 does
not appear.
42

4.3. Typical strength functions for spherical, prolate and oblate nucleus

6

6

TDHF-Total

Cb-TDHFB-Total
K=0
K=1

S(E;E1) [e2 fm2/MeV]

S(E;E1) [e2 fm2/MeV]

K=0
K=1

34

Mg-IVD

4

2

34

Mg-IVD

4

2

0

0
0

5

10

15
20
Energy [ MeV ]

25

0

30

5

10

15
20
Energy [ MeV ]25

Figure 4.11: The isovector dipole strength distri-Figure 4.12: The isovector dipole strength dis-bution for 34 Mg with HF ground state which does
not include pairing correlation. HF ground stateof 34 Mg is prolate deformed shape with β = 0.41.
Solid line means the sum of strengths for K = 0
and K = 1 channels. Dashed line indicates thestrength of K = 0 channel and dotted line indicatesthe strength of K = 1 channel. The smoothing parameter of Γ = 1 MeV is used.

tribution for 34 Mg with HF+BCS ground state
which includes pairing correlation, ∆n = 1.45 MeV.
HF+BCS ground state of 34 Mg is prolate deformed
shape with β = 0.37. The meanings of the lines
are same as Fig. 4.11. The smoothing parameterof Γ = 1 MeV is used.

4.3.3

30

Example of oblate nucleus (24 Ne)

Here, we show the results of 24 Ne as an example of oblate deformed nuclei. 24 Ne has the HF
and HF+BCS ground state of oblate shape. The deformation parameter of both ground states
are γ = 60◦ , and HF state has βHF = 0.20 and HF+BCS state has βBCS = 0.17. The shapeof nuclei closes to spherical by pairing correlation in this oblate shape case also. The diﬀerenceof total energy between HF+BCS and HF ground state is about 300 keV. HF+BCS groundstate is more stable than HF one. First, we show the results of isovector dipole mode for 24 Ne.
Fig. 4.15 shows the results of TDHF with HF ground state and Fig. 4.16 shows the results of
Cb-TDHFB with HF+BCS ground state. In the case of oblate deformed nuclei, the isovector
giant dipole resonance has also separated peaks. The vibrations along major and minor axes
generate two giant resonances. In this case also, we regard usually that oblate deformed nuclei
have two axes which a minor axis is rotational symmetric one and a major axes is vertical tothe minor axis. In the case of 24 Ne, however the peaks do not obvious separate ones, because its
deformation parameter is relatively small. In the both Fig. 4.15 and Fig. 4.16, solid lines mean
the sum of strength distributions for major and minor axes, dashed lines labeled K = 0 mean
the distribution produced by the vibration along minor axis and dot lines labeled K = 1 meanthat along major axis. The energy peak of giant resonance of K = 0 is higher energy than thatof K = 1. This character of GDR is opposite in prolate deformed nuclei. The eﬀects of paring
correlation for isovector dipole mode of oblate deformed nuclei are also small.
Next, we show results of isoscalar quadrupole excitation for oblate deformed nucleus 24 Ne.
Fig. 4.17 shows the results of TDHF with HF ground state and Fig. 4.18 shows the results of
Cb-TDHFB with HF+BCS ground state. In the both Fig. 4.17 and Fig. 4.18, solid lines mean
43

Chapter 4. Results and Discussion

100

34

34

Mg - Cb-TDHFB

S(E;Q2K) [e2 fm4/MeV]

80

60

2

4

S(E;Q2K) [e fm /MeV]

Mg - TDHF

100

Q20
Q22

40

20

0

Q20
Q22

80

60

40

20

0
0

5

10

15
20
Energy [ MeV ]

25

30

0

5

10

15
20
Energy [ MeV ]25

Figure 4.13: The isoscalar quadrupole strength

Figure 4.14: The isoscalar quadrupole strength

distributions for 34 Mg with using TDHF. The
meanings of the lines are same as Fig. 4.9. The
smoothing parameter of Γ = 1 MeV is used.

distribution for 34 Mg with using Cb-TDHFB. The
meanings of the lines are same as Fig. 4.9. The
smoothing parameter of Γ = 1 MeV is used.

30

result of Q20 -vibration and dashed lines mean result of Q22 -vibration. In oblate nuclei case,the peak position of giant resonance produced by Q20 -vibration is higher than that producedby Q22 -vibration. There appear the eﬀects of pairing correlation for isoscalar quadrupole modein low-energy region. Form comparison Fig. 4.17 with Fig. 4.18, we can see that solid line
has the new mode in lower than 1.5 MeV obviously. This new mode can be expected pairing
vibration also as like prolate case. There is no signiﬁcantly diﬀerence of peak position in thedashed line. K-quantum number is also good quantum number in oblate deformed nuclei, so
the external ﬁeld Q20 can induce the pair of nucleons with K = 0 but, the external ﬁeld Q22
does not. The breaking energy of pair can be evaluated with gap energy in 24 Ne to 2∆n ≃ 1.5
MeV. The character of this pairing vibration mode is consistent with prolate results.

44

4.3. Typical strength functions for spherical, prolate and oblate nucleus

6

6

TDHF-Total

Cb-TDHFB-Total
K=0
K=1

S(E;E1) [e2 fm2/MeV]

S(E;E1) [e2 fm2/MeV]

K=0
K=1

24

Ne-IVD

4

2

0

24

Ne-IVD

4

2

0
0

5

10

15
20
Energy [ MeV ]

25

30

0

5

10

15
20
Energy [ MeV ]25

Figure 4.15: The isovector dipole strength distri-Figure 4.16: The isovector dipole strength distri-bution for 24 Ne with HF ground state which does
not include pairing correlation. HF ground state of24
Ne is oblate deformed shape with γ = 60◦ and
β = 0.20. The meanings of the lines are same as
Fig. 4.11. The smoothing parameter of Γ = 1 MeV
is used.

bution for 24 Ne with HF+BCS ground state which
includes pairing correlation, ∆p = 0.74. HF+BCS
ground state of 24 Ne is oblate deformed shape with
γ = 60◦ and β = 0.17. The meanings of the lines
are same as Fig. 4.11. The smoothing parameterof Γ = 1 MeV is used.

100

24

24

Ne - Cb-TDHFB

S(E;Q2K) [e2 fm4/MeV]

80

60

2

4

S(E;Q2K) [e fm /MeV]

Ne - TDHF

100

Q20
Q22

40

20

0

30

Q20
Q22

80

60

40

20

0
0

5

10

15
20
Energy [ MeV ]

25

30

0

5

10

15
20
Energy [ MeV ]25

Figure 4.17: The isoscalar quadrupole strength

Figure 4.18: The isoscalar quadrupole strength

distributions for oblate deformed nucleus 24 Ne withusing TDHF. The meanings of the lines are same as
Fig. 4.9. The smoothing parameter of Γ = 1 MeV
is used.

distributions for oblate deformed nucleus 24 Ne withusing Cb-TDHFB. The meanings of the lines aresame as Fig. 4.9. The smoothing parameter ofΓ = 1 MeV is used.

45

30

Chapter 4. Results and Discussion

Table 4.1: Calculated ground-state properties of the nuclei presented in this chapter; total energy
ETotal , numbers of canonical orbitals for neutrons and protons Orbitals(n, p) which are chosen
in our energy cutoﬀ given in Chapter 3, quadrupole deformation parameters (β, γ) which are
given in Appendix D, pairing gaps (3.7) for neutrons and protons (∆n , ∆p ), chemical potentials
for neutrons and protons (λn , λp ). In the case of normal phase (∆ = 0), we deﬁne the chemical
potential as the single-particle energy of the highest occupied orbital, λn = ϵ0N and λp = ϵ0Z . In
the cases of 24 Ne,34 Mg and 34 Si, both situations whose ground states are written HF or HF+BCS
approximation, are presented for comparison between results of normal and superﬂuid phase.
When pairing gap energy is zero (normal phase), the number of used orbitals is equivalent tothat of HF single-particle state. The total energy, pairing gaps and chemical potentials are given
in units of MeV.
−ETotal127.798

Orbitals (n, p)( 16, 16 )

β
0.0

γ
−

∆n
0.0

∆p
0.0

−λn
13.53

−λp
10.25

BCS
24 NeHF

195.051
194.870( 20, 16 )( 14, 10 )

0.17
0.20

60◦
60◦

0.0
0.0

0.74
0.0

10.57
10.62

13.03
13.50

24 Mg
34 MgBCS
34 MgHF

197.299
264.071263.658( 12, 12 )( 28, 20 )( 22, 12 )

0.39
0.37
0.41

0◦
0◦
0◦

0.0
1.45
0.0

0.0
0.0
0.0

14.11
4.11
5.10

9.50
20.17
20.31

34 SiBCS
34 SiHF

287.935
287.486( 32, 20 )( 20, 14 )

0.0
0.0

−
−

1.63
0.0

0.0
0.0

6.56
7.87

16.94
16.99

40 Ca

341.559( 32, 28 )

0.0

−

0.0

0.0

14.31

7.47

90 Zr

790.943( 82, 50 )

0.0

−

0.0

1.75

12.75

7.18

144 Sm62
146 Sm62
148 Sm62
150 Sm62
152 Sm62
154 Sm62
172 Yb
70
236 U
92
238 U
92
240 U
92

1202.517
1217.104
1230.314
1244.523
1259.717
1274.550( 126, 82 )( 126, 82 )( 128, 82 )( 128, 84 )( 126, 90 )( 130, 92 )

0.0
0.0
0.12
0.20
0.29
0.32

−
−
0◦
0◦
0◦
0◦

0.0
0.93
0.90
0.94
1.04
0.87

1.97
1.94
1.62
1.45
1.08
0.99

12.26
7.18
7.21
7.30
7.63
7.23

4.41
4.96
5.41
6.01
6.62
7.15

1398.567( 146, 98 )

0.32

0◦

0.76

0.55

7.60

5.62

0.267
0.273
0.274

0◦

0.59
0.52
0.46

0.47
0.58
0.70

6.19
5.91
5.53

5.20
5.64
6.04

16 O
24 Ne-

1796.988
1809.022
1820.381( 200,132 )( 200,130 )( 202,130 )

46

0◦
0◦

4.4. Computational cost (172 Yb, 236,238,240 U)

4.4

Computational cost (172 Yb, 236,238,240 U)

In this section, we show the computational costs of Cb-TDHFB calculation for linear responseof several heavy nuclei. We have already mentioned that the computational cost of Cb-TDHFB
calculation will be much smaller than that of TDHFB calculation. But, there are no results of
TDHFB with modern eﬀective interaction in three-dimensional coordinate space representation.
So ﬁrst, we compare with fully self-consistent HFB+QRPA calculation for axially deformed
nucleus 172 Yb [57] in order to show that Cb-TDHFB can reduce the computational costs.
Fig. 4.19 shows the calculation result of isovector dipole mode of 172 Yb. Left panel shows theresult of J.Terasaki and J.Engel [57], and right panel shows the result of Cb-TDHFB calculation.
In our calculation, the ground state of 172 Yb is prolate deformed with the deformation parameters
γ = 0◦ and β = 0.32, and HF+BCS state with ∆n = 0.76, ∆p = 0.55 MeV. In their calculation,
the ground state of 172 Yb is prolate deformed with the deformation parameters γ = 0◦ and
β = 0.34, and axially deformed HFB state with ∆n = 0.77, ∆p = 1.25 MeV. In each panels,
solid line means a sum of results along a rotational axis (K = 0) and vertical to axial symmetricdirection axes (K = 1). Our Cb-TDHFB calculation produce almost same result with oneof Terasaki and Engel, except for width of total strength distribution. We mention the eachcomputational costs. The result (A) needs about 136,000 CPU hours in the parallel computation
with 10,000 CPU. Our result (B) needs only about 300 CPU hours in a single processor.

Figure 4.19: Isovector dipole strength distributions calculated by (A) axially deformed
HFB+QRPA of J.Terasaki and J.Engel [57] and by (B) our Cb-TDHFB, for prolate deformed
nucleus 172 Yb. The meanings of the lines are same as Fig. 4.11.
In addition, we show, as examples of other heavy nuclei, results of uranium isotopes. Fig.
4.20 show the photo-nuclear reaction cross section computed with isovector dipole modes of236,238,240 U. Dashed line means the result of K = 0 channel, dotted line means the result of
K = 1 channel and solid line means the sum of those. Black points and error bar show theexperimental data taken from [67]. These ground states of 236−240 U are all HF+BCS state andall prolate deformed. The calculation result of 238 U well agree with experimental data. The
computational costs of each result need about 400 CPU hours, and if we mention only real-time
47

Chapter 4. Results and Discussion
part, the costs are only about 40 CPU hours in a single processor. We can say with above results,
that Cb-TDHFB signiﬁcantly reduces the computational costs at least for linear response, canproduce reasonable results with the practical costs.

U

400

U

β=0.267∆n=0.586∆p=0.465

300
200
100

500

Total

238

K=0
K=1

Cross Section [mb]Cross Section [mb]

500

Total

236

400

U

Exp.
β=0.273∆n=0.517∆p=0.578

300
200
100

0

0
5

10

15

20

25

E [ MeV ]

Total

240

K=0
K=1

Cross Section [mb]

500

K=0
K=1

400

β=0.274∆n=0.463∆p=0.695

300
200
100
0

5

10

15

20

25

5

E [ MeV ]

10

15

20

25

E [ MeV ]Figure 4.20: Photo-nuclear reaction cross section for uranium isotopes with N = 144,146 and148. Their ground states are prolate deformed (γ = 0◦ ) and HF+BCS states. The each deformation parameters (β) and gap energies (∆n , ∆p ) are printed in each panels. The meanings ofthe lines are same as 4.11, the black points and error bar mean experimental data taken from[67]. The smoothing parameter of Γ = 2 MeV is used for the calculations.

4.5

Comparison with experiments (16 O, 24−26 Mg, 40 Ca, 90 Zr, 144−
154
Sm, 208 Pb)

This section, let us present photo-absorption cross sections in the GDR energy region (E =10 ∼ 30 MeV) for 16 O, 24−26 Mg, 40 Ca, 90 Zr, 144−154 Sm and 208 Pb together with experimental
data[68, 69, 70, 71, 72, 73]. We can calculate photo-absorption cross section σ easily fromisovector dipole (IVD) strength function S(E; E1), as likeσn (ε) =4π 2 e2Ẽ|⟨Φn |FIVD |Φ0 ⟩|2 δ(ε − Ẽ),
~c

(4.2)

where FIV D is E1 operator given by (3.8) and Ẽ means excited energy given in Section 3.2. We
can namely get the photo-absorption cross section with only to time excited energy Ẽ and some
coeﬃcient to S(E; E1). For 24 Mg, the peak energies of the GDR are underestimated by about
3 MeV. This disagreement has been already found in Ref.[66] for 24 Mg. The present calculation
also indicates that this underestimation of the GDR peak energy is also true for 26 Mg. The E1
strength distribution for 26 Mg is very similar to that in Fig.12 (bottom panel) in Ref. [56]. In
light nuclei (16 O and 24−26 Mg), the GDR energy is systematically underestimated in most of theSkyrme functionals[66], that seems to be true for nuclei with superﬂuidity. This underestimations
on the calculations with Skyrme functional has an alleviated tendency in heavier nuclei. The
GDR of 40 Ca become close to experimental data, and the mean energy of GDR for 90 Zr and
208 Pb well agree with experimental data. Next, we show the shape transition with respect48

4.5. Comparison with experiments (16 O, 24−26 Mg, 40 Ca, 90 Zr, 144−154 Sm, 208 Pb)
60
Total
Exp.

24

Cross Section [mb]

Mg

26

Mg

40

20

0
10

15

20

25

10

15

20

25

30

E [MeV]

Figure 4.21: Photoabsorption cross sections for 24,26 Mg. Experimental data (symbols) are takenfrom Ref. [68]. The smoothing parameter of Γ = 1 MeV is used for the calculations.

60

120

Total

Total

Exp.

Exp.

100

16

O

Cross Section [mb]Cross Section [mb]

50
40
30
20
10

40

Ca

80
60
40
20

0

0
5

10

15

20

25

30

35

5

E [ MeV ]

10

15

20

25

30

35

E [ MeV ]Figure 4.22: Photoabsorption cross sections for

Figure 4.23: Photoabsorption cross sections for

16

40

O. Experimental data (symbols) are taken from
Ref. [69]. The smoothing parameter of Γ = 1 MeV
is used for the calculations.

Ca. Experimental data (symbols) are taken from
Ref. [68]. The smoothing parameter of Γ = 1 MeV
is used for the calculations.

to neutron number and experimental data for lanthanide isotopes. There are many deformed
nuclei in lanthanide isotopes, and we can see the aspects of shape transition with respect toneutron number in this isotopes. As an example of the shape transition, we show the photonuclear reaction cross section of samarium isotopes from N = 82 to N = 92, samarium is
an one of lanthanide elements. 144,146 Sm has N = 82, 84 and they are spherical. N = 82 isa magic number. The photo-nuclear reaction cross sections of 144,146 Sm are in a single peakconcentrated at giant resonance. The experimental data of 144 Sm also indicate spherical shape.
This reason has been also mentioned in Section 4.3. 148−154 Sm have prolate shape for their
ground state. Their characters of prolate shape become clear gradually in photo-nuclear reaction
cross section with increasing neutron number. This character also appears in other lanthanide
isotopes (neodymium etc.)[2, 74]. The experimental data have also same character that the
GDR is not separated peak but the width of GDR becomes wide in 148,150 Sm, and the GDR
49

Chapter 4. Results and Discussion

400

800

Total

Total

90

300

Exp.

Zr

Cross Section [mb]Cross Section [mb]

Exp.

200

100

208

Pb

600

400

200

0

0
5

10

15

20

25

E [ MeV ]

5

10

15

20

25

E [ MeV ]Figure 4.24: Photo-nuclear reaction cross sec-

Figure 4.25: Photo-nuclear reaction cross sec-

tions for 90 Zr. Experimental data (symbols) aretaken from Ref. [71]. The smoothing parameter ofΓ = 1 MeV is used for the calculations.

tions for 208 Pb. Experimental data (symbols) aretaken from Ref. [72]. The smoothing parameter ofΓ = 1 MeV is used for the calculations.

is separated peak in 152,154 Sm. We must note that the results of 148,150 Sm do not have highaccuracy in low-energy region. This is because it is very diﬃcult with high accuracy to solve thereal-time calculation for 148,150 Sm which have some local minimums around true ground state.
In future work, we need to improve solving a true ground state in any case.

50

4.5. Comparison with experiments (16 O, 24−26 Mg, 40 Ca, 90 Zr, 144−154 Sm, 208 Pb)

500
400

144

Sm

150

146

Sm

152

Sm

300
200

Cross Section [mb]

100
500
0
400

Sm

300
200
100
500
0
148

400

Total

154

Sm

Sm

300

K=0
K=1

Exp.

200
100
0
5

10

15

20

5

10

15

20

25

E [ MeV ]Figure 4.26: Photo-nuclear reaction cross sections for 144−154 Sm. In each panels, solid, dashedand dotted lines indicate Cb-TDHFB calculation results with smoothing parameter Γ = 2 MeV,
and black points and error bar indicate experimental data taken from [73].

51Chapter 4. Results and Discussion

4.6

Pygmy resonance (C,O,Ne,Mg,Si,S,Ar,Ca isotopes)

In this section, we focus on the low-energy E1 strength which is often called pygmy dipole
resonance (PDR). The PDR is the one of hot topics in nuclear physics. There are two major
reasons on the PDR. One of the reason is a special structure of PDR which is indicated as a newcollective motion in neutron-rich nuclei, but has not been clariﬁed. The other reasons is the roleof PDR for nucleosynthesis. The PDR has big eﬀects for (γ, n) cross section which is important
for r-process [42, 43]. So, we discuss the appearance of the low-energy E1 strengths with using
the data of our systematic calculations. The ground state information and total E1 strength
distribution 1 of nuclei presented in this section, print to Appendix A, because a number ofthe nuclei is near 80, it is so large. In order to investigate how the low-energy pygmy strengthchanges with increasing the neutron number, we deﬁne the low-energy E1 ratio by m1 (Ec )/m1
with Ec = 10 MeV, where∫ E
m1 (E) ≡E ′ S(E ′ ; E1)dE ′ ,(4.3)
0

and m1 ≡ m1 (∞). Fig. 4.27 shows the part of strength related with this ratio for 26 Ne.
The ﬁlled part is deﬁned as “low-energy” part in this study of course, we need to conﬁrm thejustiﬁcation of this deﬁnition, but this deﬁnition gives us reasonable aspects of systematic lowenergy E1 strength at later reader can see2 . This low-energy peak structure has been recently
measured at RIKEN [76]. The calculated pygmy position is around 8 − 9 MeV, which agrees
with experimental data[76] and with the other QRPA calculations[77, 55]. The ratio is shownfor calculated even-even C, O, Ne, Mg, Si, S, Ar and Ca isotopes on some ﬁgures in this section.
First, we mention the ratio deﬁned by (4.3) for C and O isotopes. Fig. 4.28 shows theratios for C and O isotopes as function of neutron number, and Fig. 4.30 shows the ratio for C
and O isotopes as function of the diﬀerence between root mean square radii (Rr.m.s ) of neutron
and proton density distribution. In order to evaluate the width of neutron-skin, we use thediﬀerence of Rr.m.s in Fig. 4.30. In each Fig. 4.28 and 4.30, solid lines indicate the results of
Cb-TDHFB calculation and dotted lines indicate the results of HF+RPA calculation producedby FAM [59, 66]. Cross symbols show the results for C isotopes and ﬁlled circle symbols show theresults for O isotopes. These all HF+RPA calculations in this section are results of T.Inakura.
We must note that a computational space of this HF+RPA is lager than our space, they use
also the spherical box represented three-dimensional coordinate but their radius of the box is
15 fm. Some eﬀects from this diﬀerence of space-size occur in the case of heavier nuclei. If
we meet the eﬀect, we mention it. We discuss Fig. 4.28. We can see in both the ratio of C
and O isotopes, abrupt jumps between N = 14 and 16. This jump between N = 14 and 16
seems to be due to occupation of neutron s1/2 orbital. In the poor-neutron region, the ratio isalso increasing, it indicates the pygmy resonance produced by proton-rich nuclei, but we do notmention the detail of the appearance because our results for proton-rich side are not enough todiscuss the dependence on proton number. So, we only mention that there is the appearance
1

E1 strength distribution is equivalent to the strength function of isovector dipole mode.
Basically, we know that PDR appears around a separation energy [76], then we expect that the strengths
appear around a separation energy which can be evaluate with Bethe-Weisäcker mass formula [25, 75], so theaverage of separation energy is around 8.5 MeV. We do not have other signiﬁcant reason for this 10 MeV choice.
252

4.6. Pygmy resonance (C,O,Ne,Mg,Si,S,Ar,Ca isotopes)

S(E;E1) [e2 fm2/MeV]

7

Total

6

26

Ne

5
4
3
2
1
0
0

5

10

15 20 25 30 35 40
E [ MeV ]Figure 4.27: E1 strength distribution for 26 Ne, the ground state of 26 Ne is spherical and
HF+BCS state with ∆p =0.0 and ∆p =1.0 MeV. The black dotted line indicates the thresholddeﬁned as “low-energy” in this study.
of low-energy E1 strength in proton-rich nuclei also. There is a jump between N = 6 and 8 on
only C isotopes. We can expect that the reason for the jump is occupation of p1/2 orbital. We
can see the almost constant ratio in N = 10 − 14. This slow down of a increasing the ratio iscaused by start to occupy d3/2 orbital. We can interpret the eﬀect of d3/2 orbital for the ratio,as that the occupation of d3/2 orbital increase the total energy weighted sum (EWS), but does
not contribute for the low-energy E1 strength. We have discussed the ratio of pygmy resonancewith the orbitals deﬁned in spherical mean ﬁeld, however C isotopes of N = 10 − 14 are all
deformed nucleus and it is diﬃcult to separate GDR and PDR, so we can not conﬁrm that itis enough or not to explain the pygmy resonance with the orbitals. Fig. 4.30 is also a trial inorder to interpret the appearance of PDR as function of neutron-skin width, this trial is similar
to [78]. Actually, the ratio is in proportion to the width of neutron-skin has linear relation
roughly. We can see the relation between the constant ratio of C isotopes and the changes of
the neutron-skin width. In N = 10 − 14, the changes of the width are less than other neutron
number. From comparison with HF+RPA calculations, There is no remarkable eﬀects of thepairing correlation on the low-energy E1 strength for C and O isotopes.
Next, we show the ratios for Ne and Mg isotopes with Fig. 4.29 and 4.31. Fig. 4.29 shows
the ratio as function of neutron number and Fig. 4.31 shows that as function of diﬀerence
between Rr.m.s of neutron and proton. In each Fig. 4.29 and 4.31, solid lines indicate theresults of Cb-TDHFB calculation and dotted lines indicate the results of HF+RPA calculation
produced by FAM [59, 66]. Filled triangle symbols show the results for Ne isotopes and ﬁlled
square symbols show the results for Mg isotopes. In Fig. 4.29, we can see abrupt jumps between
N = 14 and 16, and between N = 20 and 22 on Ne and Mg isotopes. The ﬁrst jump between
N = 14 and 16 seems to be due to occupation of neutron s1/2 orbital. In contrast, the second
jump between N = 20 and 22 may be due to the onset of the deformation and the neutron
pairing. The second jump can not be seen in HF+RPA calculations. The ground states of Ne
53

Chapter 4. Results and Discussion
and Mg isotopes with N = 16, 18, 20 are spherical and then, those with N = 22, 24, 26, 28 are
deformed in our HF+BCS calculation. But, HF calculation yields the ground state with somedeformation for Ne and Mg isotopes in N = 16, 18, 20. The pairing correlation aﬀects the ratioof the low-energy E1 strength through the properties of ground state. In Fig. 4.31, we can not
see the eﬀects of pairing correlation for the ratio. The picture of the ratio with the diﬀerence of
Rr.m.s masks a information of nuclear shape, but we can see the same proportional relation on
the results of Ne and Mg isotopes with this picture.
Next, we show the ratios for Si and S isotopes with Fig. 4.32 and 4.34. Fig. 4.32 shows
the ratio as function of neutron number and Fig. 4.34 shows that as function of diﬀerence
between Rr.m.s of neutron and proton. In each Fig. 4.32 and 4.34, solid lines indicate theresults of Cb-TDHFB calculation and dotted lines indicate the results of HF+RPA calculation
presented by FAM. Filled diamond symbols show the results for Si isotopes and ﬁlled pentagon
symbols show the results for S isotopes. The pairing correlation aﬀects strongly the behaviorof the ratio of Si and S isotopes. In the Fig. 4.32, the ratio of Si isotopes increase gently
from N = 14 to 30 in both HF+RPA and Cb-TDHFB calculations. There is a small jump
between N = 28 and 30 in Cb-TDHFB calculation. This small jump is caused by changing
the shape of Si isotope of N = 30 to spherical. In contrary on the S isotopes, there is a
big diﬀerence between the results of HF+RPA and Cb-TDHFB calculation. The HF+RPA
calculation indicates an abrupt jump between N = 28 and 30, but the Cb-TDHFB calculationshows the smooth increase of the ratio for S isotopes from N = 16 to 34. The jump between
N = 28 and 30 in HF+RPA calculation seems to be due to occupation of neutron p3/2 orbital, thesmooth increase in Cb-TDHFB calculation is caused the pairing correlation which can describe
the continuous occupation probabilities around N = 28. In the Fig. 4.34, the relation between
the ratio and diﬀerence of Rr.m.s is not proportional but parabolic relation. We can not see theeﬀects of pairing correlation for the ratio in this Fig. 4.34 except for result of magic number.
Next, we show the ratios for Ar and Ca isotopes with Fig. 4.33 and 4.35. Fig. 4.33 shows theratio as function of neutron number and Fig. 4.35 shows that as function of diﬀerence between
Rr.m.s of neutron and proton. In each Fig. 4.32 and 4.34, solid lines indicate the results of
Cb-TDHFB calculation and dotted lines indicate the results of HF+RPA calculation presented
by FAM. Filled anti-triangle symbols show the results for Ar isotopes and star symbols show
the results for Ca isotopes. In the Fig. 4.33, there are some eﬀects of pairing correlation in the
behavior of the ratio of Ar isotopes. The smooth increase of the ratio appears in Ar isotopes
similarly with results of Si and S isotopes. We can expect that the reason is also the continuous
occupation probability of orbitals, namely the eﬀects of pairing correlation. For Ca isotopes,
there is no remarkable diﬀerence between the results of HF+RPA and Cb-TDHFB calculation.
There is an abrupt jump between N = 28 and 30 in both calculations. This jump between
N = 28 and 30 seems to be due to occupation of neutron p3/2 orbital. The results of CbTDHFB calculation behave same as results of HF+RPA calculation around N = 28, because
the pairing correlation does not work at N = 28. There is a slow down of the increase of
the ratio of low-energy E1 strength from Ca isotope with N = 34. This is expected to dueto the occupation of f5/2 orbitals. These aspects are similar to the case of C isotopes with
N = 10, 12, 14. In this Ca isotopes case, the occupation of f5/2 orbital increase the total EWS,
but does not contribute for the low-energy E1 strength. Namely, the denominator part of (4.3)
become large, but the numerator part of that increase a few, so the ratio does not increase. The
54

4.6. Pygmy resonance (C,O,Ne,Mg,Si,S,Ar,Ca isotopes)
space-size eﬀects appear in these Ca isotopes with over N = 34. The ratio decrease from N = 36
in Cb-TDHFB results but this reason is the space-size eﬀects that become diﬃcult to separate
between GDR and PDR. Fig. 4.35 indicates that f5/2 orbital does not contribute to Rr.m.s , thespacing of the symbols becomes small from N = 34. In addition, there is a proportional relation
between the ratio and the diﬀerence of Rr.m.s only while low angular momentum orbitals areoccupied.
We summarize the important things for the appearance of low-energy E1 strength fromthe above discussion about the results of C, O, Ne, Mg, Si, S, Ar and Ca isotopes. Whenthe low angular momentum orbital as like s or p orbital are occupied on the unstable nuclei,
the low-energy E1 strength appears. The ratio deﬁned by (4.3) is sensitive to the neutron
number especially around magic number, and the behavior of the ratio reﬂects the properties ofnuclear shape. The most eﬀect of pairing correlation for the appearance is continuous occupation
probabilities of orbitals around magic number, and there is the eﬀect of paring for the appearance
through the change of ground state. Pairing is also important to deﬁne the shape of nucleus. In
addition, we proposed the ratio of low-energy E1 strength, and try to investigate the appearancewith neutron number and diﬀerence between the Rr.m.s of neutron and proton density. The
picture of the ratio with neutron number is good for reﬂection of orbital and nuclear shapeproperties. In contrary, the picture with the diﬀerence of Rr.m.s gives simple aspects for a
relation between the appearance and nucleon-skin width, but only when around magic number.

55

Chapter 4. Results and Discussion

665m1(Ec) / m1 [%]m1(Ec) / m1 [%]

5
4
3
2

C
O

1

4
3
2

0
0

4

8

12

16

20

24

0

4

8

12

16

20

24

28

Neutron Number

Figure 4.28: The ratio deﬁned as Eq.(4.3) as

Figure 4.29: The ratio deﬁned as Eq.(4.3) asfunction of neutron number for C and O isotopes.
Solid line indicates the results of Cb-TDHFB calculation, dotted line indicates the results of HF+RPA
calculation produced by FAM based on [59, 66].
Cross symbols show the results for C isotopes andﬁlled circle symbols show the results for O isotopes.

function of neutron number for Ne and Mg isotopes.
The meanings of the lines are same as Fig. 4.28.
Filled triangle symbols show the results for Ne isotopes and ﬁlled square symbols show the results for
Mg isotopes.

6

6

5

m1(Ec) / m1 [%]

5
4
3
2

HF+RPA

-0.6

-0.4

-0.2
0
0.2
0.4
Difference of Rr.m.s [fm]

4
3
2

Ne
Mg

1

C
O

1
0
-0.8

HF+RPA

028

Neutron Number

m1(Ec) / m1 [%]

Ne
Mg

1

HF+RPA

0.6

0.8

0
-0.8

HF+RPA

-0.6

-0.4

-0.2
0
0.2
0.4
Difference of Rr.m.s [fm]0.6

Figure 4.31: The ratio deﬁned as Eq.(4.3) as

Figure 4.30: The ratio deﬁned as Eq.(4.3) asfunction of diﬀerence between Rr.m.s of neutron andproton density for Ne and Mg isotopes. The meanings of the lines are same as Fig. 4.28. Filled triangle symbols show the results for Ne isotopes andﬁlled square symbols show the results for Mg isotopes.

function of diﬀerence between Rr.m.s of neutron andproton density for C and O isotopes. The meaningsof the lines are same as Fig. 4.28. Cross symbols
show the results for C isotopes and ﬁlled circle symbols show the results for O isotopes.

56

0.8

6

6

5

5

m1(Ec) / m1 [%]m1(Ec) / m1 [%]4.6. Pygmy resonance (C,O,Ne,Mg,Si,S,Ar,Ca isotopes)

4
3
2

Si
S

1

8

12

16

20

24

28

32

36

3
2

Ar
Ca

1

HF+RPA

0

4

HF+RPA

0
40

44

8

12

16

Neutron Number

20

24

28

32

36

40

44

Neutron Number

Figure 4.33: The ratio deﬁned as Eq.(4.3) asfunction of neutron number for Si and S isotopes.
The meanings of the lines are same as Fig. 4.28.
Filled diamond symbols show the results for Si isotopes and ﬁlled pentagon symbols show the resultsfor S isotopes.

function of neutron number for Ar and Ca isotopes.
The meanings of the lines are same as Fig. 4.28.
Filled anti-triangle symbols show the results for Ar
isotopes and star symbols show the results for Ca
isotopes.

6

6

5

5

m1(Ec) / m1 [%]m1(Ec) / m1 [%]Figure 4.32: The ratio deﬁned as Eq.(4.3) as

4
3
2

Si
S

1
0
-0.4

0
0.2
0.4
Difference of Rr.m.s [fm]

3
2

Ar
Ca

1

HF+RPA

-0.2

4

0
-0.4

0.6

HF+RPA

-0.2

0
0.2
0.4
Difference of Rr.m.s [fm]0.6

Figure 4.34: The ratio deﬁned as Eq.(4.3) as

Figure 4.35: The ratio deﬁned as Eq.(4.3) asfunction of diﬀerence between Rr.m.s of neutron andproton density for Si and S isotopes. The meaningsof the lines are same as Fig. 4.28. Filled diamondsymbols show the results for Si isotopes and ﬁlled
pentagon symbols show the results for S isotopes.

function of diﬀerence between Rr.m.s of neutron andproton density for Ar and Ca isotopes. The meanings of the lines are same as Fig. 4.28. Filled antitriangle symbols show the results for Ar isotopes
and star symbols show the results for Ca isotopes.

57Chapter 5

Summary and Future work
We have developed a new theory to study excited states and dynamics of atomic nuclei, which
we call the Canonical-basis time-dependent Hartree-Fock-Bogoliubov (Cb-TDHFB). The CbTDHFB is a simpliﬁed theory of the full TDHFB theory, treating pairing energy functional witha BCS-like approximation. Implementing the theory with a real-space and real-time method, we
can describe nuclear excitations and dynamics taking account of deformation eﬀects and pairingcorrelations, even for heavy nuclei with a reasonable computational costs. In the thesis, we ﬁrst
presented a derivation of the Cb-TDHFB theory and its practical implementation. We thenpresented examples of the linear response calculations using the Cb-TDHFB theory.
In Chapter 2, we derived the Cb-TDHFB equations starting with the full TDHFB equations,introducing a simpliﬁed form for the pairing energy functional. In the derivation, we start withthe ordinary TDHFB equation written in the time-dependent normal density matrix ρ(t) andthe pairing tensor κ(t). We then rewrite the equation in the canonical basis, which diagonalize
ρ(t). The derived TDHFB equation in the canonical basis has a very simple form. However, it is
not useful practically for a general case since it costs much to obtain the canonical basis itself at
each time. Solving the time-dependent equation for the canonical basis is as diﬃcult as solving
the full TDHFB equation. This diﬃculty disappears if we employ a simpliﬁed pairing energy
density functional, which is equivalent to that employed in the BCS approximation for groundstate calculations. With the simpliﬁed pairing functional, the evolution equation for the timedependent canonical-basis function becomes nothing but that of time-dependent Hartree-Fock(TDHF) so that the computational cost of Cb-TDHFB is almost the same as that of TDHF
calculation.
In Chapter 3, we shown the concrete form of pairing functional, the relation strength function with linear response and the procedure of linear-response calculation with time-dependent
scheme. We explained also the detail and setup of numerical calculation in present study.
In Chapter 4, we presented results of the Cb-TDHFB calculations for excited states of somespherical and deformed nuclei to show performance of the method.
In Section 4.3, we showed strength functions of isovector dipole and isoscalar quadrupole
excitations for some spherical and deformed nuclei. In the isoscalar quadrupole response, there
appear Q20 - and Q22 -vibrational modes which converge into a single identical mode in thespherical limit. In prolate nuclei, the excitation energy of the Q20 -vibration is lower than that of59

Chapter 5. Summary and Future work
Q22 -vibration, while opposite in the oblate nuclei. The pairing correlation eﬀect is rather weakfor the isoscalar quadrupole excitations in spherical nuclei. In deformed superﬂuid nuclei, theQ20 -vibrational mode induces the pair vibrational mode, while the Q22 -vibrational mode doesnot. This is because the K-quantum number is a good quantum number in axially symmetricnuclei and the pairing correlation is coupled with K = 0 mode of isoscalar quadrupole model.
In Section 4.1, we next compared results of Cb-TDHFB calculation with those of HFB+QRPA
which is equivalent to the small amplitude limit of the full TDHFB. We compared our resultswith axially-symmetric deformed HFB+QRPA calculations for isoscalar quadrupole excitationsof neutron-rich magnesium isotope 34 Mg. Our calculation show a good agreement with existing
HFB+QRPA calculations. However, we observed a substantial diﬀerence for the height of lowest peak. We consider that this diﬀerence for the lowest peak is caused by the diﬀerence of thepairing energy functionals.
In Section 4.2, our calculation is achieved fully self-consistent, employing the same energyfunctional as that in the ground state. We mentioned the importance of the self-consistency in
our Cb-TDHFB calculation. We conﬁrm the accuracy of our time-dependent scheme by comparing our results with other fully self-consistent calculation performed with the ﬁnite amplitude
method (FAM) for isovector dipole excitation of 24 Mg. Furthermore, we found that there aresome diﬀerences in the eﬀects of the spin-orbit and Coulomb residual interactions for each excitation mode. These residual interactions are found to shift the peak energies of isoscalar giant
quadrupole resonance to higher energy region and to move low-energy peaks to lower energy
region simultaneously. The eﬀect of Coulomb residual interaction is rather weak for the isoscalar
quadrupole strength distribution. In isovector dipole excitation, the spin-orbit residual interaction shifts the peak position of the isovector giant dipole response to higher energy region.
However, the eﬀect is much weaker than that in the isoscalar quadrupole case. The eﬀect ofresidual Coulomb force is again very small for isovector dipole excitation of 24 Mg. We also
examined the signiﬁcance of the full self-consistency for the isovector dipole strength functionsof heavy nuclei taking 208 Pb and 154 Sm nuclei as examples. The eﬀect of the spin-orbit residualinteraction was found to be signiﬁcant for 208 Pb as in the case of 24 Mg. The eﬀect of Coulomb
residual interaction is found to be much larger in the heavy nuclei.
In Section 4.4, in order to assess the computational cost of our Cb-TDHFB approach, weachieved calculations of 172 Yb and compared our computational costs with those by J. Terasaki
and J. Engel who achieved a fully self-consistent HFB+QRPA calculations for axially deformed
nuclei. The calculated isovector dipole strengths for 172 Yb coincide each other with high accuracy. Terasaki and Engel reported that their computational cost was about 136,000 CPU hours
in the parallel computation with 10,000 CPU. In our calculation, we need only 300 CPU hours ina single processor. We also achieved calculations for 236−240 U. It costs 400 CPU hours, showing
approximately linear scaling for the mass number. We also note that the calculated spectra
agree well with measurements. Thus, we may conclude that the Cb-TDHFB method shows aquite nice aspect in the computational costs and in the feasibility for large nuclei.
In Section 4.5, we next compared results of Cb-TDHFB calculation with experimental datafor photo-nuclear reaction cross section. We calculated the photo-nuclear reaction cross sectionof 16 O, 24,26 Mg, 40 Ca, 90 Zr, 208 Pb and Sm isotopes. From the results of 16 O, 24,26 Mg, and
40 Ca, we saw that our results underestimate the experimental data. This is a properties of thecalculation with Skyrme interaction which has been also reported by T.Inakura. But in heavy
60

nuclei (90 Zr, 208 Pb), the calculation results become to agree well with the experimental data
for the centroid energy. In Sm isotopes, we could see the shape transitional properties withincreasing neutron number from 82 to 92 in both experimental data and calculation results.
144,146 Sm are spherical, but from 148 Sm, the isotopes are gradually becoming prolate shape withrespect to neutron number. Then, experimental data has one concentrated giant resonance peakin 144,146 Sm, but from 148 Sm, the peak are becoming wide and at last becomes separated one at
152,154 Sm.
In Section 4.6, we presented an application of our approach for low-energy isovector dipole(E1) strength distribution below 10 MeV. In the low-energy E1 response, it has been well known
that there appears the so-called pygmy resonance. There are two major interests on the pygmy
resonance. It may be regarded as a new kind of collective excitations which appear in neutronrich nuclei. It may also aﬀect the (n, γ) cross section which is important in nucleosynthesis.
Recently, systematic calculations of the E1 strength function has been achieved in the
Hartree-Fock plus RPA (HF+RPA) approach. Although the HF+RPA approach does not include pairing correlation eﬀect, it has brought an interesting result for the systematic behavior of
the low-energy E1 strength distribution. Namely, the appearance of the low-energy E1 strength
is intimately related to the occupation of the low partial wave orbitals such as s and p orbitals. We have made a similar systematic calculation and investigated the eﬀect of the pairingcorrelation.
In Carbon (Z = 6) and Oxygen (Z = 8) isotopes, the low-energy E1 strengths are found toincrease in the neutron-rich isotopes when the neutron number exceeds N = 16. Since neutron
occupies s1/2 orbital for N = 16 and more nuclei, it is considered that the occupation of the s1/2
orbital is the origin of the low-energy E1 strength. This result suggests that the appearance ofthe pygmy resonance is intimately related to the occupation of low angular momentum orbitals.
We also ﬁnd a rise of the low-energy E1 strength at N = 8 for Carbon isotopes while it isalmost constant for isotopes in N = 10 − 14. The increase at N = 8 is considered to be causedby the occupation of p1/2 orbital. In Neon (N = 10) and Magnesium (Z = 12) isotopes, weagain ﬁnd that the low-energy E1 strengths increase beyond N = 16. We consider that it
originates from the occupation of s1/2 orbital, as in the Carbon and Oxygen isotopes. For theseisotopes, we ﬁnd a signiﬁcant eﬀect of the pairing correlation through the change of the groundstate deformation. In our approach, the ground state calculation is achieved by the HF+BCS
theory. In the HF+BCS theory, the ground state of Ne and Mg isotopes with neutron number
N = 16 − 20 are all spherical. The ground states are deformed for neutron number N ≥ 22.
However, in the previous HF+RPA calculation, isotopes with N = 16 − 20 are deformed since
the pairing correlation is not included. Our calculation shows a small but abrupt increase of
the low-energy E1 strength at N = 22 reﬂecting the appearance of the deformation, which theprevious HF+RPA calculation does not. For Silicon (Z = 14) and Sulfur (Z = 16) isotopes,
the pairing eﬀect on the low-energy E1 strength becomes more appreciable. In the HF+RPA
calculation showed an abrupt increase of the low-energy E1 strength at N = 30, which is causedby the occupation of the p3/2 orbital. However, in our Cb-TDHFB calculation, the change israther smooth reﬂecting a continuous change of the p3/2 occupation. In the cases of Argon(Z = 18), the change of low-energy E1 strength is also rather smooth reﬂecting a continuouschange of the p3/2 occupation. But in Calcium (Z = 20) isotopes, an abrupt increase of thelow-energy E1 strength appeared at N = 30 in our Cb-TDHFB calculation also, because the61Chapter 5. Summary and Future work
pairing does not work for 48 Ca which is a double magic nucleus. The low-energy E1 strengthof calcium isotopes show a special character from N = 36 which can be expected eﬀects of f5/2
orbital. The increase of the low-energy E1 strength becomes slow down. It is supposed that the
f5/2 orbital does not contribute to the low-energy E1 strength.
Finally, we mention a future direction of the present work. We consider that there are two
directions which the present Cb-TDHFB approach will be of signiﬁcance. One is the systematic investigation of the nuclear response properties for a whole mass region. The other is anapplication towards large amplitude collective motion such as fusion and ﬁssion dynamics.
In the present thesis, we have reported an investigation of the photo-nuclear reactions ofneutron-rich nuclei up to Z = 20. The photo-nuclear reaction cross sections of nuclei heavier
than Iron (Z = 26) are signiﬁcant since they are crucial to understand the nucleosynthesis in
the r-process.
Regarding the large amplitude collective motion, the present Cb-TDHFB approach will allow
us to investigate the eﬀect of pairing correlation. Microscopic descriptions of the large amplitude
collective motion such as fusion and ﬁssion phenomena have been a long-standing issue in nuclear
theory. The fusion and ﬁssion dynamics are important not only in fundamental sciences but also
applications such as nuclear energy generation and transmutation of radioactivity. In order to
execute these studies, we consider the present Cb-TDHFB theory provides a promising starting
point. In view of rapid development of computational techniques, parallel computations withmassively parallel supercomputer will be a key issue to realize the above mentioned studies.

62

Acknowledgement
I would like to express my gratitude to Prof. T.Nakatsukasa and Prof. K.Yabana for patient
guidance, encouragements and supports during my study. I would like to express my acknowledgments to Prof. Y.Hashimoto for his continuous discussion as one of my collaborators. I would
like to thank Dr. T.Inakura and Dr. K.Yoshida for their useful advice as my collaborators.
I am also grateful to Prof. K.Matsuyanagi, Prof. J.A.Maruhn, Prof. D.Lacroix, Prof.
J.Dobaczewski and Prof. V.Oberacker for their encouragements for my study. I want to express
my deep appreciation to my colleagues in Theoretical Nuclear Group in University of Tsukuba,
and my colleagues in Theoretical Nuclear Physics Laboratory RIKEN Nishina Center.
The numerical calculations were performed on the computers in Theoretical Nuclear Groupof University of Tsukuba.
Finally but not least, I would like to express to thank my families for their understanding,
supports and encouragements to me in all time.

Appendixes
A

Data Table and E1 strength distributions

This section presents the ground state data and E1 strength distributions in this study. Following
tables show data of ground state: total energy ETotal , numbers of canonical orbitals for neutronsand protons Orbitals (n, p), quadrupole deformation parameters (β, γ), pairing gaps (3.7) for
neutrons and protons (∆n , ∆p ), chemical potentials for neutrons and protons (λn , λp ). In thecase of normal phase (∆ = 0), we deﬁne the chemical potential as the single-particle energyof the highest occupied orbital, λn = ϵ0N and λp = ϵ0Z , and the number of used orbitals is thesame of particle number for instance that in the case of 16 C, it has (16, 10) canonical orbitals
but with no pairing gap of neutron, practically, we use (16, 6) orbitals in which proton orbitals
are equivalent to HF single-particle states. The (16, 10) orbitals are obtained only in our cut-oﬀ
energies Eq.(3.3). The total energy ETotal , pairing gap energy ∆ and chemical potential λ are
given in units of MeV.

A.1

Ground state of Carbon (Z = 6) isotopes

8C
10 C
12 C
14 C
16 C
18 C
20 C
22 C

−ETotalOrbitals (n, p)

β

γ

∆n

∆p

−λn

−λp

31.15362.33993.567108.606
115.600
122.596128.954133.426( 6, 6 )( 10, 6 )( 8, 8 )( 16, 8 )( 16, 10 )( 16, 10 )( 16, 8 )( 16, 8 )

0.0
0.22
0.0
0.0
0.14
0.27
0.23
0.0

−
0◦
−
−
0◦
0◦
60◦
−

0.0
0.0
0.0
0.0
0.99
0.57
0.0
0.0

0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0

31.21
17.00
16.81
8.94
4.55
3.90
4.84
3.41

2.96
7.85
14.07
18.23
21.22
23.70
27.53
30.24

67

A.2

Ground state of Oxygen (Z = 8) isotopes
14 O
16 O
18 O
20 O
22 O
24 O
26 O
28 O

A.3

Orbitals (n, p)

β

γ

∆n

∆p

−λn

−λp

102.805
127.942142.569
156.168
169.172
178.613179.703180.673( 8, 14 )( 16, 14 )( 16, 14 )( 16, 14 )( 16, 14 )( 20, 14 )( 20, 14 )( 20, 14 )

0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0

−
−
−
−
−
−
−
−

0.0
0.0
1.05
1.10
0.0
0.0
0.8
0.0

0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0

20.81
13.53
7.94
7.45
8.06
5.17
1.13
1.78

5.68
10.25
13.90
17.34
20.65
22.62
24.87
27.05

Ground state of Neon (Z = 10) isotopes
16 Ne
18 Ne
20 Ne
22 Ne
24 Ne
26 Ne
28 Ne
30 Ne
32 Ne

A.4

−ETotal−ETotalOrbitals (n, p)

β

γ

∆n

∆p

−λn

−λp

102.421135.283
157.672
178.516195.051208.529213.884
219.303222.944( 10, 12 )( 14, 16 )( 16, 16 )( 20, 16 )( 20, 16 )( 20, 16 )( 20, 16 )( 20, 20 )( 24, 16 )

0.25
0.22
0.37
0.37
0.17
0.0
0.0
0.0
0.36

0◦
0◦
0◦
0◦
60◦
−
−
−
0◦

0.0
0.0
0.0
0.0
0.0
0.0
0.79
1.37
0.95

0.54
1.07
0.0
0.0
0.74
1.00
1.01
1.01
0.0

23.05
17.15
13.07
11.03
10.57
7.17
3.21
2.40
2.15

0.91
4.00
9.18
12.38
13.03
14.92
17.05
19.08
23.60

Ground state of Magnesium (Z = 12) isotopes

In 26 Mg, the ground state of HF and HF+BCS are indicated in order to express the remarkableeﬀect of pairing. Pairing changes the shape of 26 Mg from prolate to oblate.
18 Mg
20 Mg
22 Mg
24 Mg
26 MgHF
26 MgBCS
28 Mg
30 Mg
32 Mg
34 Mg
36 Mg
38 Mg
40 Mg

−ETotalOrbitals (n, p)

β

γ

∆n

∆p

−λn

−λp

101.562140.702169.987
197.299217.879217.945235.542246.237255.205
264.071270.395
275.459278.420( 10, 16 )( 14, 16 )( 16, 20 )( 20, 20 )( 14, 12 )( 22, 16 )( 28, 16 )( 28, 16 )( 28, 16 )( 28, 20 )( 28, 20 )( 30, 20 )( 30, 20 )

0.31
0.0
0.38
0.39
0.24
0.20
0.0
0.0
0.0
0.37
0.33
0.30
0.29

0◦
−
0◦
0◦
8◦
54◦
−
−
−
0◦
0◦
0◦
0◦

0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.31
0.0
1.45
1.43
1.47
0.91

0.0
1.13
0.0
0.0
0.0
0.86
1.03
1.03
1.03
0.0
0.0
0.0
0.0

25.59
20.53
16.305
14.11
11.37
13.08
9.21
5.48
5.82
4.11
3.21
2.37
1.31

0.20
2.83
6.42
9.50
11.67
11.27
13.29
15.48
17.54
20.17
21.95
23.68
35.27

68

A.5

Ground state of Silicon (Z = 14) isotopes−ETotal
24 Si
26 Si
28 Si
30 Si
32 Si
34 Si
36 Si
38 Si
38 Si
40 Si
42 Si
44 Si
46 Si

A.6

176.693208.177
237.846259.619274.482287.935
299.363309.073309.122317.662324.681330.254333.798

Orbitals (n, p)( 16, 18 )( 16, 24 )( 18, 20 )( 28, 20 )( 32, 20 )( 32, 20 )( 32, 20 )( 32, 20 )( 32, 22 )( 32, 22 )( 32, 22 )( 32, 20 )( 34, 20 )

β

γ

∆n

∆p

−λn

−λp

0.182
0.20
0.23
0.0
0.0
0.0
0.0
0.09
0.12
0.17
0.19
0.0
0.0

60◦

0.71
0.85
0.0
0.0
1.31
1.63
2.15
2.217
2.23
2.03
1.55
1.86
1.28

0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0

16.93
15.88
15.71
11.25
7.63
6.56
5.61
4.95
5.00
4.34
3.39
2.60
1.64

5.25
7.68
10.35
12.74
14.91
16.94
18.81
19.82
20.16
21.82
23.40
24.98
26.15

52◦

60◦
−
−
−
−
8◦
60◦
60◦
60◦
−
−

Ground state of Sulfur (Z = 16) isotopes

26 S
28 S
30 S
32 S
34 S
36 S
38 S
40 S
42 S
44 S
46 S
48 S
50 S

−ETotalOrbitals (n, p)

β

γ

∆n

∆p

−λn

−λp

179.410214.730
248.398
275.059294.106
311.394326.339339.607351.682361.733
369.791375.996380.092( 16, 18 )( 16, 24 )( 20, 28 )( 28, 28 )( 32, 28 )( 34, 28 )( 34, 28 )( 34, 26 )( 34, 26 )( 34, 28 )( 34, 26 )( 34, 28 )( 34, 28 )

0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.145
0.150
0.11
0.0
0.0
0.0

−
−
−
−
−
−
−
0◦
0◦
0◦
−
−
−

1.01
1.03
0.0
0.0
1.28
1.50
2.15
2.08
2.07
1.97
1.71
0.93
0.0

0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
1.70

18.62
17.84
18.06
13.16
9.61
8.35
7.34
6.73
5.82
4.87
3.97
2.77
2.21

1.63
3.49
5.39
7.30
9.25
11.08
12.76
14.94
16.65
18.00
18.97
20.24
20.98

69

A.7

Ground state of Argon (Z = 18) isotopes
34 Ar
36 Ar
38 Ar
40 Ar
42 Ar
44 Ar
46 Ar
48 Ar
50 Ar
52 Ar
54 Ar

A.8

−ETotalOrbitals (n, p)

β

γ

∆n

∆p

−λn

−λp

281.631306.247
328.262347.083364.247379.713393.338
404.114414.205
419.686422.555( 28, 20 )( 28, 28 )( 32, 28 )( 34, 28 )( 34, 28 )( 34, 28 )( 34, 28 )( 36, 28 )( 38, 28 )( 40, 28 )( 40, 28 )

0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0

−
−
−
−
−
−
−
−
−
−
−

0.0
1.22
1.11
1.99
2.16
2.01
1.78
1.41
0.83
1.20
1.41

0.80
1.26
1.24
1.22
1.20
1.18
1.17
1.15
2.04
1.96
1.85

15.11
2.03
0.48
9.27
8.43
7.61
6.49
5.31
3.95
2.36
1.57

3.25
5.77
8.02
9.93
11.79
13.57
15.19
16.55
17.61
18.75
20.12

Ground state of Calcium (Z = 20) isotopes
34 Ca
36 Ca
38 Ca
40 Ca
42 Ca
44 Ca
46 Ca
48 Ca
50 Ca
52 Ca
54 Ca
56 Ca
58 Ca
60 Ca
62 Ca
64 Ca

−ETotalOrbitals (n, p)

β

γ

∆n

∆p

−λn

−λp

251.438
286.034315.095
341.559364.455
385.027404.074420.771435.251446.220455.556462.589468.144473.074477.383
481.223( 20, 20 )( 28, 20 )( 28, 28 )( 32, 28 )( 34, 32 )( 34, 30 )( 34, 30 )( 28, 20 )( 40, 28 )( 50, 28 )( 50, 28 )( 50, 28 )( 50, 30 )( 50, 28 )( 50, 28 )( 50, 28 )

0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0

−
−
−
−
−
−
−
−
−
−
−
−
−
−
−
−

0.0
0.0
1.19
0.0
1.92
2.05
1.83
0.0
1.46
1.49
2.28
2.45
2.52
2.53
2.52
2.51

0.0
0.0
1.10
0.0
0.82
0.0
0.0
0.0
1.90
0.0
0.0
1.41
0.0
0.0
0.0
0.75

21.99
16.91
14.29
14.31
11.09
10.21
9.35
10.41
6.53
5.18
4.26
3.53
3.00
2.54
2.09
1.61

1.13
2.92
3.64
7.47
7.45
11.12
12.88
14.68
13.76
16.99
18.36
17.49
21.12
22.48
23.81
22.44

70

E1 strength distributions of Carbon (Z = 6) isotopes (N = 2 − 16)

A.9
4
3

(a) 8C(b) 10C(d) 14C(e) 16C(g) 20C(h) 22C

Total
K=0
K=1

(c) 12C

2S(E;E1) [e2fm2/MeV]

103(f) 18C

2
1
0
3

5 10 15 20 25 30 35 40

2
1
0
0

5 10 15 20 25 30 35 0

5 10 15 20 25 30 35 40

E [ MeV ]Figure A.1: E1 strength distributions of carbon isotopes. Solid line means the sum of strengths
for K = 0 and K = 1 channels. Dashed line indicates the strength of K = 0 channel and dottedline indicates the strength of K = 1 channel. The smoothing parameter of Γ = 1 MeV is used.

71

A.10

E1 strength distributions of Oxygen (Z = 8) isotopes (N = 6 − 20)

6
Total
4

(a)

14

(d)

20

(g)

26

O

(b) 16O(c) 18O

S(E;E1) [e2fm2/MeV]

2
0
4

O

(e)

22

(h)

28

O

(f)

24

O

2
0
5 10 15 20 25 30 35 40
4

O

O

2
0
0

5 10 15 20 25 30 35 0

5 10 15 20 25 30 35 40

E [ MeV ]Figure A.2: E1 strength distributions of oxygen isotopes. The meanings of the lines are sameas ﬁgure A.1. The smoothing parameter of Γ = 1 MeV is used.

72

A.11

E1 strength distributions of Neon (Z = 10) isotopes (N = 6 − 22)

6
Total
4

(a)

16

(d)

22

Ne

(e)

24

(g) 28Ne

(h)

30

Ne

K=0
K=1

(b) 18Ne(c) 20Ne

S(E;E1) [e2fm2/MeV]

2
0
4

Ne

(f)

26

Ne

2
0
4

(i) 32Ne

Ne

2
0
0

5 10 15 20 25 30 35 0

5 10 15 20 25 30 35 0

5 10 15 20 25 30 35 40

E [ MeV ]Figure A.3: E1 strength distributions of neon isotopes. The meanings of the lines are same as
ﬁgure A.1. The smoothing parameter of Γ = 1 MeV is used.

73

A.12

E1 strength distributions of Magnesium (Z = 12) isotopes (N = 6 − 28)

6
Total
4

(a)

18

(d)

24

Mg

K=0
K=1

(b) 20Mg(c) 22Mg

20

S(E;E1) [e2fm2/MeV]

4

Mg

(e)

26

Mg

(f)

28

Mg

2
0
4

(g) 30Mg(h) 32Mg(i) 34Mg(j) 36Mg(k) 38Mg(l) 40Mg

2
0
4
2
0
0

5 10 15 20 25 30 35 0

5 10 15 20 25 30 35 0

5 10 15 20 25 30 35 40

E [ MeV ]Figure A.4: E1 strength distributions of magnesium isotopes. The meanings of the lines aresame as ﬁgure A.1. The smoothing parameter of Γ = 1 MeV is used.

74

A.13

6

E1 strength distributions of Silicon (Z = 14) isotopes (N = 10 − 32)(a) 24Si

Total(b) 26Si(c) 28Si(d) 30Si(e) 32Si(f) 34Si(g) 36Si(h) 38Si(i) 40Si(j) 42Si(k) 44Si(l) 46Si

K=0
K=1

4
2
0
6

S(E;E1) [e2fm2/MeV]

4
2
0
6
4
2
0
6
4
2
0
0

5 10 15 20 25 30 35 0

5 10 15 20 25 30 35 0

5 10 15 20 25 30 35 40

E [ MeV ]Figure A.5: E1 strength distributions of silicon isotopes. The meanings of the lines are same as
ﬁgure A.1. The smoothing parameter of Γ = 1 MeV is used.

75

A.14

6

E1 strength distributions of Sulfur (Z = 16) isotopes (N = 10 − 32)(a) 26S(b) 28S(c) 30S(d) 32S(e) 34S(f) 36S(g) 38S(h) 40S(j) 44S(k) 46S

420
6S(E;E1) [e2fm2/MeV]

4
2
0
6

Total
K=0
K=1

(i) 42S4
206(l) 48S

4
20
0

5 10 15 20 25 30 35 0

5 10 15 20 25 30 35 0

5 10 15 20 25 30 35 40

E [ MeV ]Figure A.6: E1 strength distributions of sulfur isotopes. The meanings of the lines are same as
ﬁgure A.1. The smoothing parameter of Γ = 1 MeV is used.

76

A.15

E1 strength distributions of Argon (Z = 18) isotopes (N = 16 − 36)
Total

9

34

Ar

(b) 36Ar(c) 38Ar(d) 40Ar(e) 42Ar(f) 44Ar(g) 46Ar(h) 48Ar(i) 50Ar(j) 52Ar(k) 54Ar

(a)

6
3
0
9

S(E;E1) [e2fm2/MeV]

6
3
0
9
6
3
0
9

5 10 15 20 25 30 35 40

6
3
0
0

5 10 15 20 25 30 35 0

5 10 15 20 25 30 35 40

E [ MeV ]Figure A.7: E1 strength distributions of argon isotopes. The meanings of the lines are same as
ﬁgure A.1. The smoothing parameter of Γ = 1 MeV is used.

77

A.16

E1 strength distributions of Calcium (Z = 20) isotopes (N = 14 − 44)
Total

9

34

Ca

(b) 36Ca(c) 38Ca(d) 40Ca(e) 42Ca(f) 44Ca(g) 46Ca(h) 48Ca(a)

6

S(E;E1) [e2fm2/MeV]

3
0
9
6
3
0
9

5 10 15 20 25 30 35 40

6
3
0
0

5 10 15 20 25 30 35 0

5 10 15 20 25 30 35 40

E [ MeV ]Figure A.8: E1 strength distributions of calcium isotopes (N = 14 − 28). The meanings of thelines are same as ﬁgure A.1. The smoothing parameter of Γ = 1 MeV is used.

78

15

(a)

(d)

56

Total

50

Ca

(b)

52

Ca

(c)

54

(e)

58

Ca

(f)

60

Ca

10

S(E;E1) [e2fm2/MeV]

5
0

Ca

Ca

10
5
0

(g) 62Ca

5 10 15 20 25 30 35 40(h) 64Ca10
50
0

5 10 15 20 25 30 35 0

5 10 15 20 25 30 35 40

E [ MeV ]Figure A.9: E1 strength distributions of calcium isotopes (N = 30 − 44). The meanings of thelines are same as ﬁgure A.1. The smoothing parameter of Γ = 1 MeV is used.

79

B

Imaginary-Time method

In this section, we show the Gradient Method which is used in the present study for calculationof ground state. And we use also constraint imaginary-time method. We introduce the relation between gradient method and imaginary-time method, and then show the constraints inimaginary-time method

B.1

Gradient method and Imaginary-time method

We consider the gradient method for an energy functional E of a system with respect toϕ, ϕ∗ , v, v ∗ . The variation of the E is written by
δE =

δE
δE
δE
δE
δϕ + ∗ δϕ∗ +δv + ∗ δv ∗ .
δϕ
δϕ
δv
δv

(B-1)

When the parameters ϕ, ϕ∗ , v, v ∗ have some indexes(⃗r or l), δE is sum for the space expressed
each index. Next, we consider forms of δϕ, δv asδϕ ≡ −ηδEδϕ∗δv ≡ −η

,

δE
,
δv ∗

(B-2)

where η is assumed real and positive number. If we can chose δϕ∗ , δv ∗ as conjugates of δϕ, δv,
the variation δE can be written as
}
{
δE 2
δE 2
< 0.
(B-3)
δE = −2η+δϕδv
In this case, δE is always minus number. Then the energy E decreases monotonically in thisvariation, and ﬁnally the E achieves a ground state. If we can get the gradient of δϕ, δv likeabove, we obtain the state of ground state. These procedure is called gradient method, but the
points of gradient method is the search of good gradient δϕ, δv. In general, the gradient can befound through Hamiltonian of system. The time-evolution of eigenstate |ϕ⟩ of Hamiltonian h is
described as|ϕ(t + ∆t)⟩ = eih∆t |ϕ(t)⟩ ≃ |ϕ(t)⟩ + ih∆t|ϕ(t)⟩(B-4)

where ∆t is small value. If we chose the ∆t of complex value, the variation of |ϕ⟩ with respect
to ∆t is equivalent to the variation δϕ demanded on gradient method. These procedure is called
imaginary-time method.
As a concrete example, we show the case of Hartree-Fock calculation. Let us consider two
Slater determinants, |Φ⟩ composed single-particle wave functions ϕi (i = 1, · · · , A), and |Φ + δΦ⟩composed ϕi + δϕi . δϕi is the variation with respect to ∆t. The energy expectation value with
|Φ + δΦ⟩ can be expressed by
E[Φ + δΦ] = E[Φ] +

A ∫
∑
i=1

∑
δE
d⃗r
δϕi (⃗r) +
δϕi
A

i=1

81

∫
d⃗r

δE ∗
δϕ (⃗r) + O[δϕ2i , δϕ∗2i ],δϕ∗i i(B-5)

where E[Φ] ≡ ⟨Φ|H|Φ⟩. If we chose δϕi as −∆τ δE/δϕ∗i , where ∆τ is an inﬁnitesimal positiveand real value, (B-5) is rewritten as,
E[Φ + δΦ] = E[Φ] − 2∆τ

A ∫
∑
i=1

δE 2
d⃗r
δϕi

< E[Φ].

(B-6)

The (B-6) means that this gradient procedure produce the energy minimum state. And notethat the variation δE/δϕ∗i is a deﬁnition of single-particle Hamiltonian, and the coeﬃcient ∆τis equivalent to the time-step ∆t replaced complex value in HF equation.

B.2

Constraint imaginary-time method

In present study, we add some constraints to ground state in order to remove a center of mass
motion and rotations of principal axes. We show the concrete form of constraint imaginary-time
method in present study.
The constraints added by us are(i) for center of mass,A∑⟨Ψ|⃗ri |Ψ⟩ = 0, ⃗ri = ⃗xi + ⃗yi + ⃗zi ,
(B-7)
i=1

(ii) for principal axes,
⟨Ψ|

A
∑
i=1

⃗xi ⃗yi |Ψ⟩ = ⟨Ψ|A∑⃗yi⃗zi |Ψ⟩ = ⟨Ψ|

i=1

A
∑

⃗zi ⃗xi |Ψ⟩ = 0,

(B-8)

i=1

(iii) for orthogonal-normalization of single-particle states,⟨ϕk |ϕl ⟩ = δkl .

(B-9)(iv) for expectation value of particle number,⟨Ψ|N̂ |Ψ⟩ = A,

(B-10)We set the Lagrange multipliers for the constraint conditions. ξµ , ekl are correspond to eachconstraints. The variation δϕl is|ϕ′l ⟩ = |ϕl ⟩ + |δϕl ⟩


6


∑
∑
→ |ϕl ⟩ − η ĥ|ϕl ⟩ −ξµ q̂µ |ϕl ⟩ −elk |ϕk ⟩ ,µ=1(B-11)kwhere η is inﬁnitesimal positive and real value and |ϕ′l ⟩ means the state developed one imaginarytime step. q̂µ are corresponding to spatial constraints, µ is index for each centers of mass82

{⃗xi , ⃗yi , ⃗zi } and principal axes {⃗xi ⃗yi , ⃗yi⃗zi , ⃗zi ⃗xi } (µ = 1 ∼ 6). ekl is correspond to the Lagrange multiplier for orthogonal-normalization of single-particle states. ekl is derived from with⟨ϕ′k |ϕ′l ⟩ = δkl .
{⟨ϕ′k |ϕ′l ⟩= ⟨ϕk |ϕl ⟩ − η†⟨ϕk |ĥ|ϕl ⟩ + ⟨ϕk |ĥ |ϕl ⟩ −)(ξµ ⟨ϕk |q̂µ |ϕl ⟩ + ⟨ϕk |q̂µ |ϕl ⟩

6
∑
µ=1

−

∑

eml ⟨ϕk |ϕm ⟩ −l,m⇒ δkl − η⟨ϕk |gl ⟩ + ⟨gk |ϕl ⟩ − 2

6
∑

∑
k,m

ξµ ⟨ϕk |q̂µ |ϕl ⟩ − 2 ekl

µ=1



ekm ⟨ϕm |ϕl ⟩






,

(B-12)|gl ⟩ ≡ ĥ|ϕl ⟩.

When ⟨ϕ′k |ϕ′l ⟩ = δkl , ekl is
6
) ∑
1(
ekl =
⟨ϕk |gl ⟩ + ⟨gk |ϕl ⟩ −ξµ ⟨ϕk |q̂µ |ϕl ⟩.
2(B-13)µ=1

But ekl needs the value of ξµ . For the spatial constraints q̂µ , the expectation value of ξµ with
BCS wave function is written by∫⟨Ψ|q̂µ |Ψ⟩ =d⃗r qµ (⃗r)ρ(⃗r))∑ (vl2 ⟨ϕl |q̂µ |ϕl ⟩ + ⟨ϕl̄ |q̂µ |ϕl̄ ⟩ = 0,
(B-14)=
l>0

here q̂µ is one-body operator, the expectation value can be expressed by the density of a system.
We can consider the value of ξµ with aboves then.
∑vl2 ⟨ϕ′l |q̂µ |ϕ′l ⟩ = 0
l>0

=

∑

{
vl2

(

⟨ϕl |q̂µ − η ⟨ϕl |ĥ† q̂µ −

l>0

6
∑

ξν ⟨ϕl |q̂ν q̂µ −

ν=1

{

⇒

{
vl2

)}
ekl ⟨ϕk |q̂µk(

× |ϕl ⟩ − η ĥ|ϕl ⟩ −
∑

∑

6
∑

ξτ q̂τ |ϕl ⟩ −)}elm |ϕm ⟩

m

τ =1

(

∑

⟨ϕl |q̂µ |ϕl ⟩ − η ⟨gl |q̂µ |ϕl ⟩ + ⟨ϕl |q̂µ |gl ⟩ − 26∑ξν ⟨ϕl |q̂ν q̂µ |ϕl ⟩

ν=1

l>0

−

∑
k

83

ekl ⟨ϕk |q̂µ |ϕl ⟩ −∑kelk ⟨ϕl |q̂µ |ϕk ⟩)}.

(B-15)

And substitute ekl for (B-15),∑
vl2 ⟨ϕ′l |q̂µ |ϕ′l ⟩ = 0
l>0
[
{
6
∑
∑
2
⇒
vl ⟨ϕl |q̂µ |ϕl ⟩ − η ⟨gl |q̂µ |ϕl ⟩ + ⟨ϕl |q̂µ |gl ⟩ − 2ξν ⟨ϕl |q̂ν q̂µ |ϕl ⟩
l>0

ν=1



6

) ∑



−
⟨ϕk |gl ⟩ + ⟨gk |ϕl ⟩ −ξν ⟨ϕk |q̂ν |ϕl ⟩ ⟨ϕk |q̂µ |ϕl ⟩ − h.c.


2

µ=1
k
∑





1(

(B-16)S

For this S, we can use the relation A + A∗ = 2Re[A] in general.
[
]
6
(
)
∑
∑∑
Re ⟨ϕk |gl ⟩ + ⟨gk |ϕl ⟩ ⟨ϕk |q̂µ |ϕl ⟩ − 2S=−ξν ⟨ϕk |q̂ν |ϕl ⟩⟨ϕk |q̂µ |ϕl ⟩ ,
k ν=1

k

The coeﬃcient of ξν can be rewritten to simple one with Aνµ as,
2η

6
∑
ν=1

ξν Aνµ ≡ 2η

6
∑
ν=1

ξν

∑

{
[
]}
vl2 ⟨ϕl |q̂ν q̂µ |ϕl ⟩ − Re ⟨ϕk |q̂ν |ϕl ⟩⟨ϕk |q̂µ |ϕl ⟩ .

(B-17)

k,l>0

The other parameter can be replace to the parameter ζµ with one index µ as
[
{
}]
)
]
∑
∑ [(
2
ζ̃µ ≡
vl ⟨ϕl |q̂µ |ϕl ⟩ − η 2Re[⟨ϕl |q̂µ |gl ⟩] −Re ⟨ϕk |gl ⟩ + ⟨gk |ϕl ⟩ ⟨ϕk |q̂µ |ϕl ⟩(B-18)
l>0

k>0

Finally, the Lagrange multipliers of spatial constraints ξµ , ekl can be obtained from
ζµ =

6
∑

ξν Aνµ ,(B-19)

ν=1

where ζµ ≡ −ζ̃µ /2η. In addition, we can solve these coeﬃcients to need occupation probabilities
vl2 which does not need to be deﬁned at same time with single-particle state. But they of course,
need the information of basis (single-particle state) which is a single-particle energy.

84

The constraints of particle number is due to the deﬁnition of occupation probabilities v 2 , achemical-potential λ and the gap energy ∆. In order to deﬁne these parameters, we solve the gap
equation and the equation of number which they are written by εl , λ and ∆. The parameters are
λ and ∆ for in these two equations. In present study, we solve simultaneous equations following
F(λ, ∆) and G(λ, ∆).
{}∑εk − λF(λ, ∆) ≡1− √− Nτ .
(B-20)(εk − λ)2 + ∆2
k>0G(λ, ∆) ≡

Gτ ∑
∆
√
− ∆.
2
(εk − λ)2 + ∆2

(B-21)

k>0

where index τ shows isospin. The simultaneous equations for F(λ, ∆) = 0, G(λ, ∆) = 0 canbe solved with two-dimensional Newton-Raphson method. In the calculation, we need following
derivative values with respect to λ, ∆.
{
}
}
∑
∑{
1
∂F
∆(εk − λ)∂F(εk − λ)2
√
,
,
=
=
−
∂λ
∂∆
{(εk − λ)2 + ∆2 }3/2
(εk − λ)2 + ∆2 {(εk − λ)2 + ∆2 }3/2
k>0
k>0
∂G
∂λ

=

∂G
∂∆

=

Gτ ∑
∆(εk − λ),2{(εk − λ)2 + ∆2 }3/2
k>0
{
}
1
Gτ ∑
∆2
√
−
− 1.
2
(εk − λ)2 + ∆2 {(εk − λ)2 + ∆2 }3/2
k>0

The algorithm of two-dimensional Newton-Raphson method for f (x, y) = 0, g(x, y) = 0 is
(
) (
) (
) (
)
f (xA , yA )f (xi , yi )∂x f (xi , yi ) ∂y f (xi , yi )
dx
=
+
=0
g(xA , yA )g(xi , yi )∂x g(xi , yi ) ∂y g(xi , yi )dyAwhere xA , yA are solutions and xi+1 = xi + dx, i shows iteration number. (dx, dy) can be
expressed with inversion of A as,
(
)
(
)
1
dx
−∂x f (xi , yi ) ∂x g(xi , yi )=−.
(B-22)dy
detA ∂y f (xi , yi ) − ∂y g(xi , yi )When (dx, dy) ∼ (0, 0), we can get the solutions.

85

C

BCS equation

This section shows the BCS equations which can describe Gap energy and express occupation
probabilities of orbitals in the system interacted pairing correlation. First, we set Hamiltonian
H written by the creation and annihilation operators of Fermion as∑
1 ∑tαβ c†α cβ +V̄αβγδ c†α c†β cδ cγ ,
(C-1)
H=
4
α,β

α,β,γ,δ

V̄αβγδ ≡ Vαβγδ − Vαβδγ .

(C-2)

We set the BCS trial wave function |ΦBCS ⟩ as like,∏|ΦBCS ⟩ =(ul + vl c†l c†l̄ )|0⟩,

(C-3)

l>0

where ¯l means pair of l-state which are usually assumed the time-reversal relation between thepair and vl can be chosen as real number in static system. The BCS wave function does not aneigenstate of particle number, breaks the number conservation. So, usually we need the number
constraint to Hamiltonian in the BCS theory. H is modiﬁed as,
H ′ ≡ H − λN̂ ,∑⟨ΦBCS |N̂ |ΦBCS ⟩ = 2
vl2 ≡ A,

(C-4)

l>0

where N̂ is number operator, λ is a Lagrange multiplier and A is particle number of the system.
H ′ is a standard Hamiltonian in BCS theory. The expectation value of this H ′ with |ΦBCS ⟩ is,{)∑ (⟨ΦBCS |H ′ |ΦBCS ⟩ =tkk + tk̄k̄ − 2λ vk2
k>0

}
∑
) 2 2
1 ∑(
+
V̄kk′ kk′ + V̄k̄k¯′ k̄k¯′ + V̄kk¯′ kk¯′ + V̄k̄k′ k̄k′ vk vk′ +V̄kk̄k′ k¯′ uk vk uk′ vk′
2 ′
k >0
k′ >0
}
{
)
∑
∑
∑ (
1
V̄kk̄k′ k¯′ uk vk uk′ vk′V̄kk′ kk′ vk2 vk2′ +(C-5)=
tkk − λ vk2 +
2 ′
′
k

k

k,k >0

The BCS equation can be derived from the variation δ⟨ΦBCS |H ′ |ΦBCS ⟩ = 0, as like
(
)
∂
∂ul ∂
′
δ⟨ΦBCS |H |ΦBCS ⟩ =+⟨ΦBCS |H ′ |ΦBCS ⟩∂vl∂vl ∂ul
{
}
(
) ∑(
)
2
= 2vl
tll + tl̄l̄ − 2λ +V̄lklk + V̄l̄kl̄k vk
k+ 2ul⇒ ul vl{(∑V̄ll̄kk̄ uk vk − 2

k>0

)

)
∑(
tll + tl̄l̄ − 2λ +V̄lklk + V̄l̄kl̄k vk2
k

}

vl ∑
vl
V̄ll̄kk̄ uk vk = 0.
ul
k>0

+ (u2l − vl2 )

∑
k>0

87

V̄ll̄kk̄ uk vk = 0. (C-6)

We can rewrite (C-6) to simple one with the following values.
{
}
(
) ∑(
)
)
∑(
1
2
2
tll + tl̄l̄ − 2λ +V̄lklk + V̄l̄kl̄k vk =tll + tl̄l̄ +V̄lklk + V̄l̄kl̄k vk′ − λ ≡ ε̃l (C-7)
2
k

−

k

∑

V̄ll̄k′ k¯′ uk′ vk′ ≡ ∆l

(C-8)

k′ >0

Finally, the BCS equation is described with a normalization condition u2l + vl2 = 1, as
2ul vl ε̃l + ∆l (2vl2 − 1) = 0.

(C-9)

But, we note that this BCS equation is assumed that there is the time-reversal relation between
l- and ¯l-state. If we do not assume the time-reversal relation between the pair, BCS equation is
modiﬁed as
ul vl (ε̄l + ε̄l̄ ) + ∆l (2vl2 − 1) = 0,

(C-10)∑
where ε̄l ≡ tll + k Vlklk vk2 − λ.
vl2 , u2l can be formed with the normalization condition as,
{
}
{
}
1
ε̃
1
ε̃
l
l
vl2 =
1± √
,
u2l =
1± √
.
2
2
2
2
2
ε̃ + ∆ε̃ + ∆2lll(C-11)l

In the case of no pairing correlation, one has δ = 0 and vl2 = 1, u2l = 0 for occupied orbital(ε̃l < 0). The only possible the solutions of (C-11) are therefore:
}
}
{
{
1
1
ε̃
ε̃
l
l
, u2l =
(C-12)
vl2 =
1− √
1+ √
2
2
2
2
2
ε̃ + ∆ε̃ + ∆2llll

We can obtain the Gap equation with (C-8) and (C-12),∆l = −

1∑∆kV̄ll̄kk̄ √.
2ε̃2 + ∆2k>0kk

In general, vl2 , u2l , ∆l are nonlinear equations and have to be solved by iteration.

88(C-13)

D

Intrinsic frame and Deformation parameters

In this section, we show the deformation parameters β, γ which characterize quadrupole deformation, and the expression for β- and γ-vibration.

D.1

Intrinsic frame

We ﬁrst deﬁne the general deformation parameter with a function described by the length of
the radius vector pointing from the origin to the surface. The function is


∑ µ=λ
∑
R(θ, φ) = R0 1 +a∗λµ Yλ,µ (θ, φ) .
(D-1)
λ=0 µ=−λR(θ, φ) shows the length of the radius vector pointing from the origin to the surface. R0 is
the radius of the sphere with a same nuclear volume. The residual part means the diﬀerence
from sphere shape with the spherical harmonic function Yλ,µ (θ, φ). The complex parameter a∗λµshows the weight of Yλ,µ in R(θ, φ), so then we can consider the various shape of nucleus withthis variation of a∗λµ . These θ, φ has the relation with Cartesian coordinate, as likex = r sinθ cosφ,y = r sinθ sinφ,z = r cosθ.

We redeﬁne θ, φ from z-axis in the Cartesian coordinate whose origin agree with the center ofmass of nucleus as like Fig. D.10. The coordinate system redeﬁned from the center of mass iscalled intrinsic frame or body-ﬁxed frame. In adding, a∗λµ has a following symmetry,a∗λµ = (−1)µ aλ−µ .

(D-2)∗ (θ, φ) =This is because that R(θ, φ) is real and the spherical harmonic function Yλ,µ (θ, φ) has Yλ,µµ(−1) Yλ,−µ (θ, φ). The R can express various shape, but a00 , a1µ are excluded from nuclear deformation parameter because nuclear system has the saturation property and the center of massof nucleus is ﬁxed. We should consider the nuclear deformation from quadrupole case (λ ≥ 2).

D.2

Quadrupole deformation parameter (β, γ)We consider only the nuclear quadrupole deformation (ellipsoid shape). Now we chose theintrinsic frame as like Fig. D.10. We call this frame K ′ -frame. In K ′ -frame, angle parameters
is written as θ′ , φ′ and nuclear surface R(θ′ , φ′ ) is also written by


µ=2
∑
R(θ′ , φ′ ) = R0 1 +a∗2µ Y2,µ (θ′ , φ′ ) .
(D-3)
µ=−2

This R(θ′ , φ′ ) is assumed ellipsoid shape, then R(θ′ , φ′ ) is reﬂection symmetric about x′ y ′ -, y ′ z ′ and z ′ x′ -plane. In this case, the deformation parameter a∗2µ are restricted and can be reduced.
Now we chose the z ′ -axis to axial symmetric one. There are some restrictions for a∗2µ , from thereﬂection symmetry about x′ y ′ -plane, a∗21 = a∗2−1 = 0 and from the symmetry about y ′ z ′ - and90

z ′ x′ -plane, a∗22 = a∗2−2 . The ﬁve coeﬃcients a2µ reduce to two real independent variables a20
and a22 = a2−2 which together with the three Euler angles, give a complete description of thissystem. Usually, we introduce instead of a20 and a22 the so-called Hill-Wheeler coordinates β, γthrough the relation{
a20 = β cosγ,(D-4)
a22 = a2−2 = √12 β sinγ,
R can be expressed in the directions of K ′ -frame as(π )(π π )
R, 0 ≡ R1 , R,≡ R2 , R (0, 0) ≡ R3 ,
2
2 2
{
√
(
)}
5
2π
Rk = R0 1 +β cos γ − k
.
4π
3

(D-5)

k = 1, 2, 3 ≡ {x′ , y ′ , z ′ }The variation δRk (≡ Rk − R0 ) can indicate several shapes with Eq.(D-6). In the case of β = 0,
δRk = 0 namely the shape of nucleus is spherical. In the case of β > 0 with γ = 0, δR1 =δR2 < 0, R3 > R0 namely the shape is prolate which is shown in the left side of Fig. D.11. In
the case of β > 0 with γ = π, δR1 = δR2 > 0, R3 < R0 these values indicate the oblate shapewhich is shown in the right side of Fig. D.11, but there is a same oblate shape in the case of

Figure D.11: Nuclear shapes
Figure D.10: K ′ -frame
β < 0 with γ = 0 namely the case of δR1 = δR2 > 0, R3 < R0 . In general, we can express
the axial symmetric deformation in a restricted (β, γ) region. We commonly chose the regionof β ≥ 0 , 0 ≤ γ ≤ π/3 as (β, γ)-plane in order to express the axial symmetric deformation. In
adding, the case of β > 0, γ > 0 indicate R1 ̸= R2 ̸= R3 , this means this case does not belong
to both prolate and oblate case, so we call this case triaxial deformation.

91

D.3

β- and γ-vibration

We can consider the vibration of nuclear surface with deformation parameters {a20 , a22 } or{β, γ}. This means deformation parameters can describe some nuclear state called collective
motion. If we can consider the collective motion in their space, they are often called “collective
parameter” also. When a nucleus has deformation with (β0 , γ0 ) in the ground state, the nucleus
is an equilibrium point on the energy surface expressed (β, γ)-plane. This means that we canconsider the nuclear potential surface with respect to (β, γ). So, the vibration of (β, γ) values
can be regarded the vibration around the equilibrium point on (β, γ)-plane as like Fig. D.12.

Figure D.12: Aspects of V (β, γ). From left side, oblate, spherical and prolate cases are shown.
Now we consider some inﬁnitesimal vibration (δβ, δγ) around (β0 , γ0 ). in both expressions {β, γ}and {a20 , a22 }. The inﬁnitesimal vibration can be expressed as(β0 , γ0 )→ ( β0 + δβ, γ0 + δγ ),
(
)
(
)
1
(0) (0)a20 , a22 → (β0 + δβ)cos(γ0 + δγ), √ (β0 + δβ)cos(γ0 + δγ)2()
1≈ β0 + δβ , √ β0 δγ = (a20 , a22 ),δβ · δγ ≈ 0,
2
(0)

(0)

where (a20 , a22 ) is the expression of (β0 , γ0 ) in {a20 , a22 } and using Eq.(D-4). And we can see
the R(θ′ , φ′ ) with these values as{()}
R(θ′ , φ′ ) = R0 1 + a∗20 Y2,0 (θ′ , φ′ ) + a∗22 Y2,2 (θ′ , φ′ ) + Y2,−2 (θ′ , φ′ ).
(D-6)
We can understand the corresponds between a20 ⇔ Y2,0 and a22 ⇔ Y2,2 . The operators,)
1 (Q20 ≡ r2 Y2,0 and Q22 ≡ r2 √ Y2,2 + Y2,−2 ,(D-7)
2produce the vibration of δβ and δγ. They are called β-vibration and γ-vibration.

D.4

Definition of quadrupole deformation parameter (β, γ)In this thesis, we use following quadrupole deformation parameters (β, γ),√()
4π ⟨Q20 ⟩2 + ⟨Q22 ⟩2
180◦⟨Q22 ⟩β≡, γ≡arctan,(D-8)
5A⟨r2 ⟩π⟨Q20 ⟩√
where ⟨r2 ⟩ is root mean square radius of nucleus. These parameters are given also in [79, 80].
92

Bibliography[1] A.Bohr and B.R.Mottelson, “N uclear Structure, Vol.I”(Benjamin, New York,1969)[2] A.Bohr and B.R.Mottelson, “N uclear Structure, Vol.II”(Benjamin, New York,1975)[3] P. Ring and P. Schuck, “T he N uclear M any−Body P roblem” (Springer-Verlag, 1980).
[4] J.-P. Blaizot and G. Ripka, Quantum Theory of Finite Systems(MIT Press, Cambridge,
MA, 1986).
[5] S.C.Pieper , Nucl. Phys. A751(2005), 516.
[6] M.G.Mayer, Phys. Rev. 75(1949)1969.
[7] C.H.Townes, H.M.Foley and W.Low, Phys. Rev. 76(1949)1415.
[8] S.G.Nilsson, Dan. Mat . Fys. Medd. 29, no.16 (1955).
[9] W.Nazarawicz, J.Dudek, R.Bengtsson, T.Bengtsson and I.Ragnarsson, Nucl. Phys.
A435(1985), 397.
[10] A.Bohr, Dan. Mat . Fys. Medd. 26, no.14 (1952).
[11] A.Bohr and B.R.Mottelson, Dan. Mat . Fys. Medd. 27, no.16 (1953).
[12] J.W.Negele, Rev. Mod. Phys. 54, 913 (1982).
[13] P.Bonche, S.Koonin and J.W.Negele, Phys. Rev. C13, 1226 (1976).
[14] R.Y.Cusson, R.K.Smith and J.A.Maruhn, Phys. Rev. Lett. 36, 1166 (1976).
[15] H.Flocard, S.E.Koonin and M.S.Weiss, Phys. Rev. C17, 1682 (1978).
[16] R.Y.Cusson, J.A.Maruhn and H.W.Meldner, Phys. Rev. C18, 2589 (1978).
[17] J.W.Negele, S.E.Koonin, P.Möller, J.R.Nix, and A.J.Sierk, Phys. Rev. C17, 1098 (1978).
[18] J.A.Maruhn, P.G.Reinhard, P.D.Stevenson, J. R.Stone, and M.R.Strayer, Phys. Rev. C71,
064328 (2005).
[19] A.S.Umar and V.E.Oberacker, Phys. Rev. C73, 054607 (2006).
95

[20] A.S.Umar and V.E.Oberacker, Phys. Rev. C76, 014614 (2007).
[21] K.Washiyama and D.Lacroix, Phys. Rev. C78, 024610 (2008).
[22] C.Golabek and C.Simenel, Phys. Rev. Lett. 103, 042701 (2009).
[23] T.Nakatsukasa and K.Yabana, J. Chem. Phys. 114, 2550 (2001)[24] T.Nakatsukasa and K.Yabana, Phys. Rev. C71, 024301 (2005).
[25] H.A.Bethe and R.F.Bacher, Rev. Mod. Phys. 8(1936)82.
[26] A.Bohr, B.R.Mottelson and D.Pines, Phys. Rev. 110, 936 (1958).
[27] J.Bardeen, L.N.Cooper and J.R.Shrieﬀer, Phys. Rev. 108, 1175 (1957).
[28] W.J.Swiatecki, Phys. Rev. 100(1955)937.
[29] Z.Ren and C.Xu, Nucl. Phys. A759(2005), 64.
[30] J.H.E.Mattuch, W.Thiele and A.H.Wapstra, Nucl. Phys. 67(1965), 1.
[31] J.L.Wood, K.Heyde, W.Nazarawicz, M.Huyse and P.Van Duppen, Phys. Rep. 215 (1992)
101.
[32] S.G.Nilsson and O.Prior, Dan. Mat . Fys. Medd. 32, no.16 (1961).
[33] I.Tanihata, H.Hamagaki, O.Hashimoto, Y.Shida, N.Yoshikawa, K.Sugimoto,O.Yamakawa, T.Kobayashi and N.Takahashi, Phys. Rev. Lett. 55, 2676 (1985).
[34] I.Tanihata, J. Phys. G 22, 157 (1996).
[35] A.Ozawa, T.Kobayashi, T.Suzuki, K.Yoshida and I.Yanihata, Phys. Rev. Lett. 84, 5493
(2000).
[36] T.Otsuka, R.Fujimoto, Y.Utsuno, B.A.Brown,M.Honma and T.Mizusaki, Phys. Rev. Lett.
87, 082502 (2001).
[37] M.Matsuo, K.Mizuyama and Y.Serizawa, Phys. Rev. C71, 046326 (2005).
[38] M.Matsuo, Phys. Rev. C73, 044309 (2006).
[39] E.Litvinova, P.Ring and V.Tselyaev, Phys. Rev. C78, 014312 (2008).
[40] N.Paar, T.Nikšić, T.Marketin, D.Vretenar and P.Ring, Eur. Phys. J. A 25, 531 (2006).
[41] D.Vretenar, N.Paar a, P.Ringa, G.A.Lalazissis, Nucl. Phys. A692(2001), 496.
[42] S.Goriely, Phys. Lett. B 436, 10 (1998).
[43] S.Goriely and E.Khan, Nucl. Phys. A706(2002), 217.
96

[44] E.K.Warburton, J.A.Becker and B.A.Brown, Phys. Rev. C41, 1147 (1990).
[45] O.Sorlin, C.Donzaud, F.Nowacki, J.C.Angélique, F.Azaiez, C.Bourgeois, V.Chiste,
Z.Dlouhy, S.Grévy, D.Guillemaud-Mueller, F.Ibrahim, K.-L Kratz, M.Lewitowicz,
S.M.Lukyanov, J.Mrasek, Yu.-E.Penionzhkevich, F.de Oliveira Santos, B.Pfeiﬀer,
F.Pougheon, A.Poves, M.G.Saint-Laurent, M.Stanoiu, Eur. Phys. J. A 16, 55 (2003).
[46] N.Aoi, E.Takeshita, H.Suzuki, S.Takeuchi, S.Ota, H.Baba, S.Bishop, T.Fukui,
Y.Hashimoto, H.J.Ong, E.Ideguchi, K.Ieki, N.Imai, M.Ishihara, H.Iwasaki, S.Kanno,
Y.Kondo, T.Kubo, K.Kurita, K.Kusaka, T.Minemura, T.Motobayashi, T.Nakabayashi,
T.Nakamura, T.Nakao, M.Niikura, T.Okumura, T.K.Ohnishi, H.Sakurai, S.Shimoura,
R.Sugo, D.Suzuki, M.K.Suzuki, M.Tamaki, K.Tanaka, Y.Togano, and K.Yamada, Phys.
Rev. Lett. 102, 012502 (2009).
[47] O.B.Tarasov, D.J.Morrissey, A.M.Amthor, T.Baumann, D.Bazin, A.Gade, T.N.Ginter,
M.Hausmann, N.Inabe, T.Kubo, A.Nettleton, J.Pereira, M.Portillo, B.M.Sherrill, A.Stolz,
and M.Thoennessen, Phys. Rev. Lett. 102, 142501 (2009).
[48] M.Stoitsov, J.Dobaczewski, W.Nazarewicz, and P.Ring, Comput. Phys. Commun. 167, 43(2005).
[49] B.Avez, C.Simenel and Ph.Chomaz, Phys. Rev. C78, 044318 (2008).
[50] Y.Hashimoto and K.Nodeki, arXiv:0707.3083[51] J.Terasaki and J.Engel, Phys. Rev. C74, 044301 (2006).
[52] M.Matsuo, Nucl. Phys. A696(2001), 371.
[53] E.Khan, N.Sandulescu, M.Grasso, and Nguyen Van Giai, Phys. Rev. C66, 024309 (2006).
[54] M.Bender, J.Dobaczewski, J.Engel, and W.Nazarewicz, Phys. Rev. C65, 054322 (2006).
[55] K.Yoshida and Nguyen Van Giai, Phys. Rev. C78, 064316 (2008).
[56] C.Losa, A.Pastore, T.Døssing, E.Vigezzi, and R.A.Broglia, Phys. Rev. C81, 064307(2010).
[57] J.Terasaki and J.Engel, Phys. Rev. C82, 034326 (2010).
[58] J.Blocki and H.Flocard, Nucl. Phys. A273(1976), 45.
[59] T.Nakatsukasa, T.Inakura, and K.Yabana, Phys. Rev. C76, 024318 (2007).
[60] N.Tajima, S.Takahara, and N.Onishi, Nucl. Phys. A603, 23 (1996).
[61] M.Brack, J.Damgaard, A.S.Jensen, H.C.Pauli, V.M.Strutinsky, and C.Y.Wong, Rev. Mod.
Phys. 44, 320 (1972).
[62] K.T.R.Davies, H.Flocard, S.Krieger, and M.S.Weiss, Nucl. Phys. A342, 111 (1980).
97[63] J.Bartel, P.Quentin, M.Brack, C.Guet, and H.Håkansson, Nucl. Phys. A386, 79 (1982).
[64] P.Bonche, H.Flocard, and P.H.Heenen, Nucl. Phys. A467,115 (1987).
[65] E. Chabanat, P.Bonche, P.Haensel, J.Mayer and R.Schaeﬀer, Nucl. Phys. A627,710(1997).
[66] T.Inakura, T.Nakatsukasa and K.Yabana, Phys. Rev. C80, 044301 (2009).
[67] A.Veyssière, H.Beil, R.Bergère, P.Carlos, A.Leprêtre and K.Kernbath, Nucl. Phys. A199,
45 (1973).
[68] V.V.Varlamov, M.E.Stepanov, and V.V.Chesnokov, Bull. Rus. Acad. Sci. Phys. 67, 724(2003).
[69] J.Ahrens, H.Borchert, K.H.Czock, H.B.Eppler, H.Gimm, H.Gundrum, M.Kröing, P.Riehn,
G.Sita Ram, A.Zieger and B.Ziegler, Nucl. Phys. A251, 479 (1975).
[70] A.Veyssière, H.Beil, R.Berge, P.Carlos, A.Leprêtre and A.De Miniac, Nucl. Phys. A227,
513 (1974).
[71] H.B.A.Leprêtre, R.Bergère, P.Carlos, A.Veyssière and M.Sugawara, Nucl. Phys. A175,
609 (1971).
[72] A.Veyssière, H.Beil, R.Berge, P.Carlos and A.Leprêtre, Nucl. Phys. A159, 561 (1970).
[73] P.Carlos, H.Beil, R.Berge, A.Leprêtre, A.De Miniac and A.Veyssière, Nucl. Phys. A225,
171 (1974).
[74] P.Carlos, H.Beil, R.Berge, A.Leprêtre, A.De Miniac and A.Veyssière, Nucl. Phys. A172,
437 (1971).
[75] C.F.von Weizsäcker, Z. Phys. 96(1935)431.
[76] J.Gibelin, D.Beaumel, T.Motobayashi, Y.Blumenfeld, N.Aoi, H.Baba, Z.Elekes, S.Fortier,
N.Frascaria, N.Fukuda, T.Gomi, K.Ishikawa, Y.Kondo, T.Kubo, V.Lima, T.Nakamura,
A.Saito, Y.Satou, J.-A.Scarpaci, E.Takeshita, S.Takeuchi, T.Teranishi, Y.Togano,
A.M.Vinodkumar, Y.Yanagisawa, and K.Yoshida, Phys. Rev. Lett. 101, 212503 (2008).
[77] L.-G.Cao and Z.-Y.Ma, Phys. Rev. C71, 034305 (2005).
[78] J.Piekarewicz, Phys. Rev. C73, 044325 (2006).
[79] M.Yamagami and K.Matsuyanagi, Nucl. Phys. A672, 123 (2000).
[80] K.Hagino, N.W.Lwin and M.Yamagami, Phys. Rev. C74, 017310 (2006).

98