Applications of Time-Dependent Density-Matrix Approach

The equations of motion for reduced density matrices form a coupled chain known as the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy. To close the coupled chain at the two-body level, approximations for a three-body density matrix with one-body and two-body density matrices are needed. The time-dependent density-matrix theory (TDDM) assumes that the three-body density matrix is given by the antisymmetrized products of the one-body and two-body density matrices. In this review the truncation schemes of the BBGKY hierarchy beyond TDDM are discussed and a formulation for the study of excited states which is derived from the time-dependent approach is explained. The truncation schemes and the formulation for excited states are applied to the Lipkin model and the Hubbard model to corroborate their validity. Two realistic applications of the TDDM approaches are also presented. One is the dipole and quadrupole excitations of 40Ca and 48Ca and the other the fusion reactions of 16O + 16O.


INTRODUCTION
The time-dependent Hartree-Fock theory (TDHF) is the basis of the mean-field theories such as the Hartree-Fock theory (HF) and the random-phase approximation (RPA): The HF ground state is given as a stationary solution of the TDHF equation and RPA can be formulated as the small amplitude limit of the TDHF equation. HF and RPA have extensively been used as standard theories to study nuclear structure problem [1]. Extensive TDHF simulations have also been performed for heavy-ion collisions [2,3]. However, most experimental data suggest that beyondmean field theories are required for a more realistic description of nuclear structure and reactions. In this paper an approach to extend the mean-field theories based on the equations of motion for reduced density matrices is reviewed. The equations of motion for reduced density matrices form a coupled chain known as the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy [4] where the time evolution of an n-body density matrix depends on n-body and n + 1-body density matrices. The advantage of such a time-dependent density-matrix approach (TDDMA) is that it is directly connected to TDHF in the lowest-level approximation: The truncation of the BBGKY hierarchy at the level of the one-body density matrix by approximating the twobody density matrix with the antisymmetrized product of the one-body density matrices gives TDHF. A beyond TDHF is obtained by the truncation of the BBGKY hierarchy at the two-body level and it needs approximations for the three-body density matrix. A few truncations schemes have been proposed. The simplest truncation scheme is to replace the three-body density matrix with the antisymmetrized products of the one-body and two-body density matrices, neglecting the correlated part of the three-body density matrix [the three-body correlation matrix (C 3 )] [5,6]. This truncation scheme has been called the time-dependent density-matrix theory (TDDM). TDDM has been applied to heavy-ion collisions [7][8][9] and collective excitations [6,[10][11][12]. A simplified TDDM called TDDM P [13] where two-body correlations are consider only for a pair of time reversed single-particle (s.p.) states have also been employed to simulate heavy-ion collisions [13,14]. It was found that TDDM sometimes overestimates ground-state correlations in a solvable model [15], gives unphysical s.p. occupation probabilities in dynamical simulations [16,17] and causes divergent dynamical behaviors in highly excited and (or) strongly interacting cases [18,19]. Obviously the problems originate in the neglect of C 3 and are related to the loss of N-representability [20] which refers to the properties of reduced density matrices derived from an N-body total wavefunction and is completely fulfilled only in the case of the untruncated BBGKY hierarchy. To remedy the difficulties of the naive truncation scheme in the description of ground states, an approximation for C 3 has been proposed based on perturbative consideration [21], where C 3 is given by the traced products of the correlated part of the two-body density matrix (C 2 ). This truncation scheme is referred to as TDDM1 hereafter. The correlation matrices C 2 and C 3 are also called the two-body and three-body cumulants and the above approximation for C 3 corresponds to taking the leading-order terms of the three-body cumulant [22,23]. It has been demonstrated that TDDM1 improves the TDDM results for the ground states of model Hamiltonians [21,23] and also 16 O [24]. In the case of the Lipkin model [25], however, TDDM1 was found to overestimate C 3 in strongly interacting regions. There, another truncation scheme [26] where C 3 in TDDM1 is divided by a reduction factor was introduced. This truncation scheme is referred to as TDDM2. In this paper three truncation schemes TDDM, TDDM1 and TDDM2 are explained and their applications are presented.
The small amplitude limit of the TDDMA equations gives an extended RPA (ERPA) which is used for the study of excited states, as is the case of RPA which is formulated as the small amplitude limit of TDHF. RPA and ERPA are also formulated by using the equation of motion approach [27,28]. ERPA consists of the coupled equations for the one-body and two-body amplitudes and include the effects of ground-state correlations through the fractional occupation probability n α of a s.p. state α and C 2 . ERPA is related to so far proposed beyond-RPA theories. When the coupling to the two-body amplitude is omitted, the ERPA equation for the one-body amplitude is the same as the selfconsistent RPA (SCRPA) [29,30] equation which includes both n α and C 2 . The neglect of C 2 in the SCRPA equation corresponds to the renormalized RPA (rRPA) [27,28] which includes the ground-state correlation effect via n α . When the HF ground state is assumed, the equations in ERPA are reduced to those in the second RPA (SRPA) [31]. If particle-hole correlations included in the two-body amplitudes were expressed by phonons, ERPA would be connected the particle-vibration coupling or quasiparticle-phonon models [32]. ERPA has been applied to solvable models [23,33]. Realistic cases have also been studied in ERPA [34][35][36]. Main results of the ERPA applications are presented in this paper.
The TDDM truncation scheme has been used for simulations of heavy-ion collisions [7][8][9] where the TDDM equations are formulated by using the time-dependent s.p. states which obey a TDHF-like equation. The application of such a TDDM approach to the fusion reactions of 16 The paper is organized as follows: The equations of motion for the one-body and two-body density matrices formulated by using a time-independent s.p. basis are given in section 2 and the truncation schemes of the BBGKY hierarchy and the formulation of ERPA are discussed. The applications are made to two model Hamiltonians, the Lipkin model [25] and the onedimensional Hubbard model [37] in section 2 and the obtained results are compared with the exact solutions. The ERPA results for the dipole and quadrupole excitations in 40 Ca and 48 Ca are presented in section 2 as realistic applications of ERPA. The TDDM formulation using a time-dependent s.p. basis and its application to the fusion reactions of 16 O + 16 O are given in section 3. Section 4 is devoted to summary and outlook.

FORMULATION IN TIME-INDEPENDENT SINGLE-PARTICLE BASIS
The TDDMA equations are formulated for an N fermion system described by the Hamiltonian H consisting of a one-body part t (the kinetic energy term) and a two-body interaction v where a + α and a α are the creation and annihilation operators of a particle at a s.p. state α and the s.p. states are assumed time-independent.

Time-Dependent Density-Matrix Theory and Truncation Schemes
The TDDMA equations given in Tohyama and Schuck [21] are explained below. They consist of the coupled equations of motion for the one-body density matrix (the occupation matrix) n αα ′ and the correlated part of the two-body density matrix C αβα ′ β ′ (C 2 ). These matrices are defined as where | (t) is the time-dependent total wavefunction ) and A is an operator which properly antisymmetrizes n αα ′ n ββ ′ under the exchange of the s.p. indices such as α ↔ β and α ′ ↔ β ′ . Unitsh = 1 are used hereafter. The equations of motion for n αα ′ and C αβα ′ β ′ are derived from by evaluating the commutation relations. They are written as where ǫ αα ′ is the s.p. energy including the mean field and is given by Here the subscript A means that the corresponding matrix is antisymmetrized. The term B αβα ′ β ′ in Equation (6) consists of only the occupation matrices and describes 2 particle (p) -2 hole (h) and 2h-2p excitations, while P αβα ′ β ′ and H αβα ′ β ′ contain C 2 and express p-p (and h-h) and p-h correlations to infinite order, respectively [6]. The T αβα ′ β ′ term gives the coupling to the three-body correlation matrix (C 3 ) where C αβγ α ′ β ′ γ ′ (C 3 ) is given by Approximations for C 3 are needed to close the equations of motion within n αα ′ and C 2 . Three truncation schemes TDDM, TDDM1 and TDDM2 have been proposed. In TDDM, C 3 is simply omitted [5,6]. In TDDM1, C 3 is given by [21] C p 1 p 2 h 1 p 3 p 4 h 2 = h C hh 1 p 3 p 4 C p 1 p 2 h 2 h , where p and h refer to particle and hole states, respectively. These expressions were derived from perturbative consideration using the Coupled-Cluster-Doubles (CCD)-like ground state wavefunction [38]. In a time-independent density-matrix approach in quantum chemistry, known as the contracted Schrödinger equation [20], Mazziotti [22] has proposed a method for constructing the three-body cumulant (C 3 ) with n αα ′ and C 2 . Equations (10) and (11) describe the leading-order terms in the three-body cumulant [23]. TDDM2 [26] is the most effective in the large N and strongly interacting limits of the Lipkin model and gives where N is The factor N was introduced to simulate many-body effects which reduce C 3 in large N systems and (or) strongly interacting regions of the Lipkin model. In the perturbative region where the second term on the right-hand side of Equation (26) is smaller than unity, N has the meaning of the normalization factor of the total wavefunction. The conservation of the total energy and total particle number is not affected by the truncation schemes for C 3 as long as its symmetry and anti-symmetry properties under the exchange of s.p. indices is respected. However, the trace relation between the one-body and two-body density matrices n αα ′ = λ ρ αλα ′ λ /(N − 1) is not conserved when any approximation is made for C 3 . This is an example of the loss of N-representability [20]. It was pointed out [21] that the fulfillment of the trace relation is drastically improved by going from TDDM to TDDM1. In an attempt to conserve the trace relation, Cassing and Pfitzner [39] proposed an approximation for C 3 which also contains quadratic terms of C 2 . However, C 3 is not uniquely determined only by the requirement of the trace relation conservation. In contrast to TDDM1 their quadratic terms do not have the leading-order terms (Equations 10, 11) of the three-body cumulant [22,23] and the dynamical effect of C 3 was found small in one-dimensional heavy-ion simulations [39]. C 3 in Reference [39] is not anti-symmetric under the exchange of s.p. indices, which may violate even the conservation of the trace relation as was pointed out by Gherega et al. [17]. There is another attempt [40] to conserve the trace relation, where the equation motion for C 3 was solved by truncating the BBGKY hierarchy at the three-body level. However, the application of such an approach was limited to model Hamiltonians [40].

