Density functional approaches to collective phenomena in nuclei: Time-dependent density-functional theory for perturbative and non-perturbative nuclear dynamics

# Density functional approaches to collective phenomena in nuclei: Time-dependent density-functional theory for perturbative and non-perturbative nuclear dynamics

Takashi Nakatsukasa
###### Abstract

We present the basic concepts and our recent developments in the density functional approaches with the Skyrme functionals for describing nuclear dynamics at low energy. The time-dependent density-functional theory (TDDFT) is utilized for the exact linear response with an external perturbation. For description of collective dynamics beyond the perturbative regime, we present a theory of a decoupled collective submanifold to describe for a slow motion based on the TDDFT. Selected applications are shown to demonstrate the quality of their performance and feasibility. Advantages and disadvantages in the numerical aspects are also discussed.

## 1 Introduction

### 1.1 Nucleus as a quantum object

The nuclei provide the matters with mass, the stars with fuel, and the universe with a variety of elements. It was discovered by Ernest Rutherford and coworkers about 100 years ago Rut1911 [], that explains the large-angle scattering of alpha particles by a gold foil GM1909 []. This discovery was also a beginning of the era of the quantum mechanics. Rutherford estimated the upper limit of the nuclear size which turned out to be much smaller than that of the atom. In the atomic scale (Å), the nuclear size (fm) could be regarded as just a point! According to the classical mechanics, the atoms must collapse into nuclear size, because the attractive Coulomb potential eventually brings all the electrons into the nucleus. This mystery stimulated Niels Bohr to develop his idea on the quantum mechanics Bohr1913 [].

#### 1.1.1 Atoms and molecules

With a knowledge of the quantum mechanics, it is easy to understand why the atoms do not collapse to the nuclear scale. If an electron was confined in the nuclear scale of femtometer, the uncertainty principle immediately tells us that its zero-point kinetic energy would become gigantic (order of GeV). In order to decrease this kinetic energy to a reasonable magnitude, the electron’s wave function must have an atomic size of Å. Thus, the atomic size is a consequence of the quantum effect.

The molecules (and the solids) are a good contrast to atoms. The atom is bound by the Coulomb interaction, whose range is infinite (), between a positively charged nucleus and electrons with a negative charge. Since the molecules consist of these charge-neutral atoms, the interaction between a pair of neutral atoms does not have the long-range tail of , but normally has a short-range repulsive part and an intermediate-range attractive part (Fig. 1(a)). Because of this charge neutrality, the molecule is easy to disintegrate into smaller units (molecules and atoms). The atomic size is approximately constant and independent from the atomic number, while the molecular size varies depending on the number of atoms and their kinds. Last but not the least, in the zero-th order approximation, the ground states of the molecules can be classically described as atoms located at fixed relative positions. The hindered quantum fluctuation in molecules is simply due to the fact that the atomic mass, which is approximately identical to the nuclear mass, is about 2,000 times larger than the electronic mass.

#### 1.1.2 Nuclei

The nucleus has a number of properties analogous to the molecules, except for its strong quantum nature. It is a self-bound system composed of fermions of spin and isospin with approximately equal masses, called nucleons (protons and neutrons). The nuclear species are classified by the numbers of neutrons () and protons (). Rutherford discovered that the size of the nucleus is as tiny as the femtometer, but later it was found that the nuclear size varies, as its volume is roughly proportional to the mass number (). Each nucleon is a color-singlet (neutral) object. The interaction between a pair of nucleons (nuclear force) has a finite range of a scale of the pion Compton wave length (). Similar to the molecules, it has a short-range repulsive part and an intermediate-range attractive part (Fig. 1(b)). The nucleus can be disintegrated into small pieces with a small separation energy. In fact, heavy nuclei can have “negative” separation energies because of the repulsive Coulomb energy among protons.

The quantum nature is an important difference between the nucleus and the molecule. In Fig. 1, we show schematic pictures of the atomic and nucleonic potentials. The ground state of the diatomic molecule (panel (a)) is formed at the bottom of the potential, . In the length scale of Å, the atomic mass is heavy enough to localize the wave function at the location of the bottom of the potential, . Thus, the relative distance between a pair of atoms is fixed at . This property allows us to describe the atomic motion in the classical mechanics, such as in the molecular dynamics. In contrast, the nuclear interaction is not strong enough to bind nucleons at the bottom of the potential. In other words, the nucleon’s mass is too light to localize the wave function in the sub-femtometer scale. In deuteron, the zero-point kinetic energy cancels the negative potential energy (), leading to a bound state at approximately zero energy (). The deuteron wave function spatially extends beyond the range of nuclear force (), which reduces the kinetic energy . This shows a striking contrast to the diatomic molecule, and is even analogous to the atoms, that the large size of the deuteron is a consequence of the quantum effect. This strong quantum nature also tells us that the infinite nuclear matter will not be crystallized even at zero temperature, but will stay as the liquid. The nuclear system provides us with unique opportunities to study femto-scale quantum liquids.

