Role of small-norm components in extended random-phase approximation

The role of the small-norm amplitudes in extended RPA theories such as the particle-particle and hole-hole components of one-body amplitudes and the two-body amplitudes other than two particle - two holes components are investigated for the one-dimensional Hubbard model using an extended RPA derived from the time-dependent density-matrix theory. It is found that these amplitudes cannot be neglected in strongly interacting regions where the effects of ground-state correlations are significant.


Introduction
The random-phase approximation (RPA) is formulated based on the Hartree-Fock (HF) ground state where it is assumed that the lowest single-particle states (hole states) are completely occupied and other single-particle states (particle states) are empty. Therefore, RPA consists of particle (p)-hole (h) and h-p amplitudes. When the effects of ground-state correlations are explicitly taken into account in extended RPA approaches such as the renormalized rRPA [1], the self-consistent RPA (SCRPA) [2,3], the extended second RPA (ESRPA) [4], its response function version [5] and the small amplitude limit of the time-dependent densitymatrix theory (TDDM) [6][7][8], the single-particles states become partially occupied. As a consequence, the p-p and h-h amplitudes in addition to the p-h and h-p amplitudes should also be included as one-body transition amplitudes to satisfy sum rules and properties of spurious states [9][10][11]. The two-body amplitudes other than the 2p-2h and 2h-2p amplitudes may also play some roles when the effects of ground-state correlations are included. In most applications, however, these additional amplitudes which have small eigenvalues of a norm matrix have usually been neglected. In this paper we investigate how the p-p and h-h amplitudes and the 2p-2p, 2h-2h and ph-ph amplitudes are incorporated into physical states using an extended RPA approach (ERPA) based on the ground state in TDDM. We use the one-dimensional Hubbard model at half filling which in momentum space has an equal number of particle and hole states and thus allows us to investigate the effects of the above-mentioned additional amplitudes. The paper is organized as follows the formulation is given in sect. 2, the results are shown in sect. 3 and sect. 4 is devoted to summary.

Formulation
We consider the Hamiltonian H consisting of a one-body part and a two-body interaction where a † α and a α are the creation and annihilation operators of a particle at a timeindependent single-particle state α.

Time-dependent density-matrix theory and ground state
The formulation of TDDM has been presented in Refs. [6,7,12]. To be self-contained, we briefly explain it below. TDDM gives 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 αβα ′ β ′ . These matrices are defined as where |Φ(t) is the time-dependent total wavefunction |Φ(t) = exp[−iHt]|Φ(t = 0) and . We use units such that = 1. The equations of motion for n αα ′ and C αβα ′ β ′ are derived from and are written as where ǫ αα ′ is given by Here the subscript A means that the corresponding matrix is antisymmetrized. The matrix B αβα ′ β ′ in Eq. (7) does not contain C αβα ′ β ′ and describes 2p-2h and 2h-2p excitations, while P αβα ′ β ′ and H αβα ′ β ′ contain C αβα ′ β ′ and describe p-p (and h-h) and p-h correlations to infinite order, respectively [7]. These matrices are given in Ref. [7]. The matrix T αβα ′ β ′ is given by 2/11 where C αβγα ′ β ′ γ ′ is the correlated part of a three-body density-matrix, which is given by Here, A is an operator which properly antisymmetrizes n αα ′ ρ βγβ ′ γ ′ under the exchange of the single-particle indices such as α ↔ β and α ′ ↔ β ′ . The three-body correlation matrix is neglected in the original version of TDDM [6,7]. Instead of neglecting C αβγα ′ β ′ γ ′ we use the approximation where p and h refer to particle and hole states, respectively, and N is given by The above expression can be derived by perturbative consideration [12] using the coupledcluster doubles wavefunction. Equations (6) and (7) satisfy the conservation laws of the total energy and the total number of particles [6,7]. The ground state in TDDM is given as a stationary solution of the time-dependent equations (Eqs. (6) and (7)). In this work we use the following adiabatic method to obtain a nearly stationary solution [13,14]: Starting from the Hartree-Fock (HF) configuration, we solve Eqs. (6) and (7) gradually increasing the strength of the residual interaction such as v(r − r ′ ) × t/T . This method is motivated by the Gell-Mann-Low theorem [15] and has often been used to obtain approximate ground states [16,17]. To suppress oscillating components which come from the mixing of excited states, we must take large T .

