Truncation scheme of time-dependent density-matrix approach
aa r X i v : . [ nu c l - t h ] D ec Truncation scheme of time-dependent density-matrix approach
Mitsuru Tohyama
Kyorin University School of Medicine, Mitaka, Tokyo 181-8611, Japan
Peter Schuck
Institut de Physique Nucl ´ e aire, IN2P3-CNRS, Universit ´ e Paris-Sud, F-91406 Orsay Cedex, France andLaboratoire de Physique et de Mod´elisation des Milieux Condens´es, CNRS et,Universit´e Joseph Fourier, 25 Av. des Martyrs, BP 166, F-38042 Grenoble Cedex 9, France
A truncation scheme of the Bogoliubov-Born-Green-Kirkwood-Yvon hierarchy for reduced densitymatrices, where a three-body density matrix is approximated by the antisymmetrized products oftwo-body density matrices, is proposed. This truncation scheme is tested for three model hamilto-nians. It is shown that the obtained results are in good agreement with the exact solutions.
PACS numbers: 21.60.Jz
I. INTRODUCTION
The equations of motion for reduced density matri-ces have a coupling scheme known as the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy wherean n -body density matrix couples to n -body and n + 1-body density matrices. An approximation of a two-body density matrix with the antisymmetrized productsof one-body density matrices in the equation of motionfor the one-body density matrix gives the time-dependentHartree-Fock theory (TDHF). If we approximate a three-body density matrix with the antisymmetrized productsof one-body and two-body density matrices in the equa-tion of motion for the two-body density matrix, we cantruncate the BBGKY hierarchy and obtain the coupledequations of motion for the one-body and two-body den-sity matrices. This is the time-dependent density-matrixtheory (TDDM) [1, 2]. A stationary solution of TDDMgives the correlated ground state and its small ampli-tude limit (STDDM) is an extension of the random-phaseapproximation (RPA) where the effects of ground-statecorrelations and the coupling to two-body amplitudesare included. TDDM and STDDM have been appliedto model hamiltonians [3, 4] and realistic ones [4–7]. Insome cases TDDM overestimates ground-state correla-tions [3] and shows instabilities of the obtained solutions[7] for strongly interacting cases. Obviously the prob-lems of TDDM originate in the truncation scheme of theBBGKY hierarchy where the three-body density matrixis approximated by the one-body and two-body densitymatrices and true three-body correlations are completelyneglected. One way of overcoming the problems is toinclude the equation of motion for the three-body den-sity matrix approximating a four-body density matrixwith lower-level density matrices. Such an attempt hasbeen made for a model hamiltonian [8] and showed thatthe results are drastically improved by including genuinethree-body correlations. However, it is impractical totreat explicitly the equation of motion for the three-bodydensity matrix in realistic cases. In this paper we showthat an approximation for the three-body density matrix with the antisymmetrized products of the two-body den-sity matrices can simulate well the time-evolution of thethree-body density matrix without explicitly solving itsequation of motion. The paper is organized as follows;the formulation and the model hamiltonians are given insect. 2, the results are presented in sect. 3 and sect. 4 isdevoted to a summary. II. FORMULATIONA. Equations of motion for reduced densitymatrices
We consider a system of N fermions and assume thatthe hamiltonian H consisting of a one-body part and atwo-body interaction H = X α ǫ α a † α a α + 12 X αβα ′ β ′ h αβ | v | α ′ β ′ i a † α a † β a β ′ a α ′ , (1)where a † α and a α are the creation and annihilation oper-ators of a particle at a single-particle state α .TDDM gives the coupled equations of motion for theone-body density matrix (the occupation matrix) n αα ′ and the two-body density matrix ρ αβα ′ β ′ . These matricesare defined as n αα ′ ( t ) = h Φ( t ) | a † α ′ a α | Φ( t ) i , (2) ρ αβα ′ β ′ ( t ) = h Φ( t ) | a † α ′ a † β ′ a β a α | Φ( t ) i , (3)where | Φ( t ) i is the time-dependent total wavefunction | Φ( t ) i = exp[ − iHt/ ~ ] | Φ( t = 0) i . The equations inTDDM are written as i ~ ˙ n αα ′ = ( ǫ α − ǫ α ′ ) n αα ′ + X λ λ λ [ h αλ | v | λ λ i ρ λ λ α ′ λ − ρ αλ λ λ h λ λ | v | α ′ λ i ] , (4) i ~ ˙ ρ αβα ′ β ′ = ( ǫ α + ǫ β − ǫ α ′ − ǫ β ′ ) ρ αβα ′ β ′ + X λ λ [ h αβ | v | λ λ i ρ λ λ α ′ β ′ − h λ λ | v | α ′ β ′ i ρ αβλ λ ]+ X λ λ λ [ h αλ | v | λ λ i ρ λ λ βα ′ λ β ′ + h λ β | v | λ λ i ρ λ λ αα ′ λ β ′ − h λ λ | v | α ′ λ i ρ αλ βλ λ β ′ − h λ λ | v | λ β ′ i ρ αλ βλ λ α ′ ] , (5)where ρ αβγα ′ β ′ γ ′ is a three-body density-matrix. InTDDM the BBGKY hierarchy is truncated by approx-imating the three-body density matrix with the antisym-metrized products of the one-body and two-body densitymatrices [1, 2]. The TDDM equations eqs. (4) and (5)conserve the total number of particles N = P α n αα andthe total energy [1, 2] E tot = X α ǫ α n αα + 12 X αβα ′ β ′ h αβ | v | α ′ β ′ i ρ α ′ β ′ αβ . (6) B. Approximation for the three-body densitymatrix
To separate the mean-field contributions, it is conve-nient to introduce the two-body and three-body corre-lation matrices C αβα ′ β ′ and C αβγα ′ β ′ γ ′ which describegenuine two-body and three-body correlations such that ρ αβα ′ β ′ = A ( n αα ′ n ββ ′ ) + C αβα ′ β ′ , (7) ρ αβγα ′ β ′ γ ′ = A ( n αα ′ n ββ ′ n γγ ′ + n αα ′ C βγβ ′ γ ′ )+ C αβγα ′ β ′ γ ′ , (8)where A ( ·· ) means that the quantities in the parenthesisare appropriately antisymmetrized. In TDDM the three-body correlation matrix C αβγα ′ β ′ γ ′ is neglected.First we consider a perturbative expression for thethree-body correlation matrix using the following groundstate | Z i [9] of coupled cluster theory | Z i = e Z | HF i (9)with Z = 14 X pp ′ hh ′ z pp ′ hh ′ a † p a † p ′ a h ′ a h , (10)where | HF i is the Hartree-Fock (HF) ground state and z pp ′ hh ′ is antisymmetric under the exchanges of p ↔ p ′ and h ↔ h ′ . Here, p and h refer to particle and holestates, respectively. In the lowest order of z pp ′ hh ′ thetwo-body correlation matrices are given by C pp ′ hh ′ ≈ z pp ′ hh ′ , (11) C hh ′ pp ′ ≈ z ∗ pp ′ hh ′ , (12) and the three-body correlation matrices by C p p h p p h ≈ X h z ∗ p p hh z p p h h , (13) C p h h p h h ≈ X p z ∗ p ph h z p ph h . (14)These relations give the perturbative expressions for C αβγα ′ β ′ γ ′ in terms of C pp ′ hh ′ C p p h p p h ≈ X h C hh p p C p p h h , (15) C p h h p h h ≈ X p C h h p p C p ph h . (16)Similarly, in the second order of z pp ′ hh ′ the occupationmatrix and other elements of the two-body correlationmatrices are given by n hh ≈ − X pp ′ h ′ C pp ′ hh ′ C hh ′ pp ′ , (17) n pp ≈ X hh ′ p ′ C pp ′ hh ′ C hh ′ pp ′ , (18) C p h p h ≈ X ph C p ph h C hh p p , (19) C p p p p ≈ X hh ′ C p p hh ′ C hh ′ p p , (20) C h h h h ≈ X pp ′ C h h pp ′ C pp ′ h h . (21)Equations (15)-(16) and (19)-(21) imply that the three-body correlation matrix is of the same order of magni-tude as some of the two-body correlation matrix. Fromthe above perturbative considerations we propose a trun-cation scheme of the BBGKY hierarchy that in stead ofneglecting the three-body correlation matrix we includethem using eqs. (15) and (16) in the equation of motionfor the two-body density matrix eq. (5).There are identities which are satisfied by exact re-duced density matrices: n αα ′ = 1 N − X λ ρ αλα ′ λ (22) ρ αβα ′ β ′ = 1 N − X λ ρ αβλα ′ β ′ λ (23) ρ αβγα ′ β ′ γ ′ = 1 N − X λ ρ αβγλα ′ β ′ γ ′ λ (24)and so on, where ρ αβγλα ′ β ′ γ ′ λ is a four-body density ma-trix. Using the correlation matrices, the first three iden-tities are explicitly written as [10] n αα ′ = X λ ( n αλ n λα ′ − C αλα ′ λ ) , (25) C αβα ′ β ′ = 12 X λ ( n αλ C λβα ′ β ′ + n βλ C αλα ′ β ′ + n λα ′ C αβλβ ′ + n λβ ′ C αβα ′ λ − C αβλα ′ β ′ λ ) , (26) C αβγα ′ β ′ γ ′ = 13 X λ ( n αλ C λβγα ′ β ′ γ ′ + n βλ C αλγα ′ β ′ γ ′ + n γλ C αβλα ′ β ′ γ ′ + n λα ′ C αβγλβ ′ γ ′ + n λβ ′ C αβγα ′ λγ ′ + n λγ ′ C αβγα ′ β ′ λ − C αβα ′ λ C γλβ ′ γ ′ − C αβγ ′ λ C γλα ′ β ′ − C αβλβ ′ C γλα ′ γ ′ − C αγα ′ λ C λββ ′ γ ′ − C αγβ ′ λ C βλα ′ γ ′ − C αγλγ ′ C βλα ′ β ′ − C αλα ′ β ′ C βγγ ′ λ − C αλα ′ γ ′ C βγλβ ′ − C αλβ ′ γ ′ C βγα ′ λ − C αβγλα ′ β ′ γ ′ λ ) , (27)where C αβγλα ′ β ′ γ ′ λ is the four-body correlation matrix.The TDDM equations without the three-body correlationmatrix do not conserve the identity eq. (25) [2]. Supposethat the dominant two-body correlation matrices are ei-ther two particle - two hole and two hole - two particletypes and that the deviation of n αα ′ from the HF valuesare small. Then eq. (27) without C αβγλα ′ β ′ γ ′ λ is reducedto eqs. (15) and (16). Similarly, it can be shown usingeqs. (15) and (16) that the occupation matrices and thetwo-body correlation matrices given by eqs. (17)-(21)satisfy eqs. (25) and (26). We also test the followingtruncation scheme: Neglecting the four-body correlationmatrix C αβγδα ′ β ′ γ ′ δ ′ in eq. (27), we use the three-bodycorrelation matrix given by eq. (27) in the equation ofmotion for the two-body density matrix. For simplicitywe do not consider the off-diagonal elements of n αα ′ inthe applications shown below. We show that eqs. (15)and (16) give better results than eq. (27) in the caseof the Lipkin model, indicating that the use of the frac-tional occupation in eq. (27) sometimes overestimatesthe effects of the three-body correlation matrix when thefour-body correlation matrix is neglected.The ground state in TDDM is given as a stationarysolution of the TDDM equations (eqs. (4) and (5)). Weuse the following adiabatic method to obtain a nearlystationary solution [11]: Starting from a non-interactingconfiguration, we solve eqs. (4) and (5) gradually in-creasing the interaction v ( r ) × t/T . To suppress oscillat-ing components which come from the mixing of excitedstates, we must take large T . This method is based onthe Gell-Mann-Low adiabatic theorem [12] and has oftenbeen used to obtain correlated ground states [5, 6]. C. Model hamiltonians
We apply our truncation schemes to the following threecases where comparison with the exact solutions can be made.
1. Lipkin mdel
The Lipkin model [13] describes an N -fermions systemwith two N -fold degenerate levels with energies ǫ/ − ǫ/
2, respectively. The upper and lower levels are la-beled by quantum number p and − p , respectively, with p = 1 , , ..., N . We consider the standard hamiltonianˆ H = ǫ ˆ J z + V J + ˆ J − ) , (28)where the operators are given asˆ J z = 12 N X p =1 ( a † p a p − a − p † a − p ) , (29)ˆ J + = ˆ J †− = N X p =1 a † p a − p . (30)
2. Hubbard model
To test our truncation schemes for a case which in-volves more single-particle states than the Likin model,we consider the one-dimensional (1D) Hubbard modelwith periodic boundary conditions. In momentum spacethe hamiltonian is given by H = X k ,σ ǫ k a † k ,σ a k ,σ + U N X k , p , q ,σ a † k ,σ a k + q ,σ a † p , − σ a p − q , − σ , (31)where U is the on-site Coulomb matrix element, σ spinprojection and the single-particle energies are given by ǫ k = − t P Dd =1 cos( k d ) with the nearest-neighbor hop-ping potential t . We consider the case of the six sites athalf filling. In the first Brillouin zone − π ≤ k < π thereare the following wave numbers k = 0 , k = π , k = − π ,k = 2 π , k = − π . k = − π. (32)The single-particle energies are ǫ = − t , ǫ = ǫ = − t , ǫ = ǫ = t and ǫ = 2 t .
3. Dipolar fermion gas
As a more realistic case, we consider a magnetic dipolargas of fermions with spin one half, which is trapped in aspherically symmetric harmonic potential with frequency ω . The system is described by the hamiltonian H = X α ǫ α a † α a α + 12 X αβα ′ β ′ h αβ | v | α ′ β ′ i a † α a † β a β ′ a α ′ , (33) E t o t / FIG. 1. Ground-state energies in TDDM with the three-bodycorrelation matrix given by eq. (27) (dotted line) and byeqs. (15) and (16) (green (gray) solid line) as a function of χ = ( N − | V | /ǫ for N = 4. The dashed line depicts theresults in TDDM where the three-body correlation matrixis neglected, and the dot-dashed line the exact values. Theresults with the three-body correlation matrix obtained fromsolving the equation of motion are given by the black solidline. where a † α and a α are the creation and annihilation op-erators of an atom at a harmonic oscillator state α cor-responding to the trapping potential V ( r ) = mω r / ǫ α = ω ( n + 3 /
2) with n = 0 , , , .... . We as-sume that α contains the spin quantum number σ . In eq.(33) h αβ | v | α ′ β ′ i is the matrix element of a pure magneticdipole-dipole interaction [14] v ( r ) = − r (3( d · ˆ r )( d · ˆ r ) − d · d ) − π d · d δ ( r ) , (34)where d is the magnetic dipole moment, r = r − r and ˆ r = r /r . The magnetic dipole moment for spin1/2 is given by d = d σ where σ is Pauli matrix. In thecase of completely polarized gases the second term on theright-hand side of eq. (34) can be neglected because theexchange term cancels out the direct term. The contactterm (the second term on the right-hand side of eq. (34))is usually omitted in the study of dipolar gases. However,it is well-known that the contact term for the proton andelectron magnetic dipole moments is essential to explainthe hyperfine splitting of a hydrogen atom. Therefore, inthe following calculations we keep it as it is. We considerthe N = 8 case where noninteracting HF ground stateconsists of the fully occupied 1 s and 1 p harmonic oscil-lator states. For simplicity we take only the 2 s and 1 d states as particle states and keep the 1 s state frozen, i.e.,only the 1 p state is active as a hole-state. n pp FIG. 2. Occupation probability of the upper state as a func-tion of χ for N = 4. The meaning of the five lines is the sameas in fig. 1. C pp ’ - p - p ’ FIG. 3. Two-body correlation matrix C pp ′ − p − p ′ as a functionof χ for N = 4. The meaning of the five lines is the same asin fig. 1. III. RESULTSA. Lipkin model
The ground-state energies obtained from various ap-proximations for the three-body correlation matrix areshown in fig. 1 as a function of χ = ( N − | V | /ǫ for N = 4. The RPA solution becomes unstable at χ = 1.The dashed line shows the results in TDDM where thethree-body correlation matrix is neglected. The dottedand green (gray) solid lines are for the results with thethree-body correlation matrix given by eq. (27) and byeqs. (15) and (16), respectively. The results with thethree-body correlation matrix obtained from solving theequation of motion are given by the black solid line [8].The exact solutions are given by the dot-dashed line. The E t o t / t U/t
FIG. 4. Ground-state energy in TDDM (solid line) with C αβγα ′ β ′ γ ′ given by eqs. (15) and (16) as a function of U/t forthe half-filling six-site Hubbard model. The dashed line de-picts the results in TDDM where the three-body correlationmatrix is neglected, and the squares the exact values. occupation probability of the upper state and the two-body correlation matrix C pp ′ − p − p ′ are shown in figs. 2and 3, respectively. It is clear from figs. 1-3 that the ne-glect of the three-body correlation matrix overestimatesthe ground-state correlations. The perturbative treat-ment eqs. (15) and (16) seems to give better resultsthan the approximation eq. (27) which overly suppressesground-state correlations. Comparing with the resultsobtained with the HF values of n αα , we found that theuse of the fractional occupation n αα in eq. (27) causesthis over suppression. We also found that the omissionof C php ′ h ′ in eq. (5) does not bring serious overestima-tion of the ground-state correlations even if the three-body correlation matrix is excluded. This means that thethree-body correlation matrix plays an important role insuppressing the particle - hole correlations. This is truefor the other two cases of the model hamiltonians. B. Hubbard model
The total energy calculated in TDDM with the three-body correlation matrix given by eqs. (15) and (16)(solid line) is shown in fig. 4 as a function of
U/t forthe six-site Hubbard model with half-filling. The re-sults obtained with eq. (27) are similar and not shown.The RPA solution becomes unstable at
U/t = 2 . C αβα ′ β ′ , we show the correlation energy given by E cor = P αβα ′ β ′ h αβ | v | α ′ β ′ i C α ′ β ′ αβ / n h U/t
FIG. 5. Occupation probability of the initially occupied statewith k as a function of U/t . The meaning of the lines andsymbol is the same as in fig. 4. n h U/t
FIG. 6. Occupation probability of the second occupied statewith k as a function of U/t . The meaning of the lines andsymbol is the same as in fig. 4. lation matrix overestimate the ground-state correlationsand that the approximation given by eqs. (15) and (16)improves the results. The differences among the threecalculations are not so large in the ground-state energyas in the correlation energy. This is due to the fact thatthe increase in the correlation energy is compensated bythat in the mean-field energy.
C. Dipolar atomic gas
The total energy calculated in TDDM with the three-body correlation matrix given by eqs. (15) and (16)is shown by the solid line in fig. 10 as a function of C = d / ~ ωξ for the N = 8 dipolar fermion gas. Here ξ is the oscillator length given by ξ = p ~ /mω . The dashed n p U/t
FIG. 7. Occupation probability of the first unoccupied statewith k as a function of U/t . The meaning of the lines andsymbol is the same as in fig. 4. n p U/t
FIG. 8. Occupation probability of the last unoccupied statewith k as a function of U/t . The meaning of the lines andsymbol is the same as in fig. 4. line shows the results in TDDM without the three-bodycorrelation matrix. The exact solutions are given by thesquares. The RPA solution becomes unstable for a spin-orbit coupling mode at C = 0 .
34 [16]. Due to the spin-orbit coupling in the tensor force, the magnetic quantumnumber m of orbital angular momentum and σ are notgood quantum numbers, and n αα ′ has nonvanishing off-diagonal elements. In the application of eq. (27) we ne-glected the off-diagonal elements because they are small.The results with eq. (27) are similar to those with eqs.(15) and (16) and not shown here. The occupation prob-abilities of the 1 p state with m = 1 and σ = − / d state with m = 2 and σ = − / E cor isshown in fig. 13. Figures 10-13 again show that the ne-glect of the three-body correlation matrix overestimatesthe ground-state correlations and that the approximation U/t E c o r / t FIG. 9. Correlation energy E cor as a function of U/t . Themeaning of the lines and symbol is the same as in fig. 4. E t o t / C FIG. 10. Ground-state energy in TDDM (solid line) with C αβγα ′ β ′ γ ′ given by eqs. (15) and (16) as a function of C forthe N = 8 dipolar fermion gas. The dashed line depicts theresults in TDDM where the three-body correlation matrix isneglected, and the squares the exact values. eqs. (15) and (16) drastically improves the results of theoccupation probabilities and the correlation energy. D. Stability of the ground state solution
In the above we demonstrated that the inclusion of thethree-body correlation matrix improves the ground-stateproperties in TDDM. We now show that it also stabi-lizes long-time behavior of the ground state. To showthe stability of the ground-state solution in TDDM, wepresent the sum S of the absolute value of eq. (25) S = P αα ′ | n αα ′ − P λ ( n αλ n λα ′ − C αλα ′ λ ) | in fig. 14for the Hubbard model. In this model n αα ′ has no off-diagonal elements. The time step 4000 corresponds to T and for t > T the interaction strength is fixed at n h C FIG. 11. Occupation probability of the 1 p state with m = 1and σ = − / C . The meaning of the linesand symbol is the same as in fig. 10. n p C FIG. 12. Occupation probability of the 1 d state with m = 2and σ = − / C . The meaning of the linesand symbol is the same as in fig. 10. U/t = 3. The dotted and solid lines show the resultsin TDDM with the three-body correlation matrix givenby eq. (27) and by eqs. (15) and (16) while the resultswithout the three-body correlation matrix are given bythe dashed line. Since the ground-state solution obtainedin the adiabatic approach inevitably contains small time-dependent components, it has time evolution. If the gra-dient method ref. [4] is used to obtain a true stationarysolution of the TDDM equations, the occupation matrixand the two-body correlation matrix have no time evo-lution and S stays constant though the identity eq. (25)is violated. The violation of the identity eq. (25) causessuch a serious problem for strong interactions that n αα exceeds unity or becomes negative. When n αα becomeunphysical, the TDDM solution is no longer stable asshown in fig. 14. The three-body correlation matrixplays a role in approximately conserving the identity eq. E c o r / C FIG. 13. Correlation energy E cor as a function of C . Themeaning of the lines and symbol is the same as in fig. 10. S Time steps
FIG. 14. Sum of the absolute value of eq.(25) as a functionof time steps at
U/t = 3 for the six-site Hubbard model. Thedotted line and the green line show the results in TDDM withthe three-body correlation matrix given by eq. (27) and byeqs. (15) and (16) while the results without the three-bodycorrelation matrix are given by the dashed line. (25) and stabilizing the time evolution. We found thatneglect of the two-body correlation matrices of three par-ticle - one hole and three hole - one particle types suchas C pp ′ p ′′ h and C hh ′ h ′′ p also stabilizes the time evolutionin the Hubbard model even if the three-body correlationmatrix is neglected. IV. SUMMARY
We proposed truncation schemes of the time-dependent density-matrix approach in which the three-body density matrix is approximated in terms of thesquares of the two-body density matrices. We appliedthem to the ground states of the Lipkin model, the Hub-bard model and the trapped dipolar fermion gas and com-pared with the exact solutions. It was shown that thetruncation schemes give better results than the approachwithout the three-body density matrix. It was pointedout that the three-body correlation matrix plays a rolein suppressing the particle - hole correlations. It was alsoshown that such approximations for the three-body den- sity matrix also drastically improve the stability of theground-state solutions because of approximate conserva-tion of the identity for the one-body and two-body den-sity matrices. Thus it was found that our density-matrixapproach supplemented by the truncation schemes forthe three-body density matrix gives good approximationfor the total ground-state wavefunction. TDDM is nowon good grounds for realistic approximations. [1] S. J. Wang and W. Cassing, Ann. Phys. , 328 (1985).[2] M. Gong and M. Tohyama, Z. Phys. A , 153 (1990).[3] S. Takahara, M. Tohyama and P. Schuck, Phys. Rev.C , 057307 (2004).[4] M. Tohyama, Phys. Rev. C , 044310 (2007).[5] A. Pfitzner, W. Cassing, and A. Peter, Nucl. Phys. A ,753 (1994).[6] M. Assi´e and D. Lacroix, Phys. Rev. Lett. , 202501(2009).[7] M. Tohyama, J. Phys. Soc. Jpn. , 054707 (2012).[8] M. Tohyama and P. Schuck, Eur. Phys. J. A , 257(2010). [9] M. Jema¨i and P. Schuck, Atomic Nuclei , N0. 8, 1139(2011).[10] M. Tohyama, P. Schuck and S. J. Wang, Z. Phys. A ,341 (1991).[11] M. Tohyama, Phys. Rev. A (2005) 043613.[12] M. Gell-Mann and F. Low, Phys. Rev. , 350 (1951).[13] H. J. Lipkin, N. Meshkov and A. J. Glick, Nucl. Phys. , 188 (1965).[14] V. D. Barger and M. G. Olsson, Classical electricity andmagnetism (Allyn and Bacon, Boston, 1987).[15] M. Jema¨i, P. Schuck, J. Dekelsky and R. Bennaceur,Phys. Rev. B , 085115 (2005).[16] M. Tohyama, J. Phys. Soc. Jpn.82