### 1.2 Computing nucleus from scratch

The strong quantum nature in finite nuclei leads to a rich variety of unique phenomena. Remarkable experimental progress in production and study of exotic nuclei requires us to construct theoretical and computational approaches with high accuracy and a reliable predictive power. Extensive studies have been made for constructing theoretical models to elucidate basic nuclear dynamics behind a variety of nuclear phenomena. Simultaneously, significant efforts have been made in the microscopic foundation of those models.

For light nuclei, the “first-principles” large-scale computation, starting from the bare nucleon-nucleon (two-body & three-body) forces, is becoming a current trend in theoretical nuclear physics. Among them, the Green’s function Monte Carlo (GFMC) method is the most successful first-principles approach to nuclear structure calculation PW01 []. In this approach, using the Monte Carlo technique, the many-body wave function is sampled in the coordinate space with spin and isospin degrees of freedom. The success of the GFMC clearly demonstrates that we are able to construct a light nucleus from the scratch on the computer. The GFMC method has been applied to nuclei up to the mass number . Another first-principles approach is to project the nuclear Hamiltonian in a truncated Hilbert space, then diagonalize it. This is called no-core shell-model (NCSM) method NNBV06 []. The NCSM also shows successful applications up to the -shell nuclei. The GFMC and NCSM both indicate the exponential increase in computational tasks with respect to the increasing nucleon number. The third approach, the coupled-cluster method (CCM), has an advantage that the required computation increases only in power with respect to the nucleon number. The CCM, which was originally invented in nuclear physics KLZ78 [] and later became extremely successful in quantum chemistry BM07 [], has been revisited as an ab-initio computational approach to nuclear structure Wlo05 []. Especially, it is powerful to study closed-shell nuclei.

Although these first-principles approaches have recently shown a significant progress, they are still limited to nuclei with the small mass number. This may sound mysterious to physicists in other fields, because we know that the similar kinds of approaches are able to treat systems of much larger particle numbers. For instance, using the CCM, nowadays, the chemists can easily calculate a molecular structure with 100 electrons. Why is the first-principles calculation of nuclear structure so difficult? The answer is perhaps trivial for nuclear physicists, but may not be so for others. Let us pick up several important aspects leading to this answer. (1) The strong quantum nature. As we have discussed previously, the full quantum mechanical treatment is necessary for nuclear structure calculation. (2) Strong coupling nature. The nucleon-nucleon scattering length is approximately fm in channel. This is much larger than the mean distance between nucleons inside the nucleus (). (3) Singular property of the nucleonic interaction. To reproduce the phase shift for the nucleon-nucleon scattering, the interaction should have a strong repulsive core at short distance BM69 []. (4) Complexity of the nucleonic interaction. The interaction has a strong state dependence which may be represented by the spin- and isospin-dependence BM69 []. It also contains strong non-central and non-local terms, such as tensor and spin-orbit interactions. Furthermore, for a quantitative description of nuclei, it is indispensable to introduce three-body interactions in addition to the two-body force. (5) Coexistence of different interactions. In addition to the strong interaction, we need to treat the electromagnetic interactions among protons. (6) Coexistence of different energy scales. The nuclear binding energy amounts to order of GeV for heavy nuclei. Because of the strong quantum nature in nuclei, this is a consequence of the cancellation between positive kinetic energy and negative potential energy. Thus, we need to compute these enormous positive and negative components in high accuracy to obtain a correct binding energy. (7) Fermionic nature of nucleons. Needless to say, since the nucleons are fermions, the total wave function must be anti-symmetrized. (8) Finite systems without an external potential. Electronic many-body problems in molecules and solids are solved with external fields produced by nuclei with positive charges. In contrast, the nucleus is a self-bound finite system. Thus, we are usually required to obtain an intrinsic wave function without the center-of-mass degrees of freedom. This requirement often restricts our choice of the basis functions.

Despite of these difficulties, significant advances in the computer power may lead to a realistic “first-principles” construction of the -shell nuclei in near future. There have been extensive efforts toward this direction, especially developing an algorithm suitable for parallel use of the vast number of processors unedf [].

### 1.3 Density functional theory (DFT)

In contrast to the first-principles calculations, which are limited to nuclei with small number of nucleons, the density functional theory (DFT) is currently a leading theory for describing nuclear properties of heavy nuclei BHR03 [], LPT03 []. It is capable of describing almost all nuclei, including nuclear matter, with a single universal energy density functional (EDF). An argument based on the quantum many-body theory leading to nuclear EDF was also developed in 70’s80’s Neg82 [], Nak11 []. In addition, its strict theoretical foundation is given by the basic theorem of the DFT HK64 [], KS65 []. Since the nucleus is a self-bound system without an external potential, we should slightly modify the DFT theorem. This will be discussed in Sec. 2.