Extended RPA
The ERPA equation can be derived by either taking the small amplitude limit of the TDDM equations or using the equation of motion approach [1]. It is written in matrix form for the one-body and two-body amplitudes x µ αα ′ and X µ αβα where ω µ is the excitation energy of an excited state µ. The matrices A, B, C and D are given by where |Φ 0 is the ground state in ERPA but approximated by that in TDDM. The norm matrices S 1 , T 1 , T 2 and S 2 are given by The effects of ground-state correlations are included in the above matrices through the occupation matrix and the two-body and three-body correlation matrices. Although the relations A = A † and C † = B hold, D = D † is not satisfied because a condition from the Jacobi identity [19], which is given by in the above equation includes three-body operators, the stationary condition for the three-body density matrix |Φ 0 = 0 is required but it is not satisfied by the approximate threebody correlation matrix given by Eqs. (11) and (12). Thus the Hamiltonian matrix of Eq. (14) is not Hermitian. As a consequence some eigenstates of Eq. (14) can have complex eigenvalues and negative transition strength. It does not cause any serious problems in the applications shown below as long as the interaction is not so strong and the approximation for the three-body correlation matrix is meaningful.
In the following we explain how small-norm components become non-negligible. The matrix S 1 is given by In HF where n αα = 0 or 1 the diagonal elements of S 1 are 1 for the p-h component and −1 for the h-p component. When the single-particle states are fractionally occupied, the p-p and h-h elements of S 1 become non-vanishing though they are small. As a consequence the contribution of the p-p and h-h amplitudes to Eq. (14) become non-negligible. The matrix S 2 includes n αα ′ , C αβα ′ β ′ and C αβγα ′ β ′ γ ′ . When C αβα ′ β ′ and C αβγα ′ β ′ γ ′ are neglected for simplicity, the diagonal element of S 2 is given by [19] where we assume that n αα ′ = δ αα ′ n α . In HF S 2 is not vanishing only for the 2p-2h and 2h-2p configurations: S 2 is 1 (−1) for the 2p-2h (2h-2p) configurations. When the single-particle states are fractionally occupied, other two-body configurations such as ph-ph, 2p-2p and 2h-2h configurations have non-vanishing values of S 2 though they are small, and can contribute to Eq. (14). The matrices T 1 and T 2 , which can couple the p-p and h-h amplitudes to the 2p-2h and 2h-2p amplitudes and the p-h and h-p amplitudes to the ph-ph amplitudes, are given by C αβα ′ β ′ .

4/11
To explain the relation of ERPA to rRPA and SCRPA, we explicitly show the matrix A: where ǫ αα ′ and n αα ′ are assumed to be diagonal for simplicity. The first two terms on the right-hand side of Eq. (22) are of the same form as those in the RPA equation, the next two terms with C αβα ′ β ′ describe the self-energies of the α − α ′ configurations due to ground-state correlations [20], and the last four terms with C αβα ′ β ′ are interpreted as the modification of the interaction [20]. The one-body sector of Eq. (14) Ax µ = ω µ S 1 x µ is formally the same as the equation in SCRPA [2,3,20]. In rRPA C αβα ′ β ′ is neglected. In rRPA n αα is selfconsistently calculated from the p-h and h-p amplitudes, and both n αα and C αβα ′ β ′ are selfconsistently obtained from the p-h and h-p amplitudes in SCRPA. We perform calculations using rRPA-like and SCRPA-like formulations where n αα and C αβα ′ β ′ are not given by the p-h and h-p amplitudes but replaced by those of the TDDM ground state.