Ground-State Calculation
The ground state in TDDMA is given as a stationary solution of the time-dependent equations (Equations 5, 6) which satisfieṡ n αα ′ = 0 andĊ 2 = 0. Two methods have been employed to obtain the stationary solution. One is the adiabatic method: Equations (5) and (6) are solved by starting from the HF configuration and gradually increasing the strength of the residual interaction such as v( r − r ′ ) × t/T. This method is based on the Gell-Mann-Low theorem [41] and has often been used to obtain approximate ground states with various time-dependent functionals [11,13,14,42,43]. To suppress oscillating components which come from the mixing of excited states, T must be much larger than the longest period in the system considered. The other method is a usual iterative gradient method which is useful to obtain a rigorously stationary solution.
Since it involves matrix inversion, the application of the gradient method is limited to small systems: The gradient method has been employed to obtain the ground states of the Lipkin model [15] and the oxygen, calcium and tin isotopes [34][35][36] using several s.p. states around the Fermi level [34,36] or the valence neutron s.p. states [35].

Excited-States Calculation
The formulation for excited states can be derived by either taking the small amplitude limit of the TDDM equations or using the equation of motion approach [1,27]. Here the formulation based on the equation of motion approach is presented. Let us consider a generalized RPA operator with one-body and two-body sectors where there is no restriction on the s.p. indices: They can be both p and h. As usual with the equation of motion approach, the properties of the excitation operator Q + µ | 0 = | µ and Q µ | 0 = 0 are assumed and the following equations of motion satisfied by exact states are taken into account where ω µ is the excitation energy of an excited state | µ . The equations for x µ λλ ′ and X µ λ 1 λ 2 λ ′ 1 λ ′ 2 are obtained by inserting Equation (15) into the above equations. They are written in matrix form [15,40] The Hamiltonian matrices A, B, C and D are given by The norm matrices S 1 , T 1 , T 2 , and S 2 are given as These matrices are evaluated by assuming | 0 to be the ground state in TDDMA. This means that the effects of ground-state correlations are included in the above matrices through n αα ′ and C 2 . All matrix elements in Equation (18) are given in Tohyama and Schuck [40]. The one-body sector of Equation (18), Ax µ = ω µ S 1 x µ , is explicitly shown below to explain how n αα ′ and C 2 are included. The matrix S 1 is given by and A by where ǫ αα ′ and n αα ′ are assumed to be diagonal for simplicity. The first two terms on the right-hand side of Equation (22) are of the same form as the RPA and rRPA euqations, the next two terms with C 2 and the Kronecker delta δ αα ′ describe the self-energies of the α-α ′ configurations due to ground-state correlations [29,44], and the other terms with C 2 are interpreted as the vertex corrections [29,44]. Equation (18) has the most general form of beyond RPA theories: It is reduced to SRPA when the ground-state correlations are neglected and the onebody sector of Equation (18) Ax µ = ω µ S 1 x µ is formally the same as the equation in SCRPA. Equation (18) is referred to as the extended RPA (ERPA) hereafter. Although ERPA has great advantages over other extended RPA theories, it is worth pointing out its limitations. The numbers of the matrix elements of C 2 and X µ αβα ′ β ′ increase rapidly with increasing number of the s.p. states. Therefore, truncation of the s.p. space is required in realistic applications. As shown below, basic effects of two-body correlations can be described with rather small s.p. space, however. The other limitation is concerned with hermiticity of D in Equation (18), which is related to the truncation of the BBGKY hierarchy. Equation (21) contains C 3 and it is approximated depending on the truncation scheme. Hermiticity of D which is guaranteed only when all the matrix elements of C 3 satisfy the stationary condition as those of n αα ′ and C 2 do is not fulfilled [40] when any approximated is made for C 3 . The non-hermiticity has not caused serious problems in the applications thus far considered, though. In the case of the Lipkin model, an attempt [40] to obtain Hermitian D was carried out by solving the equation of motion for C 3 .