The nucleus by itself produces a potential confining nucleons which is analogous to the Kohn-Sham potential in DFT KS65 []. Nuclear physicists often call this potential “mean field”, though it is different from a naive mean-field potential directly constructed from the nucleonic interaction. There are many evidences for the fact that the mean-free path of nucleons inside the nucleus is significantly larger than the nuclear radius BM69 [], in spite of the strong (and even singular) two-body interaction. This is partially due to the Pauli exclusion principle. Since the nucleon mean-free path is roughly proportional to Gal58 [], it is significantly enhanced for nucleons whose energies are close to the Fermi energy. Another even simpler argument was given by Bohr and Mottelson BM69 [], that the nuclear normal density is much lower than the one giving the close packing limit (crystalline limit). In this argument, the quantum effect plays a primary role.

The DFT theorem guarantees the existence of generalized density functionals for every physical observable (see Sec. 2). However, to construct the exact functional, we need experimental data and other theoretical inputs. Currently, there are a variety of EDFs which predict somewhat different properties for nuclei very far from the stability line. We are in want of finding an ultimate universal energy density functional, which is capable of exact description of every nucleus in the nuclear chart. In addition to the recent progress in the first-principles calculations for light nuclei, the radioactive beam facilities in the world will give us ideal opportunities to determine the parameter set for a better functional. New data on neutron halos and skins in medium heavy nuclei may provide important information on its dependence on density and density gradient. Observation of new isotopes in a long isotopic (isotonic) chain may lead to useful constraints on the isovector parts of the energy functional.

An extension of the DFT to the time-dependent DFT (TDDFT) provides a feasible description of many-body dynamics, which contains information on excited states in addition to the ground state. The TDDFT is justified by the one-to-one correspondence between the time-dependent density and time-dependent external potential RG84 [], which will be presented in Sec. 2.4. The TDDFT has vast applications to quantum phenomena in many-body systems. Among them, the perturbative regime has been mostly studied so far. Different approaches to the linear response calculations will be presented in Sec. 3. It is of significant interest but challenging to go beyond the perturbative regime. Nuclei show numerous phenomena related to the large amplitude collective motion, such as fission, shape transition, shape coexistence, anharmonic vibrations, and so on. We present, in Sec. 4, a theory to identify an optimal collective submanifold in the classical phase space of large dimensions.

## 2 Basic formalism for particles in self-bound systems

The density functional theory (DFT) has been extremely successful for calculations of ground-state properties of atoms, molecules, and solids. It describes a many-particle system exactly in terms of its one-body density alone. The DFT is based on the original theorem of Hohenberg and Kohn (HK) HK64 [] which was proved for the ground-state of the many-particle system. Every observable can be written, in principle, as a functional of density.

In nuclear physics, however, we need to treat an isolated system without an external potential. The present nuclear EDF produces a localized density profile without an external potential. Thus, the ground state spontaneously violates the translational symmetry and seems to contain spurious excitation related to the center-of-mass motion. Furthermore, it also violates the rotational symmetry when the ground state is spontaneously deformed. It is interesting to see whether the nuclear EDF can be theoretically justified in a strict sense. We present a possible justification in Sec. 2.1, according to the recent progress Eng07 [], Gir08 [], GJB08 [].

### 2.1 DFT theorem for a wave-packet state

The HK theorem HK64 [] guarantees one-to-one mapping between a one-body density and an external potential . Then, since the ground-state wave function is a functional of density, in principle, all the observables should be functionals of the density as well. However, to describe an isolated self-bound finite system in a box of volume , it is somewhat useless to use the ground-state density in the laboratory frame, because it must be constant, (), where is the particle number. Instead, we want to use a functional of the intrinsic density , where is the center-of-mass coordinate of the total system. In this case, the original HK theorem cannot be justified, because it adopts a one-body external field coupled to the density in the laboratory frame.

Validity of the DFT for the intrinsic states has been discussed by several authors Eng07 [], Gir08 [], GJB08 []. In this section, we present a method to define the DFT for an intrinsic state, more precisely for a “wave-packet” state. The argument here essentially follows the idea by Giraud et al GJB08 [].

In fact, all the nuclear EDFs, currently available, produce a wave-packet state. The minimization of an EDF without the external potential () leads to the nucleon density distribution with a finite radius. This violates the translational symmetry of the original Hamiltonian. Therefore, it is of our interest here to justify the energy functional of the wave-packet density in the laboratory frame in a strict sense.

First, we assume that a wave-packet state in the laboratory frame can be expressed by a product wave function of intrinsic and spurious degrees of freedom, . Here, indicates the intrinsic state and defines the spurious motion. This decomposition can be exactly done for the translational motion.

 H=Hintr(ξ,π)+→P22Nm, (1) Φ(→r1,⋯,→rN)=ϕ(→ξ1,⋯,→ξN−1)⊗χ(→R), (2)