Results
To test our approach for a case where a few h and p states are involved, we chose the onedimensional (1D) Hubbard model with periodic boundary conditions and compare with the results in exact diagonalization approach (EDA). In momentum space the Hamiltonian is given by [21] H = k,σ ǫ k a † k,σ a k,σ where U is the on-site Coulomb matrix element, σ spin projection and the single-particle energies are given by ǫ α = −2t cos k α with the nearest-neighbor hopping potential t. We consider the case of the six sites at half filling. In the first Brillouin zone −π ≤ k < π there are the following wave numbers The single-particle energies are ǫ 1 = −2t, ǫ 2 = ǫ 3 = −t, ǫ 4 = ǫ 5 = t and ǫ 6 = 2t. In HF the lowest 6 states are fully occupied as shown in Fig. 1.

Ground state
The ground state energy calculated in TDDM (circles) using the adiabatic method is shown in Fig. 2

Spin mode with momentum transfer |q| = π/3
Since the h-h and p-p transitions such as k 1 → k 2 and k 6 → k 5 can contribute when the single-particle states are partially occupied, we first consider the spin mode with the momentum transfer q = π/3 which can be excited mainly by the one-body operator 7/11 a † k4↑ a k2↑ − a † k4↓ a k2↓ . The excitation energies calculated in RPA (open triangles), rRPA-like (triangles) and SCRPA-like (squares) formulations are shown in Fig. 4 as functions of U/t. These approaches do not have the coupling to the two-body amplitudes. The results in EDA are given with the solid line. Since U is a repulsive interaction, spin modes where the single-particle transitions between spin-down states destructively interfere with those between spin-up states become soft with increasing U . The results in RPA are in reasonable agreement with the exact solutions. In the rRPA-like calculations the two states appear below E/t < 2. The lower state at E/t ≈ 1 mainly consists of the p-p and h-h transitions and the higher state consists of the p-h and h-p components. The transition strength carried by the lower state increases with increasing U . Thus in the rRPA-like approach the excited states consisting of the p-p and h-h components appear as if they are physical ones. This indicates that the inclusion of the ground-state correlation effects only through the fractional occupation of the single-particle states is not appropriate. In the SCRPA-like calculations the states mainly consisting of the p-p and h-h components move to the high energy region (E/t > 10) due to the terms in Eq. (22) with C αβα ′ β ′ . This is because these terms are divided by the small values n pp − n p ′ p ′ or n hh − n h ′ h ′ when Ax µ = S 1 x µ is diagonalized. Thus in the SCRPA-like approach the p-p and h-h components are mixed with the p-h and h-p components but their main components are energetically separated from the low-lying state. However, the excitation energies calculated in the SCRPA-like approach are significantly larger than the exact values, suggesting the importance of the coupling to the two-body amplitudes.
Now we discuss the coupling to the two-body amplitudes. Since the coupling to the 2p-2h amplitudes is large, the results in SRPA (crosses) deviate strongly from the RPA values with increasing U and SRPA collapses at U/t = 3.9 for this spin mode. The results in ERPA are shown with the filled circles. For this spin mode the three-body correlation matrix used in Eq. (14) at U/t = 3.5 and 4 is reduced by 16 % from the value given by Eqs. (11) and (12). The original three-body correlation matrix gives the results shown with the open circles at U/t = 3.5 and 4. The approximation Eqs. (11) and (12) seems to overestimate the threebody correlation matrix for large U . For U/t ≤ 3 the excitation energies in ERPA are little affected by such 16 % reduction of the three-body correlation matrix. The coupling to the two-body amplitudes plays a role in bringing down the results of the SCRPA-like approach. The 2p-2h and 2h-2p amplitudes also play a dominant role in ERPA as in SRPA: The 2p-2h and 2h-2p amplitudes alone can bring the spin mode down to the position slightly below the second excited state calculated in the rRPA-like approach. To clarify the contribution of the p-p and h-h amplitudes, we calculate the transition strength in ERPA at U/t = 4 in two cases using an operatorQ q = k σa † k+q,σ a k,σ : In one case (S full ) all components of x µ αα are included and in the other case (S ph ) only the p-h and h-p components of x µ αα are included. We found S ph /S full = 0.54. Thus the small-norm components have significant contribution to the transition strength.

