Theories and calculations
Nuclear energy density functional theory

According to the basic theorems of the density functional theory (DFT) [1], all
the physical quantities at the ground state of many-body systems can be given byfunctionals of one-body density 𝜌(𝒓). In nuclei, there are two kinds of particles,
protons and neutrons, thus, they are functionals of proton and neutron densities,
𝜌% (𝒓) and 𝜌& (𝒓). We start from a given energy density functional 𝐸(𝜌% , 𝜌& *, andtry to find the density distributions, 𝜌% (𝒓) and 𝜌& (𝒓), which minimize the energy.
In a variational form, it is given by𝛿 ,𝐸(𝜌% , 𝜌& * − 𝜇% / 𝜌% (𝑟)𝑑𝑟 − 𝜇& / 𝜌& (𝑟)𝑑𝑟2 = 0 ,(1)where 𝜇% and 𝜇& are the Lagrange multipliers to fix the numbers of protons andneutrons. Equation (𝟏) is enough in principle. However, in practice, to achieve
high accuracy of the calculation, we use the Kohn-Sham scheme [2] with a
Bogoliubov extension [3] and replace 𝐸(𝜌% , 𝜌& * by 𝐸(𝜌% , 𝜌& , 𝜏% , 𝜏& , ⃡𝐽% , ⃡𝐽& ; 𝜅% , 𝜅& *where 𝜏= (𝒓) are the kinetic density and ⃡𝐽= are the spin-current density, 𝜅= are
the pair density (𝑞 = 𝑛, 𝑝). The definitions of these densities can be found inliterature, e.g., [4]. Finally, this leads to the Kohn-Sham-Bogoliubov-de-Gennes
equations:ℎ= − 𝜇=A−Δ∗==𝑈I=(=) 𝑈I∗ G A = G = 𝐸I A = G ,−Eℎ= − 𝜇= F
𝑉I
𝑉I

Δ=(2)

where 𝑞 = 𝑛, 𝑝. Here, 𝜇= are called “chemical potentials” to control the average
numbers of protons and neutrons, ℎ= (𝒓, 𝜎) are the single-particle Hamiltonians,
Δ= (𝒓, 𝜎) are the pair fields which introduce the Cooper pairs and the superfluidity
in nuclei. The pair fields are coordinate-dependent in general, and their averagevalues are often called “pairing gap energy” denoted as Δ= . Equation (2) hassolutions E𝑈I= (𝒓, 𝜎), 𝑉I= (𝒓, 𝜎)F with positive quasi-particle energies 𝐸I= > 0 andnegative ones 𝐸&= < 0. In order to construct ℎ= (𝒓, 𝜎) and Δ= (𝒓, 𝜎) in the left
hand side of Eq. (2), we need all the solutions E𝑈I= (𝒓, 𝜎), 𝑉I= (𝒓, 𝜎)F with either
positive or negative 𝐸I= . Thus, Eq. (2) must be self-consistently solved and theiterative procedure is necessary. The theory is called “Hartree-Fock-Bogoliubov”(HFB) theory in nuclear physics. See Ref. [4,5] for more details.

Time-dependent density functional theory and canonical-basis formulation

The basic theorems for the ground state can be generalized to systems under atime-dependent external field [6, 5]. This leads to time-dependent equations:=ℎ= − 𝜇=𝜕 𝑈I𝑖ℏ A = G = A
𝜕𝑡 𝑉I−Δ∗=

Δ=−Eℎ= − 𝜇= F

𝑈I=(3)∗G A = G ,𝑉I

where all the fields and wave functions now depend on time. The time-dependent
external fields are included either in ℎ= (𝒓, 𝜎; 𝑡) or Δ= (𝒓, 𝜎; 𝑡). The self-consistentsolutions of Eq. (3) describe the time evolution of the densities of a nucleus. The
chemical potential 𝜇= is arbitrary here. In fact, the dynamics does not depend on
the value of 𝜇= . The theory is called “Time-dependent Hartree-Fock-Bogoliubov”(TDHFB) theory in nuclear physics.
Solution of Eq. (3) is still computationally difficult, since we have to calculate thetime evolution of all the quasi-particle wave functions. In order to reduce thecomputational cost, we have developed an approximate method, called“Canonical-basis

time-dependent

Hartree-Fock-Bogoliubov”(Cb-TDHFB)

method [7]. Adopting a diagonal approximation for the pair fields, we end up the=equations for a pair of canonical states (𝜙U= (𝑡), 𝜙UV (𝑡)) with 𝑞 = 𝑛, 𝑝:𝑖ℏ𝜕 =W𝜙U (𝑡)X = (ℎ= (𝑡) − 𝜂U= (𝑡)*W𝜙U= (𝑡)X,𝜕𝑡(4)𝑖ℏ𝜕 ===W𝜙UV (𝑡)X = (ℎ= (𝑡) − 𝜂UV (𝑡)*W𝜙UV (𝑡)X,𝜕𝑡(5)and the equations for their occupation and pair probabilities (𝜌U= (𝑡), 𝜅U= (𝑡)):𝑑 =\=∗
𝜌U (𝑡) = 𝜅U= (𝑡)Δ=∗U (𝑡) − 𝜅U (𝑡)ΔU (𝑡),𝑑𝑡

(6)

𝑑 =\=
𝜅U (𝑡) = 𝜅U= (𝑡)Δ=∗U (𝑡) − ΔU (𝑡)(2𝜌U (𝑡) − 1*.
𝑑𝑡(7)𝑖ℏ𝑖ℏ=Here, the real parameters (𝜂U= (𝑡), 𝜂UV (𝑡)) are chosen as the single-particle energies,𝜂U= (𝑡) = ⟨𝜙U= (𝑡)Wℎ= (𝑡)W𝜙U= (𝑡)⟩. The gap parameters Δ=U (𝑡) are the diagonal part of===the pair field, ΔU (𝑡) = −⟨𝜙U (𝑡)WΔ= (𝑡)W𝜙UV (𝑡)⟩. This diagonal approximation for thepair field leads to the “Hartree-Fock+BCS” (HF+BCS) states for the calculationof ground states. The computational cost of Cb-TDHFB is significantly smaller
than the full TDHFB of Eq. (3), roughly speaking, by a few to several orders ofmagnitude. See Ref. [7,8] for details.

Linear response and finite amplitude method

When the external fields are weak, we may linearize the TDDFT equation (3)and obtain the linear response equation. We have developed a method to facilitate
the solutions of the linear response equation, called “finite amplitude method”(FAM) [9,10,11].
Applying a weak time-dependent external field of a given frequency 𝜔, 𝐹(𝑡) =𝐹(𝜔)𝑒 efgh + h. c., to the ground state,
1 no1 onrr (𝜔)𝑎 q(𝜔)𝑎Iq 𝑎mq + 𝐹Im(𝜔)𝑎m 𝑎I + 𝐹Im
𝐹(𝜔) = l , 𝐹Im
I 𝑎m 2 ,22(8)I,mwhere (𝑎Iq , 𝑎m ) are the quasi-particle creation/annihilation operators at theground state, the linear response equations are given asno (𝜔)no (𝜔),E𝐸I + 𝐸m − ℏ𝜔F𝑋Im (𝜔) + 𝛿𝐻Im= −𝐹Im

(9)on (𝜔)on (𝜔).
E𝐸I + 𝐸m + ℏ𝜔F𝑌Im (𝜔) + 𝛿𝐻Im= −𝐹Im

(10)Here, (𝑋Im , 𝑌Im ) are the forward and backward amplitudes, which we should tono(on)determine by solving Eqs. (9) and (10) . 𝛿𝐻Im(𝜔) are the residual fields

which linearly depend on (𝑋Im , 𝑌Im ). The main idea of the FAM is to calculate theresidual fields without expansion with respect to (𝑋Im , 𝑌Im ). Instead, we use an
iterative solution of the equations with the finite difference method. This will
facilitate programing computer codes and reduce the computational cost. See Refs.
[9,10] for the detailed formulation and [11] for numerical implementation.

From nuclear structure to nuclear reaction

The results of nuclear structure calculations can be used for nuclear reaction

models. For instance, the nucleon density distribution can be used for calculationof the total reaction cross sections [13]. In the optical limit approximation of theGlauber model, the phase shift function as a function of the impact parameter 𝒃
is given by𝑒 fy(𝒃) = exp }− ~ 𝑑𝒓• 𝑑𝒓 € 𝜌• (𝒓• )𝜌€ (𝒓 € )Γ‚‚ (𝒔• − 𝒔€ + 𝒃)„
where 𝒔•(€) are the transverse component of the projectile (target) coordinates,Γ‚‚ (𝒔• − 𝒔€ + 𝒃) are the nucleon-nucleon profile function, and 𝜌•(€) (𝒓) are theprojectile (target) density distributions. The total reaction cross section is given
as integration over 𝒃,
𝜎… = / 𝑑𝒃𝑃(𝒃)with the reaction probability𝑃(𝒃) = 1 − W𝑒 fy(𝒃) WnThe density distributions can be also used for calculation of the folding potentials.
𝑈(𝒓) = ~ 𝑑𝒓• 𝑑𝒓 € 𝜌• (𝒓• )𝜌€ (𝒓€ )𝑣‚‚ (𝒓 − 𝒓• + 𝒓 € )The nucleon-nucleus optical potentials can be calculated in a similar way [14].

Numerical calculations with 3D coordinate-space representation

To perform the numerical calculations with (TD)DFT, we must choose the basis
to represent each (quasi-)particle state. The choice of the basis is arbitrary, but
depends on the purpose and the adopted energy density functional. We use theSkyrme energy density functional, and it is convenient to adopt the threedimensional (3D) coordinate space discretized in mesh as the basis. Therefore,
the quasi-particle wave functions in Eq. (2) are given by vectors,
𝑈I= (1)ˆ ⋮ ‹,
𝑈I= (𝑀)

𝑉I= (1)ˆ ⋮ ‹,
𝑉I= (𝑀)(11)

where 𝑀 is twice of the number of discretized mesh points. The factor of two
comes from the spin degrees of freedom. Accordingly, the single-particle

Hamiltonian ℎ= and the pair fields Δ= are given by 𝑀 × 𝑀 matrices. However,
due to the zero-range nature of the Skyrme energy density functional, the matrix
is almost diagonal, except for kinetic terms and a few derivative terms. Therefore,
the coordinate-space representation is an effective way to perform the calculations.
The time-dependent simulation can be done with the same coordinate and realtime representation. The time is also discretized with a mesh Δ𝑡, and the timedependent quasi-particle states at 𝑡f = 𝑖 × Δ𝑡 are expressed as vectors in Eq. (3):𝑈I= (1, 𝑡f )
ˆ‹,
⋮
=(𝑀,)
𝑈I𝑡f

𝑉I= (1, 𝑡f )
ˆ‹.
⋮
=(𝑀,)
𝑉I𝑡f(12)In the canonical-basis TDHFB equations (4) and (5), the same representationin the coordinate space and the real time is adopted for the canonical basis states.
See Refs. [7,8] for details.
In the linear response equation, the forward and backward amplitudes are also
represented in the mesh space. Instead of time, they are functions of frequency 𝜔.
Thus, we also discretize 𝜔 and solve the FAM equations (9) and (10). For theiterative solver, we use either the generalized conjugate residual method or themodified Broyden method. See Refs. [9,10,11] for details.
The calculations of properties provided in InPACS and Dens have been
performed as follows:(1) Ground-state deformation and other properties shown [8,12] (InPACS)
Method: HF+BCS
Adopted energy density functional: SkM*+monopole pairing
Representation: 3D coordinate mesh representation
Symmetry restriction: None
(2) Proton and neutron density distributions [14] (InPACS)
Method: HFB
Adopted energy density functional: Sly4+mixed pairing
Representation: 1D coordinate mesh representation

Symmetry restriction: Spherical shape(3) Photoabsorption cross sections [8,12]
Method: Cb-TDHFB
Adopted energy density functional: SkM*+monopole pairing
Representation: 3D coordinate mesh representation
Symmetry restriction: None
(4) Proton and neutron density distributions [8,12] (Dens)
Method: HF+BCS
Adopted energy density functional: SkM*+monopole pairing
Representation: 3D coordinate mesh representation
Symmetry restriction: NoneReferences[1] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
[2] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
[3] N. N. Bogoliubov, Sov. Phys. JETP 7, 41 (1958).
[4] M. Bender and P. H. Heenen and P.-G. Reinhard, Rev. Mod. Phys. 75, 121(2003).
[5] T. Nakatsukasa, K. Matsuyanagi, M. Matsuo, and K. Yabana, Rev. Mod. Phys.
88, 045004 (2016).

[6] E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984).
[7] S. Ebata, T. Nakatsukasa, T. Inakura, K. Yoshida, Y. Hashimoto, and K.
Yabana, Phys. Rev. C 82, 034306 (2010).
[8] S. Ebata, T. Nakatsukasa, T. Inakura, Phys. Rev. C 90, 024303 (2014).
[9] T. Nakatsukasa, T. Inakura, and K. Yabana, Phys. Rev. C 76, 024318 (2007).
[10]

P. Avogadro and T. Nakatsukasa, Phys. Rev. C 84, 014314 (2011).

[11]

K. Washiyama and T. Nakatsukasa, Phys. Rev. C 96, 041304 (2017).

[12]

S. Ebata and T. Nakatsukasa, Phys. Scr. 92, 064005 (2017).

[13]

W. Horiuchi, T. Inakura, T. Nakatsukasa, Y. Suzuki, Phys. Rev. C 86,024614 (2012).
[14]

K. Minomo, K. Washiyama, and K. Ogata, Preprint: arXiv:1712.10121.