where denote the center-of-mass coordinates (total momenta). and are the relative Jacobi coordinates and their conjugate momenta, respectively. Since the intrinsic ground state, which is supposed to be unique, is completely independent from the spurious motion, we can adopt an arbitrary form of ; e.g., . Then, the ground wave-packet state can be obtained by the variation after the projection:

 (E0)χ≡minχ:fixed[⟨Φ|HP|Φ⟩⟨Φ|P|Φ⟩] (3)

where the projection operator does not change the intrinsic state but makes an eigenstate of the total linear momentum . The variation with respect to the full space () contradicts the uniqueness of the ground state, because states with different give the same . Therefore, the variation here should be performed with respect only to the intrinsic state . This is indicated by the subscript “:fixed” in Eq. (2.1). The wave-packet density profile is simply given by , that depends on the form of . In the followings, we always assume a fixed form of for the wave-packet state .

Next, we introduce an external potential . The following minimization with respect to the intrinsic state leads to the “minimum” energy and defines the wave packet .

 E0[v0]=minχ:fixed[⟨Φ|HP|Φ⟩⟨Φ|P|Φ⟩+⟨Φ|V0|Φ⟩]=minχ:fixed[⟨Φ|HP|Φ⟩⟨Φ|P|Φ⟩+∫ρ(→r)v0(→r)d3r]. (4)

Noted that operates on a state , not on a projected state . does not correspond to the ground-state energy of a system with the Hamiltonian , however, it reduces to Eq. (2.1) for . Now, let us show the one-to-one correspondence between the external potential and the wave-packet density . The proof proceeds by reductio ad absurdum in the same manner as the original proof of HK HK64 []. Assume that another potential , which defines the wave packet , produces the same density . Then, the energy for is given by

 E0[v′0]=⟨Φ′|HP|Φ′⟩⟨Φ′|P|Φ′⟩+⟨Φ′|V′0|Φ′⟩. (5)

If we replace the state by , the energy must increase.

 E0[v′0] < ⟨Φ|HP|Φ⟩⟨Φ|P|Φ⟩+⟨Φ|V′0|Φ⟩ (6) = ⟨Φ|HP|Φ⟩⟨Φ|P|Φ⟩+⟨Φ|V0|Φ⟩+⟨Φ|V′0−V0|Φ⟩ = E0[v0]+∫ρ(→r){v′0(→r)−v0(→r)}d3r.

Interchanging and , we also find

 E0[v0]

Addition of Eqs. (6) and (7) leads to the inconsistency, . This proves the one-to-one correspondence between the external field and the wave-packet density . Thus, both and the wave packet are functionals of the density . In order to lift restriction to -representative densities, we use the constrained search Lev79 [] in which one considers only states that produce a given density , and define the universal functional

 Fχ[ρ]≡minΦ→ρ[⟨Φ|HP|Φ⟩⟨Φ|P|Φ⟩]. (8)

Here, the subscript “” indicates the minimization with a constraint on . The density functional contains the energy of the fixed center-of-mass spurious motion, , that is trivially given by . Therefore, the energy of the ground state with (intrinsic energy) may be obtained by the minimization with a constraint on the total particle number, as

 Eintr(N)≡EP=0(N)=minρ[Fχ[ρ]−λ(∫ρ(→r)d3r−N)]−Eχ(N). (9)

Note that, although the wave function is fixed, may depend on the total mass of the system (particle number ). In principle, the expectation value of any observable , which only depends on the intrinsic degrees of freedom, is a functional of ,

 O[ρ]≡⟨ϕ|^O|ϕ⟩×⟨χ|χ⟩=⟨Φ|^O|Φ⟩, (10)

because the wave-packet state is given as a functional of . Since the state is fixed, there is a trivial correspondence between the wave-packet density and the intrinsic density . This completes the basic theorem of the DFT for the wave packet.

### 2.2 The Kohn-Sham (KS) scheme

For a many-body system of fermions, the shell effects play a major role to determine the ground state. In other words, we need a density functional which takes account of the kinetic energy properly. This is known to be difficult in the local density approximation RS80 []. At present, a scheme given by Kohn and Sham KS65 [] only provides a practical solution for this problem. Here, we follow the same argument.

We introduce a reference system which is a “virtual” non-interacting system with an external potential . This reference system is supposed to reproduce the same density of the “physical” interacting wave packet, but does not have to reproduce the center-of-mass wave function . The ground state of the reference system is trivially obtained as a Slater determinant constructed by the solution of222 Precisely speaking, the orbitals are not necessarily the eigensolutions of Eq. (11), but arbitrary as far as they give the same Slater determinant. We come back to this gauge freedom in Sec. 2.5.

 (−12m∇2+vs(→r))ϕi(→r)=ϵiϕi(→r), (11)