Spin mode with |q| = π
Next we consider a spin mode with momentum transfer |q| = π, which consists of the singleparticle transitions from the h states with k = ±π/3 to the p states with k = ∓2π/3. In the Hubbard model the spin mode with |q| = π becomes the first excited state for U > 0. 8 Since only the transitions between the single-particle states with k = ±2π/3 and k = ∓π/3 can make |q| = π, there is no contribution to this mode from the p-p and h-h transitions.
The obtained results are shown in Fig. 5. The spin modes in RPA (open triangles) and SRPA(crosses) collapse beyond U/t = 2.3. Since the number of the 2p-2h configurations with |q| = π is small, the coupling to the 2h-2h states is small for this mode. In the rRPA-like approach (filled triangles) which includes the effects of fractional occupation of the singleparticles the point where the spin mode collapses goes slightly up to U/t = 2.6. This means that fractional occupation of the single-particle states is not sufficient to suppress excess p-h correlations. In the SCRPA-like approach which includes n αα and C αβα ′ β ′ the spin mode is stable but the excitation energy deviates from the exact values with increasing U . In ERPA (filled circles) the coupling to the two-body amplitudes plays a role shifting down the result of the SCRPA-like approach but the 2p-2h and 2h-2p amplitudes are not important in contrast to the case of the spin mode with |q| = π/3: The ERPA results calculated only with the 2p-2h and 2h-2p amplitudes are shown in Fig. 5 with the open circles. Thus the two-body amplitudes which have small eigenvalues of S 2 play a role in this spin mode.

Two-phonon states with q = 0
Finally we consider a state with q = 0 excited by two-body operators such as a † k4↑ a † k5↓ a k3↓ a k2↑ . Since the two-body operator is written as , the state with q = 0 corresponds to the two-phonon states of the spin modes with |q| = π/3 and π discussed above. In HF the matrix elements of S 2 for the 2p-2p, 2h-2h and phph amplitudes vanish (see Eq. (21)) and SRPA does not have such amplitudes. In ERPA these amplitudes contribute to Eq. (14). Since the two-phonon state with q = 0 has the same quantum number as the ground state, we neglect X µ αβα ′ β ′ with odd number of p and h states,  to be consistent with the treatment of the ground state. The excitation energy of the twophonon state in ERPA (circles) is shown in Fig. 6 as a function of U/t. The results in EDA are given with the solid line. The results in ERPA slightly overestimate the excitation energy for large U but are in good agreement with the exact values. In SRPA the excited energy goes down with increasing U and deviates from the exact values. The difference between the results in ERPA and SRPA demonstrates the importance of the effects of ground-state correlations.
The strength function calculated in ERPA (solid line) at U/t = 4 for the operatorQ 2 q is shown in Fig. 7, whereQ q is the excitation operator for the spin mode with |q| = π. The dotted line depicts the ERPA result where only the 2p-2h and 2h-2p components are included 10/11 as the two-body amplitudes. The result in EDA is given with the dot-dashed line. The distributions are smoothed with artificial width Γ/t = 1/20. The full ERPA result has the strength comparable to the exact solution, while the ERPA result calculated with only the 2p-2h and 2h-2p amplitudes has quite small strength. The coupling of the 2p-2h amplitudes to the 2h-2p amplitudes is essential to increase the collectivity of two-phonon states [22]. The ph-ph amplitude plays an important role in such a coupling because the number of the ph-ph configurations is much larger than that of the 2p-2p and 2h-2h configurations.

Summary
The role of the small-norm amplitudes such as the particle-particle and hole-hole components of the one-body amplitudes and the two-body amplitudes other than the two particletwo holes amplitudes were investigated for the one-dimensional Hubbard model using an extended RPA (ERPA) derived from the time-dependent density-matrix theory. It was found that these amplitudes cannot be neglected in strongly interacting regions (U/t > 1) where the effects of ground-state correlations are significant. In realistic applications of ERPA, however, we are usually forced to truncate the two-body amplitudes because the number of their elements increases rapidly with increasing number of single-particle states.