Applications
TDDMA's in the time-independent s.p. basis have been applied to model Hamiltonians [23,33] to corroborate their validity: The Lipkin model [25] was used to compare three truncation schemes TDDM, TDDM1 and TDDM2. Comparison of ERPA with other beyond-RPA theories, rRPA, SCRPA and SRPA was performed for the one-dimensional Hubbard model [37]. As realistic applications of ERPA, the quadrupole excitations of oxygen isotopes [34], the low-lying quadrupole states in tin isotopes [35] and the dipole and quadrupole excitations of 40 Ca and 48 Ca [36] have been studied. In the following some of the results are presented.

Lipkin Model
The Lipkin model [25] has extensively been used to test theoretical models. It describes an N-fermions system with two N-fold degenerate levels. The upper (lower) levels have energies ǫ/2 (−ǫ/2) and quantum number p (−p) with p = 1, 2, ..., N. The Hamiltonian is given by where the operators are the followings For χ = |V|(N − 1)/ǫ ≤ 1 the HF ground state is given by |HF = N p=1 c + −p |0 , where |0 is the true vacuum. For χ > 1 the lowest s.p. states are obtained by the transformation where α satisfies cos 2α = 1/χ. The HF ground state in this case is often called the "deformed" HF (DHF) state and is given by |DHF(α) = N p=1 a + −p |0 . The truncation schemes TDDM, TDDM1 and TDDM2 have been applied to the Lipkin model and it was found that the simplest scheme TDDM gives the exact solutions in the limits of large N and χ [45]. This is due to the unique property of the Lipkin model that the ground-state energy in DHF becomes exact in such limits [1]. The relation between the density-matrices in DHF and TDDM is discussed below. The occupation matrix in DHF is given by n −p−p = cos 2 α, n pp = 1 − n −p−p = sin 2 α, and n p−p = cos α sin α. The two-body and three-body density matrices in DHF are given by the above elements of the occupation matrix. For example the 2p-2h and ph-ph elements of the two-body density matrix are expressed as Similarly, the three-body density matrix is given by This means that the correlated part (C p−p ′ p ′′ pp ′ −p ′′ ) of the threebody density matrix vanishes in DHF. The "spherical" total wavefunction | in DHF which does not have the mixing of the p and −p states is given by the two DHF solutions as Since the overlap between |DHF(α) and |DHF(−α) is negligibly in the large χ and N limits, the three-body density matrix calculated with | has also no correlated part. This is the reason why the results in TDDM approach the exact solutions in the large χ and N limits. The ground-state energy E 0 calculated in TDDM (solid line) for N = 12 and 50 is presented respectively in Figures 1, 2 as a function of χ. The dashed and green (gray) lines denote the results in TDDM1 and TDDM2, respectively. In the case of N = 50 the results in TDDM2 are not displayed because they lie between the TDDM results and the exact values. The dotted and dot-dashed lines depict the results in HF and DHF (χ > 1) and the exact values, respectively. As seen in Figures 1, 2, the factor N in Equations (12) and (13) plays a role in greatly reducing C 3 , making TDDM2 almost equivalent to TDDM for N = 50. In the limits of large N and χ both TDDM and DHF results become close to the exact solutions. In the transitional region χ ≈ 1, however, TDDM1 and TDDM2 were found better than TDDM and DHF [45]. In the case of N = 12 this extends to χ ≈ 2 as seen in Figure 1.
The excitation energies of the first and second excited states calculated from the small oscillations of the TDDM solutions (dots) are compared with the exact solutions (dot-dashed line) in Figures 3, 4 for N = 200 where TDDM is supposed to give the nearly exact ground state. The dotted lines in Figures 3, 4 depict the results in RPA and the "deformed" RPA (χ > 1) for the first excited state. In contrast to RPA and the deformed RPA TDDM reproduces the smoothly decreasing excitation energy of the first excited state with increasing interaction strength beyond χ = 1. Figure 4 shows that TDDM also gives good description of the second excited state. As seen in Figure 4, the excitation energies for the first excited state calculated in the deformed RPA become close to the exact values for the second excited state [1] with increasing χ. The excited states were also calculated for a small system with N = 4 by using the ground states in TDDM [15] and TDDM1 [23], and it was found that the TDDM1 ground state gives much better results.
As was pointed out above, the fact that TDDM becomes exact in large N and χ limits is due to the unique feature of the Lipkin model that the mean-field theory DHF gives the exact solutions in such limits. In the transitional region χ ≈ 1 TDDM1 and TDDM2 give better description of the exact solutions than TDDM, and the applications to the ground states of other solvable models [23,26] and 16 O [24] also showed that TDDM1 largely improves TDDM whereas the improvement  (10) and (11) are shown with the dashed line. The green (gray) line depicts the results in TDDM2 where the three-body correlation matrix is given by Equations (12) and (13). The results in HF and DHF(χ > 1) are depicted with the dotted line. The exact values are given by the dot-dashed line. Adapted from Tohyama and Schuck [45] with permission from Società Italiana di Fisica/Springer-Verlag GmbH Germany. from TDDM1 to TDDM2 is not large [26]. Therefore, TDDM1 or TDDM2 may be a useful truncation scheme to be applied to realistic cases except for strongly interacting regions.