adopting the unit , and the density is given by . The kinetic energy333 The HK theorem guarantees that of the reference system is a functional of the density. is given by

 Ts[ρ]=N∑i=1⟨ϕi|(−12m∇2)|ϕi⟩. (12)

The variation of the total energy of the reference system

 Es[ρ]=Ts[ρ]+∫vs(→r)ρ(→r)d3r, (13)

with a constraint on the particle number, , leads to the following equation:

 μ=δTs[ρ]δρ(→r)+vs(→r). (14)

Although we did not explicitly construct as a functional of , the solution of Eq. (14) must be identical to the solution of Eqs. (11) and (12).

The success of the Kohn-Sham (KS) scheme comes from a simple idea to decompose the kinetic energy in the physical interacting system into two parts; , which is a major origin of the shell effects, and the rest, which is treated as a part of “correlation energy” described by a simple functional of density,

 Fχ[ρ]=Ts[ρ]+Ec[ρ], (15)

where . Then, the variation of leads to Eq. (14) but the potential is now a functional of density, defined by . The only practical difference between the reference system and the interacting system is that, since is a functional of density in the latter, Eq. (11) must be self-consistently solved. These equations are called Kohn-Sham (KS) equations. The self-consistent solution of the KS equations provides the density of a wave-packet state with a fixed corresponding to a (local) minimum of the EDF, . The success of the KS scheme is attributed to the goodness of the local density approximation for .

### 2.3 Open issues

#### 2.3.1 Subtraction of the center-of-mass energy

In the proof of the basic theorem for the wave packet, given in Sec. 2.1, we need to fix a wave function of the center-of-mass motion . The choice of this spurious wave function is arbitrary, but the energy depends on this choice. In practice, the subtraction of is normally performed by constructing the state from the obtained result. This is somewhat inconsistent with the assumption of the fixed center-of-mass state . This could be easily corrected by taking of a given . However, as far as we know, this has not been examined yet.

#### 2.3.2 Validity of the Kohn-Sham scheme

The KS scheme is to take into account a major part of the kinetic energy as , and the rest as a correction. In other words, the KS scheme implicitly assumes that the energy is able to be well approximated by a simple functional of . In fact, this question is still an open issue, not only in the nuclear physics but also in other quantum many-body systems. In the present wave-packet theory, the kinetic energy of the wave packet, , depends on the center-of-mass state . Therefore, there may be an optimal choice for to minimize the difference . The question about the validity of the Kohn-Sham scheme remains to be answered.

#### 2.3.3 Non-spherical wave packet

The nuclear EDFs are known to produce a spontaneous symmetry breaking about the rotational symmetry. Namely, we often encounter a deformed wave-packet density, which accounts for appearance of the rotational spectra in nuclei. For instance, many experimental evidences indicate that nuclei in the rare-earth region and in the actinide region are deformed BM75 []. According to the argument in Sec. 2.1, we may separate the rotational motion from the intrinsic degrees of freedom, then, we have

 Φ(→r1,⋯,→rN)≈ϕ(ξ)⊗χ(→Θ,→R), (16)

where indicates angle variables. Then, replacing the operator by that of angular momentum projection, the DFT for deformed wave packet can be shown in the same manner. However, in contrast to the translational motion, the separation of the rotation degrees of freedom is not exact. but only approximate. Thus, there remains an ambiguity for the definition of the functional (8): Namely, the minimization must be performed in the entire space except for the subspace that accounts for translational and rotational correlations. In this sense, the use of the EDF, which produces a deformed state, can be justified only approximately.

### 2.4 Time-dependent density functional theory

The DFT is designed for calculating the ground-state properties. For excited-state properties and reactions, the time-dependent density functional theory (TDDFT) is a powerful and useful tool. In this section, we recapitulate the basic theorem for the time-dependent density functional theory (TDDFT).

Since the proof of the HK theorem is based on the Rayleigh-Ritz variational principle, its extension to the time-dependent density is not straightforward. This was done by Runge and Gross RG84 [], showing that there is one-to-one correspondence between a time-dependent density and a time-dependent external potential . The external potential is required to be expandable in a Taylor series about the initial time ,

 v(→r,t)=∞∑k=01k!vk(→r)(t−t0)k. (17)

The external potentials, and , are defined to be different if there exist some minimal nonnegative integer such that where . In other words, and differ more than a time-dependent function, . Note that the potentials differing by the time-dependent constant produce the same density because the corresponding wave functions differ by a merely time-dependent phase, as in Eq. (27) with .

Now, let us assume that starting a common initial state , two different external potentials, and , produce densities and , respectively. From this, we first prove that the current densities, and , are different. Using the current density operator , the equation of motion is written as

 i∂∂t→j(→r,t)=⟨Ψ(t)|[→^j(→r),^H(t)]|Ψ(t)⟩, (18)

where . We have the same equation for , with and replaced by and , respectively. Then, we have

 ∂∂t{→j(→r,t)−→j′(→r,t)}∣∣∣t=t0 = −i⟨Ψ0|[→^j(→r),^H(t0)−^H′(t0)]|Ψ0⟩ (19) = −1mρ(→r,t0)∇w0(→r).

If , it is easy to see that and are different at . In case that and , we need to further calculate derivative of Eq. (18) with respect to .

 i∂2∂t2→j(→r,t)∣∣∣t=t0=⟨Ψ0|[→^j(→r),∂^H(t0)∂t]|Ψ0⟩+⟨Ψ0|[[→^j(→r),^H(t0)],^H(t0)]|Ψ0⟩. (20)

The second term of Eq. (20) vanishes for , because . Thus,

 ∂2∂t2{→j(→r,t)−→j′(→r,t)}∣∣∣t=t0=−1mρ(→r,t0)∇w1(→r)≠0. (21)

Again, we can conclude that at . In general, if for and , we repeat the same argument to reach

 (∂∂t)n+1{→j(→r,t)−→j′(→r,t)}∣∣∣t=t0=−1mρ(→r,t0)∇wn(→r)≠0. (22)

Therefore, there exists a mapping from the expandable potential to the current density .

Next, we use the continuity equation

 ∂∂t{ρ(→r,t)−ρ′(→r,t)}=−∇⋅{→j(→r,t)−→j′(→r,t)}, (23)

and calculate the (n+1)-th derivative of Eq. (23) at . From Eq. (22),

 (24)

Provided , we can prove that the right-hand side of Eq. (24) does not vanish identically. This is done by using the following identity:

 (25)

The integral of both sides of Eq. (25) over the entire space leads to

 −∫d→rwn(→r)∇⋅{ρ(→r,t0)∇wn(→r)}=∫d→rρ(→r,t0){∇wn(→r)}2>0, (26)

where we assume is localized in space so that the surface integral vanishes. Therefore, from Eq. (24), we can conclude that the densities and are different at . This completes the proof.

The time-dependent density determines the time-dependent external potential except for the time-dependent constant. Therefore, the many-body time-dependent state should be a functional of density except for a time-dependent phase.

 |Ψ(t)⟩=exp(−iα(t))|Ψ[ρ](t)⟩. (27)

Any observable quantity must be independent from the global phase, , thus, a unique functional of density, . Note that these time-dependent density functionals depend on the initial many-body state .

### 2.5 Time-dependent Kohn-Sham (TDKS) equations

In practice, we use the Kohn-Sham scheme KS65 [] for numerical calculations. Assuming the -representability, the time-dependent Kohn-Sham (TDKS) equation is given by

 i∂∂tψi(→r,t)={−12m∇2+vs[ρ](→r,t)}ψi(→r,t),i=1,⋯,N. (28)

The density of a system is expressed by . In practice, we usually adopt the potential same as the one for calculation of the ground state (“adiabatic approximation”), except for the external potential .

 vs[ρ](→r,t)=v(→r,t)+δE[ρ]δρ(→r)∣∣∣ρ (29)

The density is invariant with respect to the unitary transformation among occupied KS orbitals. Therefore, there are gauge degrees of freedom to choose this transformation at any instant of time. For explicit notification of the gauge freedom, it is convenient to introduce the matrix notation as follows. Let be an arbitrary single-particle basis set, and we define a matrix of size of , as . Then, the density matrix can be written as

 ραβ(t)=∑i⟨α|ψi(t)⟩⟨ψi(t)|β⟩=(¯ψ(t)¯ψ†(t))αβ. (30)

The orthonormal property of the KS orbitals is expressed as . Denoting the TDKS Hamiltonian in Eq. (28) as , the TDKS equations (28) can be generalized into the following form.

 i∂∂t¯ψ(t)=hs(t)¯ψ(t)−¯ψ(t)ξ(t), (31)

where is an arbitrary time-dependent Hermitian matrix which represents a generator of the transformation. Equation (31) is equivalent to the well-known equation for the density matrix Eba10 [].

 i∂∂tρ(t)=[hs(t),ρ(t)]. (32)

The stationary state corresponds to the time-indenpendent density, .

### 2.6 Pairing correlations: Kohn-Sham-Bogoliubov (KSB) equations

The HK theorem (or its wave-packet version), in principle, guarantees that the energy of the system can be exactly written as a functional of density, . However, in practice, it is often difficult to take into account all the correlation energy in , solely depending on . The kinetic energy is such an example, which is resolved by the genius idea by Kohn and Sham. The pairing correlation energy , which is important for heavy nuclei in open-shell configurations, is another example difficult to be expressed by only.