One-Dimensional Hubbard Model
ERPA based on the TDDM1 ground state has been applied to the one-dimensional (1-D) Hubbard model [37] to compare with other beyond RPA theories [33]. The Hubbard model is one of the most widespread models to investigate strong electron correlations and has often been used to corroborate the validity of beyond RPA theories [46]. In momentum space the Hamiltonian is given by [46] where U is the matrix element of the on-site repulsive Coulomb interaction, σ spin projection and the s.p. energies are given by ǫ α = −2t cos k α with the nearest-neighbor hopping potential t.
In the case of the six sites at half filling considered here, there are the following six wave numbers The s.p. energies are ǫ 1 = −2t, ǫ 2 = ǫ 3 = −t, ǫ 4 = ǫ 5 = t, and ǫ 6 = 2t. In HF the lower states with ǫ 1 , ǫ 2 and ǫ 3 are fully occupied by 6 particles. The ground state energy calculated in TDDM1 (dots) with the adiabatic method is displayed in Figure 5 as a function of U/t. It was found [21] that TDDM1 gives much better ground states of this model Hamiltonian than TDDM and that the improvement from TDDM1 to TDDM2 is minor [26]. The exact values obtained in exact diagonalization approach are depicted with the solid line. The TDDM1 results agree well with the exact values. The occupation probabilities of the four s.p. states in TDDM1 (circles and squares) given in Figure 6 as functions of U/t also have reasonable agreement with the exact solutions (solid, dotted, dashed and dot-dashed lines). The deviation of the occupation probabilities from the HF values (n αα = 1 or 0) exceeds 10% at U/t = 4.
ERPA is compared with other beyond RPA theories for the spin mode with the momentum transfer q = π/3 which is excited mainly by the one-body operator a + k 4 ↑ a k 2 ↑ − a + k 4 ↓ a k 2 ↓ . Since the s.p. states are partially occupied, the h-h and p-p transitions such as k 1 → k 2 and k 6 → k 5 also contribute in rRPA, SCRPA and ERPA. In Figure 7 the excitation energies in ERPA (dots), RPA (open triangles), rRPA (filled triangles), SCRPA (squares) and SRPA (crosses) are shown as functions of U/t. The exact solutions are given with the solid line. The rRPA and SCRPA results are calculated with n αα and C αβα ′ β ′ which are not self-consistently determined by the p-h and h-p amplitudes [28,29] but given by the TDDM1 calculations. In the case of a repulsive interaction, the excitation energy of a spin mode where the s.p. transitions between spin-up states and spin-down states destructively interfere decreases with increasing U. The results in RPA agree rather well with the exact solutions. In rRPA there are two states below E/t < 2. The main components of the lower state at E/t ≈ 1 are the p-p and h-h transitions and the higher state consists of the p-h and h-p components. Thus in rRPA the configurations consisting of the p-p and h-h components appear as the lowest state as if it is a physical state. This indicates that it is not appropriate to include the ground-state correlation effects only through n α . In SCRPA the states originating from the p-p and h-h transitions gain self energies and move to the high energy region (E/t > 10). This is because the terms in Equation (26) with C 2 are divided by the small values n pp − n p ′ p ′ or n hh − n h ′ h ′ when Ax µ = S 1 x µ is solved. Thus in SCRPA the states consisting of the p-p and h-h components are energetically separated from the lowest state. The excitation energies of the lowest state calculated in SCRPA, however, exceed significantly the exact values. This is due to the neglect of the coupling to the two-body amplitudes. SRPA includes the coupling to the twobody amplitudes though the ground-state correlation effects are neglected. As shown in Figure 7 with the crosses, the coupling to the 2p-2h amplitudes considered in SRPA downwardly shift the RPA results with increasing U and SRPA collapses at U/t = 3.9. The results in ERPA have reasonable agreement with the exact solutions. The coupling to the two-body amplitudes included in ERPA plays a role in bringing down the results in SCRPA.
From the application to the 1-D Hubbard model it was clarified that when the couping to the two-body amplitude is considered, the effects of ground-state correlations should also be included, and vice versa. Therefore, rRPA, SCRPA and SRPA which only include either the ground-state correlations effects or the coupling to the two-body amplitude cannot give satisfactory results.