To overcome this difficulty, a common strategy is to extend the KS equations, according to the Bogoliubov’s quasiparticles Bog58 []. Each orbital now has two components, , and its number is basically infinite (). These are called quasiparticle (qp) orbitals. The previous KS equations are extended to the following equations, which we call Kohn-Sham-Bogoliubov (KSB) equations444 Again, this corresponds to a special choice in the gauge degrees of freedom. See Sec. 2.7. hereafter:

 (Hs−μN)Φν=EνΦν, (33)

where

 Hs≡(h[ρ,κ]Δ[ρ,κ]−Δ∗[ρ,κ]−h∗[ρ,κ]),N=(100−1). (34)

The Hamiltonian is in the same form as that in Eq. (11), , while the pair potential is introduced to describe the pairing correlations. The same form of equation is known as the Hartree-Fock-Bogoliubov equation in nuclear physics RS80 [], BR86 []. Now, the KSB Hamiltonian in Eq. (34) not only depends on the density , but also on the pair density . In Eq. (33), the matrix convention is assumed. Namely, when we adopt a single-particle basis of , we have and , and and correspond to Hermitian and anti-symmetric matrices, respectively.

The chemical potential is determined so as to satisfy the following number condition: . With this number constraint, must have a finite value. On the other hand, the values of the pair density are determined solely by the variation of the total energy. Therefore, the self-consistent solution of the KSB equations (33) spontaneously produces the finite values of and .

The solutions of the KSB equation have a “paired” property: If the qp state is a solution of Eq. (33) with eigenvalue , the qp state is also a solution with eigenvalue . We call “unoccupied” qp orbitals, and “occupied” qp orbitals DFT84 []. This naming is based on the generalized density matrix,

 R=(ρκ−κ∗1−ρ∗) (35)

which is Hermitian and idempotent: . The “unoccupied” (“occupied”) qp orbitals correspond to the eigenvectors of with eigenvalue 0 (1); and .