Dipole and Quadrupole Excitations of 40 Ca and 48 Ca
In this subsection the applications of ERPA to the dipole excitation of 48 Ca and the quadrupole excitation of 40 Ca [36] are presented. It is demonstrated that the effects of ground-state correlations which are not fully incorporated in other beyond RPA theories play a significant role in the fragmentation of transition strengths.
Since the numbers of C 2 and X µ αβα ′ β ′ increase rapidly with the number of the s.p. states, rather sever truncation of the s.p. space is required in realistic applications. The occupation probability n αα and C 2 were calculated within TDDM by using the truncated s.p. basis consisting of the 2s-1d and 1f -2p states, and only the 2p-2h and 2h-2p elements of C 2 were included to reduce the dimension size. It was pointed out in Reference [24] that TDDM with this truncation of C 2 gives as good results for the ground state of 16 O as TDDM1 with all components of C 2 . The Skyrme III force [47] was used to obtain the s.p. Observed occupation probabilities [48] are shown for the 1d 3/2 states.
wavefunctions which satisfy a HF-like n αα -dependent equation. In Pfitzner et al. [11] and Peter et al. [12] a fixed harmonic oscillator basis was chosen to facilitate the calculations of twobody matrix elements when the TDDM equations were applied to the study of giant resonances. A simplified interaction which contains only the t 0 and t 3 terms of the Skyrme III force was used as the residual interaction. The ground states were calculated with the iterative gradient method [15]. The one-body amplitudes x µ αα ′ in Equation (18) were defined with a large number of s.p. states including those in the continuum to satisfy the energy-weighted sum rule: The continuum states were discretized by confining the wavefunctions in a sphere with radius 15 fm and all the s.p. states with ǫ α ≤ 50 MeV and j α ≤ 11/2h were included. The residual interaction in Equation (18) was assumed to have the same form as that used in the ground-state calculations. Since the residual interaction differs from the effective interaction used in the calculation of the s.p. states, it is necessary to reduce the strength of the residual interaction. The reduction factors 0.66 and 0.69 for 40 Ca and 48 Ca, respectively, were determined so that the spurious mode corresponding to the center-of-mass motion has zero excitation energy in RPA. To reduce the number of the two-body amplitudes, only the 2p-2h and 2h-2p components of X µ αβα ′ β ′ were considered for the 2s-1d and 2p-1f states. The occupation probabilities calculated in TDDM for 40 Ca and 48 Ca are given in Tables 1, 2, respectively. They deviate more than 10% from the HF values (n αα =1 or 0) in 40 Ca. This indicates that the ground state of 40 Ca is highly correlated as an RPA study [49] and perturbative calculations [31,50] have already shown. Since the occupation of the neutron 1f 7/2 state in 48 Ca blocks some 2p-2h excitations, the ground-state correlations are weaker in 48 Ca than in 40 Ca. As will be discussed below, the fractional occupation of the 2p-1f states plays an important role in the fragmentation of dipole and quadrupole transition strengths. Occupation probabilities deduced from ground-state-to-groundstate (p, d) and (e, e ′ p) reactions [48] are also shown for some s.p. states in Tables 1, 2 (values in parentheses). These values also strongly deviate from the HF value (n α = 1). The TDDM results cannot be directly compared with these data, however. A more appropriate formalism such as the odd-particle number RPA [51] which deals with odd particle systems is needed to compare with experiment.
The strength function for the isovector dipole excitation in 48 Ca calculated in ERPA (solid line) is displayed in Figure 8. Observed occupation probabilities [48] are shown for the proton 2s 1/2 and neutron 1f 7/2 states.
The dotted and dot-dashed lines depict the results in RPA and SRPA, respectively. An artificial width Ŵ = 0.5 MeV is used to smooth the distributions. The strength distributions in Figure 8 exhaust about 90% of the energy-weighted-sum-rule (EWSR) value including the enhancement term given by the t 1 and t 2 parameters of the Skyrme III force. A better treatment of the residual interaction and the continuum states is required to fulfill the EWSR value. The sharp peak in RPA corresponds to the giant dipole resonance (GDR). GDR strongly couples to the 2p-2h states and it is damped both in SRPA and ERPA. The occupation of the neutron 1f 7/2 state in 48 Ca allows the 2p-2h states which include the neutron 1f 7/2 state as a hole state. Since these states have energies close to the energy of GDR, GDR is strongly damped due to the coupling to the 2p-2h states. The SRPA result in Figure 8 dose not show a strong downward shift of the dipole strength which has been reported in large scale SRPA calculations [53]. This is due to fact that a rather small number of the s.p. states are used to define the 2p-2h amplitudes. The peak position and width of GDR in ERPA are comparable with the experimental photo absorption cross section [52] as shown in the inset of Figure 8. ERPA gives 7 states below 10 MeV, which are compared with experiment [54] in Figure 9. These states involve the transitions from the partially occupied neutron 2p 1/2 , 2p 3/2 and 1f 5/2 states and the p-p transition components exhaust 15 − 39% of the transition amplitude (x µ , X µ ). The summed strength below 10 MeV is 213 ×10 −3 e 2 fm 2 , which somewhat overestimates the experimental value 61.5 ± 7.8 ×10 −3 e 2 fm 2 . SRPA gives two dipole states at 9.2 MeV and 9.3 MeV with the summed strength 21 × 10 −3 e 2 fm 2 . The study of low-lying dipole strength distribution has been attracting strong interests and various theoretical approaches such as the large scale SRPA [55,56] and the quasi-particle phonon coupling models [57,58] have been successfully applied to calcium isotopes.
The strength function for the isoscalar quadrupole excitation in 40 Ca calculated in ERPA (solid line) is shown in Figure 10.   the strength distributions in Figure 10 exceed the EWSR value by about 10% due to the simple approximations for the residual interaction and the continuum states. The main peak in RPA corresponds to the giant quadrupole resonance (GQR). ERPA brings much larger fragmentation of the quadrupole strength than SRPA especially to the low energy region, indicating the importance of the ground-state correlations effects included in ERPA. The large fragmentation of the quadrupole strength is consistent with experiment [59,60]. The p-p transitions allowed by the fractional occupation of the 1f 7/2 states and the coupling of the 2p-2h amplitudes to the h-h or p-p amplitudes through C 2 were found to play an important role in the large fragmentation of the quadrupole strength [34,36]. The importance of the coupling of the one-body amplitude to C 2 in the fragmentation of GQR in 40 Ca was also pointed out by the 1p-1h⊗phonon configurations model [61,62]. A large scale SRPA calculation [63] shows a downward shift of the quadrupole strength and larger fragmentation of GQR than the SRPA result in Figure 10. This difference again originates from the difference in the number of the 2p-2h configurations used. In the inset the ERPA strength distribution in the GQR region (lower part) is compared with the experimental data from (p, p ′ ) experiments [59] (upper part). Although the peak position in ERPA corresponds to the experimental data, ERPA cannot describe the large fragmentation of GQR. The result of the large scale SRPA calculation [63] suggests the importance of higher configurations.
There are 19 sates below 10 MeV in ERPA, which are compared with experiment [54] in Figure 11. The first 2 + state in 40 Ca cannot be described in RPA and ERPA because it mainly consists of 4p-4h states [64] as in the case of 16 O [65]. The summed strength below 10 MeV is 166 e 2 fm 4 in ERPA, which is about two thirds of the experimental value 263 ± 46 e 2 fm 4 where the first 2 + state is excluded.
From the applications of ERPA to the dipole and quadrupole excitations of 48 Ca and 40 Ca it was clarified that the ground state correlation effects should be included to explain the large fragmentation of the dipole and quadrupole strengths in doubly-magic nuclei. The ground-state correlation effects in magic nuclei have extensively been studied for spin-isospin modes [31,50,[66][67][68].  [54]. Adapted from Tohyama [36] under the Creative Commons CCBY license.

FORMULATION IN TIME-DEPENDENT SINGLE-PARTICLE BASIS
The first applications of the time-dependent density-matrix approach were based on the TDDM truncation scheme and the TDDM calculations [6,7] were performed by using the time-dependent s.p. wavefunctions obtained from the then available TDHF code with axial symmetry and without spin-orbit force [69]. More advanced TDHF codes with spin-orbit force, unconstrained symmetry and improved effective interactions have been used to solve the TDDM equations [8,13]. Since the calculation of the two-body matrix elements is time-consuming, a simpler approximation called TDDM P [13,14] has also been employed in heavy-ion collisions, where two-body correlations are considered only for a pair of time reversed s.p. states to reduce the number of matrix elements of C 2 . The TDDM approaches based on the time-dependent s.p. basis have been applied to study the particle transfers in heavy-ion collisions [7], the fusion reactions [8,9] and the damping of giant resonances at zero [6,70,71] and finite temperatures [10]. TDDM P has been applied to the particle transfers in heavy-ion collisions [13] and the fusion reactions [14]. In the following the TDDM formulation in the time-dependent s.p. basis is given and the application to the fusion reactions of 16 O + 16 O is presented in some detail as an example.

TDDM Equations
The one-body density matrix ρ and the correlated part C 2 of the two-body density matrix ρ 2 are expanded with a finite number of time-dependent s.p. states ψ α ρ(11 where the numbers denote space, spin and isospin coordinates. The equations of motion of TDDM in the time-dependent basis consist of the following three coupled equations [6]: where h is the mean-field Hamiltonian given by ρ. When the time-dependent s.p. states are chosen, the terms with the s.p. energies on the right-hand side of Equations (5) and (6) are incorporated into the equation for the s.p. wavefunctions [6]. In TDDM P [13,14] H αβα ′ β ′ is neglected and two-body correlations are considered only for a pair of time reversed s.p. states.

Fusion Reactions of 16 O + 16 O
The fusion reactions of 16 O + 16 O studied in TDDM with the time-dependent s.p. basis are explained below. This work Tohyama and Umar [8] finally solved the longstanding problem of fusion window anomaly. Early TDHF calculations showed that the colliding heavy ions do not fuse in a small impact parameter region when incident energy is higher than a certain relatively low threshold value E th [2]. This is known as the fusion window anomaly. Experimental search for the fusion window anomaly has found no evidence [72][73][74][75]. It was found that the inclusion of spin-orbit force introduced enough one-body dissipation to 16 O + 16 O collisions [76] because the degeneracy of the 1p 3/2 and 1p 1/2 states is lifted. It was also found [77] that the effects of two-body dissipation taken in TDDM resulted in the doubling of E th without incorporating spin-orbit force. This is due to the inclusion of additional unoccupied s.p. states in TDDM. In Tohyama and Umar [8] the combined effects of spin-orbit force and two-body dissipation were studied for 16 Table 3. The inclusion of either spinorbit force or two-body dissipation dramatically increases E th . However, two-body dissipation increases E th only about 10 MeV when spin-orbit force is included. It was also found that the translational motion damps faster in TDDM than in TDHF [8] below E c.m. = 69 MeV where the colliding system fuses both in TDHF and TDDM. The fusion reactions of 16 O + 16 O below E th were also studied by Wen et al. [14] using the TDDM P approach and a paring interaction as the residual interaction and it was found that extracted friction coefficients are enhanced by about 20% due to two-body dissipation.
In the case of heavy-ion collisions the TDDM (and TDDM P ) equations ostensibly do not conserve the total energy because of the truncation of the s.p. space [7,14]. Wen et al. [14] has proposed a method to recover the energy conservation within the truncated s.p. space.

SUMMARY AND OUTLOOK
An approach which extends the time-dependent Hartree-Fock theory (TDHF) based on the equations of motion for reduced density matrices was presented. The equations of motions for reduced density matrices form a coupled chain known as the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy. In this time-dependent density-matrix approach (TDDMA) the truncation of the BBGKY hierarchy is applied at the twobody level by approximating the correlated part of the threebody density matrix (C 3 ). TDDMA has great advantages that a correlated ground state is obtained as a stationary solution of the TDDMA equations and that the small amplitude limit of the TDDMA equations gives the most general form of beyond the random-phase approximation (RPA). TDDMA was applied to the Lipkin model to test the approximations for C 3 . It was found that the simplest approximation where C 3 is neglected becomes exact in the limits of large number of particles and strong interaction. The extended RPA (ERPA) derived from the TDDMA equations was applied to the one-dimensional Hubbard model to compare with other beyond RPA theories. It was pointed out that when the effects of ground-state correlations are included, the coupling to the two-body amplitudes should also be considered, and vice versa. As the realistic applications of ERPA, the dipole and quadrupole excitations of 40 Ca and 48 Ca were studied. It was found that the effects of ground-state correlations play an important role in fragmenting the dipole and quadrupole strengths. The TDDMA study for the fusion reactions of 16 O + 16 O was also presented as an application of the TDDMA formulation with a time-dependent singe-particle basis. It was pointed out that the two-body dissipation plays a role in further damping the translational motion of 16 Although the obtained results indicate that TDDMA provides a promising beyond mean-field framework to include two-body correlation effects which are missing in TDHF, TDDMA has limitations and issues to be resolved. One limitation is the fact that the number of matrix elements of the two-body density matrix rapidly increases with increasing number of the s.p. states, which forces us to use small s.p. space around the Fermi level in realistic applications, though the obtained results show that basic effects of two-body correlations can be described with rather small s.p. space. In the realistic applications, simple residual interactions of the δ function form were used to facilitate the calculations of two-body matrix elements whereas the s.p. wavefunctions were obtained from the Skyrme interactions included in TDHF codes. The consistent treatment of the effective interactions is a subject to be addressed in TDDMA.
Another point is the truncation scheme of the BBGKY hierarchy itself. In TDDMA the BBGKY hierarchy is truncated at the twobody level by making approximations for C 3 . The truncation violates the properties of reduced density matrices which should be fulfilled if they are derived from an N-body total wavefunction. It was pointed out that the TDDM1 truncation scheme where C 3 is given by the traced products of the 2p-2h elements of the two-body correlation matrix largely improves the simple scheme where C 3 is neglected. However, the validity of TDDM1 cannot be simply extended to highly excited cases such as heavy-ion collisions. The study of the truncation schemes in such cases remains a difficult but interesting subject to be investigated.