Denoting the dimension of the single-particle Hilbert space as , we may define the matrix as follows.

 Φαν(t)={⟨α|Uν⟩α=1,⋯,M⟨α−M|Vν⟩α=M+1,⋯,2M (36)

which represents “unoccupied” qp orbitals (). The “occupied” orbitals with size of are defined in the same manner, with () replaced by (). The generalized densities are expressed in terms of these matrices as . The orthonormal property of the qp orbitals is given by . Combining the ‘unoccupied” and “occupied” orbitals to construct the matrix , the matrix becomes a unitary matrix RS80 [].

In the energy functional of Skyrme type VB72 [], the pairing correlation energy is simply added to the original energy functional.

 E[ρ,κ]≡E[ρ]+Epair[ρ,κ], (37)

that depends only on the local densities. Therefore, Eq. (33) becomes local in coordinates. However, in general, the KSB Hamiltonian, and , are not necessarily local. For instance, the Gogny functional DG80 [] gives non-local KSB equations.

### 2.7 Time-dependent Kohn-Sham-Bogoliubov (TDKSB) equations

For a time-dependent description, the inclusion of the pair density leads to the time-dependent Kohn-Sham-Bogoliubov (TDKSB) equations. They can be formulated in a matrix form analogous to Eq. (31). Using an arbitrary Hermitian matrix , the TDKSB equations may be written as

 i∂∂tΨ(t)=Hs(t)Ψ(t)−Ψ(t)Ξ(t), (38)

where the TDKSB Hamiltonian is given by Eq. (34), which depends on time through the densities and . Here, represent time-dependent “unoccupied” qp orbitals (). The “occupied” orbitals are defined in the same manner, with () replaced by (). The TDKSB equation (38) holds for , as well.

Analogous to the stationary case, the generalized density can be written as , and the orthonormal property of the qp orbitals is given by the unitarity of the matrix . The “unoccupied” (“occupied” ) correspond to the subspace with eigenvalue 0 (1); and . In the generalized density matrix formalism, the TDKSB equation is written in an analogous manner to Eq. (32):

 i∂∂tR(t)=[Hs(t),R(t)]. (39)

So far, we have shown similarities between Eqs. (31)(32) and (38)(39). However, there is an important difference between Eq. (32) and Eq. (39). The stationary solution in Eq. (32) corresponds to . In contrast, it is not the case in Eq. (39), . Let us examine this difference in details. The TDKSB equation (38) can be recast into another form, convenient for taking its stationary limit. First, let us factor out the time-dependent phases as follows: and . Here and hereafter, we denote the remaining parts of the quantities as the “primed” ones. The generalized density becomes

 R(t)=¯Ψ(t)¯Ψ†(t)=exp(−iμNt)R′(t)exp(+iμNt), (40)

where

 (41)

Namely, the transformation does not change the density , but modifies the pair density as . Since the pair potential is usually a linear functional of , the same time-dependent phase should be assumed for as well: . The Hamiltonian is transformed in the same manner:

 Hs(t)=exp(−iμNt)H′s(t)exp(+iμNt). (42)

With these primed quantities, the TDKSB equation (38) can be rewritten as

 i∂∂tΨ′(t)={H′s(t)−μN}Ψ′(t)−Ψ′(t)Ξ(t), (43)

or equivalently, in the generalized density matrix,

 i∂∂tR′(t)=[H′s(t)−μN,R′(t)]. (44)

It is now clear that the stationary solution corresponds to , not to , with a proper choice for the parameter identical to the chemical potential. In Eq. (43), it corresponds to with a choice of the gauge matrix . It should be noted that the generalized density is invariant with respect to the choice of the gauge matrix . In contrast, the time-dependent phase factor in and have a physical origin and cannot be removed by the gauge choice. In fact, it is a boost transformation, , from the laboratory frame of reference to the body-fixed frame. The stationary solution with () corresponds to a time-dependent solution in the laboratory frame:

 Ψν(t)=(e−iμt00e+iμt)Ψ′ν. (45)

This is a collective motion associated with the spontaneous generation of the pair density, called pair rotation. Therefore, in terms of the TDKSB formalism, the appearance of the chemical potential in the stationary KSB equation (33) comes from the boost transformation to the body-fixed frame rotating in the gauge space. This is analogous to the appearance of the cranking term in the spatially rotating frame of reference RS80 []. In the case of pair rotation, since the particle number is finite , the system is rotating in the gauge space, even at the ground state. This rotation affects the intrinsic modes of excitation, thus, the Hamiltonian in the rotating frame, , should be utilized to calculate the intrinsic excitation spectra. This point will be discussed again in Sec. 4.5.

## 3 Perturbative regime: Linear response

The theorem of the TDDFT tells us that the functional may depend on the initial state, . This additional ambiguity can be removed by assuming that the initial state is identical to the ground state. With this assumption, the linear response theory with a weak time-dependent perturbation is formulated in this section. The formulation is basically identical to the one known as the random-phase approximation in nuclear physics RS80 [], BR86 []. However, according to the concept of the TDDFT, the theory gives the exact linear density response, with no approximation involved, in principle.555 In practice, some approximations are involved, such as the adiabatic approximation of Eq. (29).

### 3.1 Time-dependent linear density response

We consider a system subject to a time-dependent external potential

 v(→r,t)={0t<0v1(→r,t)t≥0 (46)

in addition to the static potential of the unperturbed system.666 For isolated nuclear systems, we have . In this section, we use the notation of the four vector . We assume that the system is at the ground state at times . Thus, the initial density at can be obtained from the self-consistent solution of the Kohn-Sham equations (11). The first-order density response, , is given by

 ρ1(x)=∫d4x′Π(x,x′)v1(x′) (47)

with the density-density response function

 Π(x,x′)=δρ(x)δv(x′)∣∣∣v=0. (48)

The right-hand side of Eq. (48) is a well-defined quantity, since the basic theorem of TDDFT in Sec. 2.4 guarantees that the time-dependent density is a functional of the time-dependent external potential; .

For non-interacting particles moving in an external potential of , there is a one-to-one correspondence between the time-dependent density and the potential. Therefore, we have

 ρ(x)=ρ[vs](x),vs(x)=vs[ρ](x). (49)

The density-density response function for the non-interacting system is given by

 Πs(x,x′)=δρ(x)δvs(x′)∣∣∣vs[ρ0]. (50)

The potential is written as a sum of a given external potential and the rest of the Kohn-Sham potential, . For instance, in the adiabatic approximation of Eq. (29), . Therefore, using the chain rules, Eq. (48) can be connected to its non-interacting :

 Π(x,x′) = ∫d4yδρ(x)δvs(y)∣∣∣vs[ρ0]⋅δvs(y)δv(x′)∣∣∣v=0 (51) = ∫d4yΠs(x,y){δ(y−x′)+∫d4y′δvks(y)δρ(y′)∣∣∣ρ0δρ(y′)δv(x′)∣∣∣v=0} = Πs(x,x′)+∫d4y∫d4y′Πs(x,y)w(y,y′)Π(y′,x′),

where the residual kernel is given by

 w(x,x′)≡δvks(x)δρ(x′)∣∣∣ρ0. (52)

In the adiabatic approximation, this is equal to the second derivative of the energy functional.

 w(x,x′)≡δ2E[ρ]δρ(x)δρ(x′)∣∣∣ρ0. (53)

Most of the energy functionals currently available are local in time, which leads to .

Multiplying the Dyson-type equation (51) by the perturbing external potential leads to the linear density response of Eq. (47).

 ρ1(x)=∫d4x′Πs(x,x′)vscf(x′), (54)

where the self-consistent effective field, given by

 vscf(x)=v1(x)+∫d4y w(x,y)ρ1(y)