# Microscopic Theory of Multipole Ordering inf-Electron Systems

Authors: Takashi Hotta
Journal: Physics Research International
 (2012)

---

## Abstract

A microscopic framework to determine multipole ordering inf-electron systems is provided on the basis of the standard quantum field theory. For the construction of the framework, a seven-orbital Hubbard Hamiltonian with strong spin-orbit coupling is adopted as a prototype model. A type of multipole and ordering vector is determined from the divergence of multipole susceptibility, which is evaluated in a random phase approximation. As an example of the application of the present framework, a multipole phase diagram on a three-dimensional simple cubic lattice is discussed for the case of n=2, where n denotes the average f-electron number per site. Finally, future problems concerning multipole ordering and fluctuations are briefly discussed.


---

## Body

## 1. Introduction

Recently, complex magnetism in rare-earth and actinide compounds has attracted much attention in the research field of condensed matter physics [1–3]. Since in general, spin-orbit coupling between electrons in 4f and 5f orbitals is strong, spin and orbital degrees of freedom are tightly coupled in f-electron materials. Thus, when we attempt to discuss magnetic ordering in f-electron systems, it is necessary to consider the ordering of spin-orbital complex degrees of freedom, that is, multipole. In fact, ordering of higher-rank multipole has been actively investigated both from experimental and theoretical sides in the research field of strongly correlated f-electron systems [2, 3]. Moreover, due to recent remarkable developments in experimental techniques and measurements, nowadays it has been possible to detect directly and/or indirectly the multipole ordering. Note, however, that only spin degree of freedom often remains, when orbital degeneracy is lifted, for instance, due to the effect of crystal structure with low symmetry. In order to promote the research of multipole phenomena, f-electron compounds crystallizing in the cubic structure with high symmetry are quite important. For instance, octupole ordering has been discussed in the phase IV of Ce0.7La0.3B6 [4] and NpO2 [3, 5–8] with cubic structure. As for NpO2, recently, a possibility of dotriacontapole ordering has been also pointed out [9, 10].Here we emphasize that the study of multipole phenomena has been activated due to the focusing research of filled skutterudite compoundsLnT4X12 with lanthanide Ln, transition metal atom T, and pnictogen X [11]. Since these compounds crystallize in the cubic structure of Th point group, they have provided us an ideal stage for the research of multipole physics. Furthermore, many isostructural materials with different kinds of rare-earth and actinide ions have been successfully synthesized, leading to the development of systematic research on multipole ordering. In fact, recent experiments in close cooperation with phenomenological theory have revealed that multipole ordering frequently appears in filled skutterudites. For instance, a rich phase diagram of PrOs4Sb12 with field-induced quadrupole order has been unveiled experimentally and theoretically [12–14]. Furthermore, antiferro-Γ1-type higher multipole order [2] has been discussed for PrRu4P12 [15, 16] and PrFe4P12 [17–19].Now we turn our attention to theoretical research on multipole order. Thus far, theory of multipole ordering has been developed mainly from a phenomenological viewpoint on the basis of an LS coupling scheme for multi-f-electron state. It is true that several experimental results have been explained by those theoretical studies, but we believe that it is also important to promote microscopic approach for understanding of multipole phenomena in parallel with phenomenological research. Based on this belief, the present author has developed a microscopic theory for multipole-related phenomena with the use of a j-j coupling scheme [1, 20–22]. In particular, octupole ordering in NpO2 has been clarified by the evaluation of multipole interaction with the use of the standard perturbation method in terms of electron hopping [6–8, 23]. We have also discussed possible multipole states of filled skutterudites by analyzing multipole susceptibility of a multiorbital Anderson model based on the j-j coupling scheme [24–29].On the other hand, it is still difficult to understand intuitively the physical meaning of multipole degree of freedom due to the mathematically complicated form of multipole operator defined by using total angular momentum. As mentioned above, multipole is considered to be spin-orbital complex degree of freedom. In this sense, it seems to be natural to regard multipole as anisotropic spin-charge density. This point has been emphasized in the visualization of multipole order [6–8, 23]. Then, we have defined multipole as spin-charge density in the form of one-body operator from the viewpoint of multipole expansion of electromagnetic potential from charge distribution in electromagnetism [30, 31]. Due to the definition of multipole in the form of one-electron spin-charge density operator, it has been possible to discuss unambiguously multipole state by evaluating multipole susceptibility even for heavy rare-earth compounds with large total angular momentum [30].As for the determination of the multipole state, we have proposed to use the optimization of multipole susceptibility on the basis of the standard linear response theory. We have analyzed an impurity Anderson model including sevenf orbitals with the use of the numerical renormalization group technique and checked the effectiveness of the microscopic model on the basis of the j-j coupling scheme for the description of multipoles. We have also shown the result for multipole susceptibility of several kinds of filled skutterudite compounds. With the use of the seven-orbital Anderson model, we have discussed field-induced multipole phenomena in Sm-based filled skutterudites, [32] multipole Kondo effect, [33] and multipole state of Yb- and Tm-based filled skutterudites [34]. We have also discussed possible multipole state in transuranium systems such as AmO2 [35] and magnetic behavior of CmO2 [36].From our previous investigations on the basis of the multiorbital Anderson model, it has been clarified that the multipole can be treated as spin-orbital complex degree of freedom in the one-electron operator form. However, in order to discuss the ordering of multipole, it is necessary to consider a periodic system including sevenf orbitals per atomic site with strong spin-orbit coupling. The validity of the model on the basis of the j-j coupling scheme can be also checked by such consideration. Namely, for the steady promotion of multipole physics, it is highly expected to treat the multipole ordering in a seven-orbital periodic model by overcoming a heavy task to solve the model including 14 states per atomic site.In this paper, we define a seven-orbital Hubbard model with strong spin-orbit coupling and explain a procedure to define the multipole ordering by the divergence of multipole susceptibility from a microscopic viewpoint. For the evaluation of multipole susceptibility, we introduce a random phase approximation. In principle, we can treat all the cases forn=1~13 on the same footing, but here we focus on the case of n=2 corresponding to Pr and U compounds. As a typical example of the present procedure, we show a phase diagram including quadrupole ordering in a three-dimensional simple cubic lattice. Finally, we also discuss some future problems such as superconductivity induced by multipole fluctuations near the multipole phase.The organization of this paper is as follows. In Section2, we explain each part of the seven-orbital Hubbard model with strong spin-orbit coupling. For the reference of readers, we show the list of hopping integrals among f-orbitals along x, y, and z-axes through σ, π, δ, and ϕ bonds. In Section 3, we define the multipole operator as the complex spin-charge degree of freedom in the one-electron form. Then, we explain a scheme to determine the multipole ordering from the multipole susceptibility. Here we use a random phase approximation for the evaluation of the multipole susceptibility. In Section 4, we show the results for the case of n=2 in a three-dimensional simple cubic lattice. We discuss the phase diagram of the multipole ordering. In Section 5, we discuss some future problems and summarize this paper. Throughout this paper, we use such units as ℏ=kB=1.


## 2.  Model Hamiltonian

The model HamiltonianH is split into two parts as(1)H=Hkin+Hloc⁡,
where Hkin denotes a kinetic term and Hloc is a local part for potential and interaction. The latter term is further given by(2)Hloc⁡=Hso+HCEF+HC,
where Hso is a spin-orbit coupling term, HCEF indicates crystalline electric field (CEF) potential term, and HC denotes Coulomb interaction term. We explain each term in the following.


### 2.1.  Localf-Electron Term

Among the three terms ofHloc, the spin-orbit coupling part is given by(3)Hso=λ∑i,m,σ,m′,σ′ζm,σ;m′,σ′fimσ†fim′σ′,
where fimσ is an annihilation operator of f-electron at site i, σ=+1 (−1) for up (down) spin, m is the z-component of angular momentum ℓ=3, λ is the spin-orbit interaction, and the matrix elements are expressed by(4)ζm,σ;m,σ=mσ2,ζm+σ,-σ;m,σ=l(l+1)-m(m+σ)2,
and zero for other cases.Next we consider the CEF term, which is expressed as(5)HCEF=∑i,m,m′,σBm,m′fimσ†fim′σ,
where Bm,m′ is the CEF potential for f electrons from the ligand ions, which is determined from the table of Hutchings for angular momentum ℓ=3 [37]. For the cubic structure with Oh symmetry, Bm,m′ is expressed by using three CEF parameters, B40 and B60, as(6)B3,3=B-3,-3=180B40+180B60,B2,2=B-2,-2=-420B40-1080B60,B1,1=B-1,-1=60B40+2700B60,B0,0=360B40-3600B60,B3,-1=B-3,1=6015(B40-21B60),B2,-2=300B40+7560B60.
Note the relation of Bm,m′=Bm′,m. Following the traditional notation [38], we define
(7)B40=WxF(4),B60=W(1-|x|)F(6),
where W determines an energy scale for the CEF potential, x specifies the CEF scheme for Oh point group, and F(4)=15 and F(6)=180 for ℓ=3.Finally, the Coulomb interaction termHC is given by(8)H1=∑i,m1~m4∑σ,σ′Im1m2,m3m4fim1σ†fim2σ′†fim3σ′fim4σ,
where the Coulomb integral Im1m2,m3m4 is expressed by(9)Im1m2,m3m4=∑k=06Fkck(m1,m4)ck(m2,m3).
Here Fk is the Slater-Condon parameter and ck is the Gaunt coefficient which is tabulated in the standard textbooks of quantum mechanics [39]. Note that the sum is limited by the Wigner-Eckart theorem to k=0, 2, 4, and 6. The Slater-Condon parameters should be determined for the material from the experimental results, but in this paper, for a purely theoretical purpose, we set the ratio among the Slater-Condon parameters as physically reasonable values, given by(10)F0:F2:F4:F6=10:5:3:1.
Note that F6 is considered to indicate the scale of Hund’s rule interaction among f orbitals.


### 2.2.  Kinetic Term

Next we consider the kinetic term off electrons. When we discuss magnetic properties of f-electron materials as well as the formation of heavy quasiparticles, it is necessary to include simultaneously both conduction electrons with wide bandwidth and f electrons with narrow bandwidth, since the hybridization is essentially important for the formation of heavy quasiparticles. In this sense, it is more realistic to construct orbital-degenerate periodic Anderson model for the theory of multipole ordering in heavy-electron systems.However, if we set the starting point of the discussion in the periodic Anderson model, the calculations for multipole susceptibility will be very complicated. Thus, we determine our mind to split the problem into two steps: namely, first we treat the formation of heavy quasiparticles and then, we discuss the effective model for such heavy quasiparticles. If we correctly include the symmetry off-electron orbital, we believe that it is possible to grasp qualitatively correct points concerning the multipole ordering by using an effective kinetic term for f electrons.Based on the above belief, we consider the effective kinetic term in a tight-binding approximation forf electrons. Then, Hkin is expressed as(11)Hkin=∑i,a,m,m′,σtm,m′afimσ†fi+am′σ,
where tm,m′a indicates the f-electron hopping between m- and m′-orbitals of adjacent atoms along the a direction. The hopping amplitudes are obtained from the table of Slater-Koster integrals, [40–42] but, for convenience, here we show explicitly tm,m′a on the three-dimensional cubic lattice.The hopping integrals along thez-axis are given in quite simple forms as(12)t0,0z=(ffσ),t-1,-1z=t1,1z=(ffπ),t-2,-2z=t2,2z=(ffδ),t-3,-3z=t3,3z=(ffϕ),
and zeros for other cases. Here (ffℓ) denotes the Slater-Koster integral through ℓ bond between nearest neighbor sites. Note that the above equations are closely related to the definitions of (ffσ), (ffπ), (ffδ), and (ffϕ).On the other hand, hopping integrals along thex- and y-axes are given by the linear combination of (ffσ), (ffπ), (ffδ), and (ffϕ). We express tm,m′a as(13)tm,m′a=∑l(ffl)Em,m′al,
where the coefficient Em,m′aℓ indicates the two-center integral along a direction between m and m′ orbitals and ℓ runs among σ, π, δ, and ϕ. In Table 1, we show the values of Em,m′aℓ. Other components are zeros unless they are obtained with the use of relation of Em,m′aℓ=Em′,maℓ=E-m,-m′aℓ.Table 1
CoefficientsEm,m′aℓ along the x- and y-axes between f orbitals of nearest neighbor sites. Note that in double signs, the upper and lower signs correspond to the value along the x- and y-axes, respectively.
mm′σπδϕ−3−35/1615/323/161/32−3−1∓15/16∓15/32±15/16±15/32−3115/16-15/32-15/1615/32−33∓5/16±15/32∓3/16±1/32−2−205/161/23/16−200∓30/160±30/16−2205/16-1/23/16−1−13/161/325/165/32−11∓3/16±1/32∓5/16±5/320003/805/8By using the experimental results concerning the Fermi-surface sheets for actual materials, it is possible to determine the Slater-Koster parameters,(ffσ), (ffπ), (ffδ), and (ffϕ), so as to reproduce the experimental results. Namely, the hopping integrals should be effective ones for quasiparticles, as mentioned above. Here it is important to include correctly the symmetry of local f orbitals in the evaluation of hopping amplitudes, although the whole energy scale will be adjusted by experimental results and band-structure calculations.


## 2.1.  Localf-Electron Term

Among the three terms ofHloc, the spin-orbit coupling part is given by(3)Hso=λ∑i,m,σ,m′,σ′ζm,σ;m′,σ′fimσ†fim′σ′,
where fimσ is an annihilation operator of f-electron at site i, σ=+1 (−1) for up (down) spin, m is the z-component of angular momentum ℓ=3, λ is the spin-orbit interaction, and the matrix elements are expressed by(4)ζm,σ;m,σ=mσ2,ζm+σ,-σ;m,σ=l(l+1)-m(m+σ)2,
and zero for other cases.Next we consider the CEF term, which is expressed as(5)HCEF=∑i,m,m′,σBm,m′fimσ†fim′σ,
where Bm,m′ is the CEF potential for f electrons from the ligand ions, which is determined from the table of Hutchings for angular momentum ℓ=3 [37]. For the cubic structure with Oh symmetry, Bm,m′ is expressed by using three CEF parameters, B40 and B60, as(6)B3,3=B-3,-3=180B40+180B60,B2,2=B-2,-2=-420B40-1080B60,B1,1=B-1,-1=60B40+2700B60,B0,0=360B40-3600B60,B3,-1=B-3,1=6015(B40-21B60),B2,-2=300B40+7560B60.
Note the relation of Bm,m′=Bm′,m. Following the traditional notation [38], we define
(7)B40=WxF(4),B60=W(1-|x|)F(6),
where W determines an energy scale for the CEF potential, x specifies the CEF scheme for Oh point group, and F(4)=15 and F(6)=180 for ℓ=3.Finally, the Coulomb interaction termHC is given by(8)H1=∑i,m1~m4∑σ,σ′Im1m2,m3m4fim1σ†fim2σ′†fim3σ′fim4σ,
where the Coulomb integral Im1m2,m3m4 is expressed by(9)Im1m2,m3m4=∑k=06Fkck(m1,m4)ck(m2,m3).
Here Fk is the Slater-Condon parameter and ck is the Gaunt coefficient which is tabulated in the standard textbooks of quantum mechanics [39]. Note that the sum is limited by the Wigner-Eckart theorem to k=0, 2, 4, and 6. The Slater-Condon parameters should be determined for the material from the experimental results, but in this paper, for a purely theoretical purpose, we set the ratio among the Slater-Condon parameters as physically reasonable values, given by(10)F0:F2:F4:F6=10:5:3:1.
Note that F6 is considered to indicate the scale of Hund’s rule interaction among f orbitals.


## 2.2.  Kinetic Term

Next we consider the kinetic term off electrons. When we discuss magnetic properties of f-electron materials as well as the formation of heavy quasiparticles, it is necessary to include simultaneously both conduction electrons with wide bandwidth and f electrons with narrow bandwidth, since the hybridization is essentially important for the formation of heavy quasiparticles. In this sense, it is more realistic to construct orbital-degenerate periodic Anderson model for the theory of multipole ordering in heavy-electron systems.However, if we set the starting point of the discussion in the periodic Anderson model, the calculations for multipole susceptibility will be very complicated. Thus, we determine our mind to split the problem into two steps: namely, first we treat the formation of heavy quasiparticles and then, we discuss the effective model for such heavy quasiparticles. If we correctly include the symmetry off-electron orbital, we believe that it is possible to grasp qualitatively correct points concerning the multipole ordering by using an effective kinetic term for f electrons.Based on the above belief, we consider the effective kinetic term in a tight-binding approximation forf electrons. Then, Hkin is expressed as(11)Hkin=∑i,a,m,m′,σtm,m′afimσ†fi+am′σ,
where tm,m′a indicates the f-electron hopping between m- and m′-orbitals of adjacent atoms along the a direction. The hopping amplitudes are obtained from the table of Slater-Koster integrals, [40–42] but, for convenience, here we show explicitly tm,m′a on the three-dimensional cubic lattice.The hopping integrals along thez-axis are given in quite simple forms as(12)t0,0z=(ffσ),t-1,-1z=t1,1z=(ffπ),t-2,-2z=t2,2z=(ffδ),t-3,-3z=t3,3z=(ffϕ),
and zeros for other cases. Here (ffℓ) denotes the Slater-Koster integral through ℓ bond between nearest neighbor sites. Note that the above equations are closely related to the definitions of (ffσ), (ffπ), (ffδ), and (ffϕ).On the other hand, hopping integrals along thex- and y-axes are given by the linear combination of (ffσ), (ffπ), (ffδ), and (ffϕ). We express tm,m′a as(13)tm,m′a=∑l(ffl)Em,m′al,
where the coefficient Em,m′aℓ indicates the two-center integral along a direction between m and m′ orbitals and ℓ runs among σ, π, δ, and ϕ. In Table 1, we show the values of Em,m′aℓ. Other components are zeros unless they are obtained with the use of relation of Em,m′aℓ=Em′,maℓ=E-m,-m′aℓ.Table 1
CoefficientsEm,m′aℓ along the x- and y-axes between f orbitals of nearest neighbor sites. Note that in double signs, the upper and lower signs correspond to the value along the x- and y-axes, respectively.
mm′σπδϕ−3−35/1615/323/161/32−3−1∓15/16∓15/32±15/16±15/32−3115/16-15/32-15/1615/32−33∓5/16±15/32∓3/16±1/32−2−205/161/23/16−200∓30/160±30/16−2205/16-1/23/16−1−13/161/325/165/32−11∓3/16±1/32∓5/16±5/320003/805/8By using the experimental results concerning the Fermi-surface sheets for actual materials, it is possible to determine the Slater-Koster parameters,(ffσ), (ffπ), (ffδ), and (ffϕ), so as to reproduce the experimental results. Namely, the hopping integrals should be effective ones for quasiparticles, as mentioned above. Here it is important to include correctly the symmetry of local f orbitals in the evaluation of hopping amplitudes, although the whole energy scale will be adjusted by experimental results and band-structure calculations.


## 3.  Multipole Ordering

In order to discuss the multipole ordered phase from the itinerant side, we evaluate the multipole susceptibilityχ by following the standard quantum field theory. The multipole susceptibility is defined by(14)χ(q,iνn)=∫01/Tdτeiνnτ〈X̂q(τ)X̂-q†(0)〉,
where X̂q denotes the multipole operator with momentum q, ν=2πTn is the boson Matsubara frequency with an integer n, T is a temperature, Xq(τ)=eHτXqe-Hτ, and 〈⋯〉 indicates the thermal average by using H. In the following, we introduce the multipole operator and explain a method to evaluate the susceptibility.


### 3.1.  Multipole Operator

In any case, first it is necessary to define multipole. As for the definition of multipole, readers should consult with [30, 31], but here we briefly explain the definition in order to make this paper self-contained. We define X in the one-electron density-operator form as(15)X̂q=∑k,γpk,γ(q)T̂γ(k)(q),
where k denotes the rank of multipole, γ indicates the irreducible representation for cubic point group, and T̂γ(k)(q) indicates the cubic tensor operator, expressed in the second-quantized form as(16)T̂γ(k)(q)=∑k,m,σ,m′γ′Tmσ,m′σ′(k,γ)fkmσ†fk+qm′σ′.
Here the matrix elements of the coefficient T̂(k,γ) are calculated from the Wigner-Eckert theorem as [43]
(17)Tmσ,m′σ′(k,γ)=∑j,μ,μ′,qGγ,q(k)〈j‖T(k)‖j〉2j+1〈jμ∣jμ′kq〉×〈jμ∣lmsσ2〉〈jμ′∣lm′sσ′2〉,
where ℓ=3, s=1/2, j=ℓ±s, μ runs between -j and j, q runs between -k and k, Gγ,q(k) is the transformation matrix between spherical and cubic harmonics, 〈JM|J^′M′J′′M′′〉 denotes the Clebsch-Gordan coefficient, and 〈j∥T(k)∥j〉 is the reduced matrix element for spherical tensor operator, given by(18)〈j‖T(k)‖j〉=12k(2j+k+1)!(2j-k)!.
Note that k≤2j and the highest rank is 2j. When we define multipoles as tensor operators in the space of total angular momentum J on the basis of the LS coupling scheme, there appear multipoles with k≥8 for the cases of J≥4, that is, for 2≤n≤4 and 8≤n≤12, where n is local f-electron number. If we need such higher-rank multipoles with k≥8, it is necessary to consider many-body operators beyond the present one-body definition.Note that when we express the multipole moment as (16) and (17), we normalize each multipole operator so as to satisfy the orthonormal condition [44](19)Tr⁡{T̂(k,γ)T̂(k′,γ′)}=δkk′δγγ′,
where δkk′ denotes the Kronecker’s delta.


### 3.2.  Multipole Susceptibility

Now we move to the evaluation of multipole susceptibility. In order to determine the coefficientpk,γ(q) in (15), it is necessary to calculate the multipole susceptibility in the linear response theory. The multipole susceptibility is expressed as(20)χ(q,iνn)=∑kγ,k′γ′pk,γχkγ,k′γ′(q,iνn)pk′,γ′*,
where the susceptibility matrix is given by(21)χkγ,k′γ′(q,iνn)=∑m1~m4∑σ1~σ4Tm1σ1,m3σ3(k,γ)×χm1σ1m2σ2,m3σ3m4σ4(q,iνn)Tm2σ2,m4σ4(k,γ)*.
Then, χ and pk,γ are determined by the maximum eigenvalue and the corresponding normalized eigenstate of the susceptibility matrix equation (21).In order to calculate actually the multipole susceptibility, it is necessary to introduce an appropriate approximation. In this paper, we use a random phase approximation (RPA) for the evaluation of multipole susceptibility. For the purpose, we redivide the HamiltonianH into two parts as(22)H=H0+H1,
where H0 indicates the one-electron part given by H0=Hkin+Hso+HCEF and H1 is the interaction part, which is just equal to HC in the present case. Then, we consider the perturbation expansion in terms of the Coulomb interaction.The susceptibility diagrams are shown in Figure1 and they are expressed in a compact matrix form as
(23)χ̂=χ̂(0)[1̂-Ûχ̂(0)]-1+χ̂(0)[1̂+Ĵχ̂(0)]-1-χ̂(0),
where Û and Ĵ are, respectively, given by(24)Um1σ1m2σ2,m3σ3m4σ4=Im1m2,m3m4δσ1σ4δσ2σ3,Jm1σ1m2σ2,m3σ3m4σ4=Im1m2,m4m3δσ1σ3δσ2σ4,
and the dynamical susceptibility χ̂(0) is given by(25)χm1σ1m2σ2,m3σ3m4σ4(0)(q,iνn)=-T∑n′∑kGm1σ1,m4σ4(0)(k,iωn′)×Gm2σ2,m3σ3(0)(k+q,iωn′+iνn).
Here G(0) is the one-electron Green’s function defined by the noninteracting part H0.Figure 1
Feynman diagrams for multipole susceptibility in the RPA. The solid curve and broken line denote the noninteracting Green’s functionG(0) and Coulomb interaction, respectively.In order to determine the multipole ordering, it is necessary to detect the divergence ofχ at νn=0. We cannot evaluate the susceptibility just at a diverging point, but we find such a critical point by the extrapolation of 1/χmax as a function of U, where U indicates the energy scale of the Slater-Condon parameters and χmax denotes the maximum eigenvalue of susceptibility matrix equation (21) for νn=0. When we increase the magnitude of U, 1/χmax is gradually decreased from the value in the weak-coupling limit. In actual calculations, we terminate the calculation when 1/χmax arrives at a value in the order of unity. By using the calculated values of 1/χmax, we make an extrapolation of 1/χmax as a function of U. Then, we find a critical value of U at which 1/χmax becomes zero. As for the type of multipole and ordering vector in the ordered phase, we extract such information from the eignevectors of the susceptibility matrix corresponding to the maximum eigenvalue. By performing the above calculations, it is possible to find the multipole ordered phase from a microscopic viewpoint in principle.


## 3.1.  Multipole Operator

In any case, first it is necessary to define multipole. As for the definition of multipole, readers should consult with [30, 31], but here we briefly explain the definition in order to make this paper self-contained. We define X in the one-electron density-operator form as(15)X̂q=∑k,γpk,γ(q)T̂γ(k)(q),
where k denotes the rank of multipole, γ indicates the irreducible representation for cubic point group, and T̂γ(k)(q) indicates the cubic tensor operator, expressed in the second-quantized form as(16)T̂γ(k)(q)=∑k,m,σ,m′γ′Tmσ,m′σ′(k,γ)fkmσ†fk+qm′σ′.
Here the matrix elements of the coefficient T̂(k,γ) are calculated from the Wigner-Eckert theorem as [43]
(17)Tmσ,m′σ′(k,γ)=∑j,μ,μ′,qGγ,q(k)〈j‖T(k)‖j〉2j+1〈jμ∣jμ′kq〉×〈jμ∣lmsσ2〉〈jμ′∣lm′sσ′2〉,
where ℓ=3, s=1/2, j=ℓ±s, μ runs between -j and j, q runs between -k and k, Gγ,q(k) is the transformation matrix between spherical and cubic harmonics, 〈JM|J^′M′J′′M′′〉 denotes the Clebsch-Gordan coefficient, and 〈j∥T(k)∥j〉 is the reduced matrix element for spherical tensor operator, given by(18)〈j‖T(k)‖j〉=12k(2j+k+1)!(2j-k)!.
Note that k≤2j and the highest rank is 2j. When we define multipoles as tensor operators in the space of total angular momentum J on the basis of the LS coupling scheme, there appear multipoles with k≥8 for the cases of J≥4, that is, for 2≤n≤4 and 8≤n≤12, where n is local f-electron number. If we need such higher-rank multipoles with k≥8, it is necessary to consider many-body operators beyond the present one-body definition.Note that when we express the multipole moment as (16) and (17), we normalize each multipole operator so as to satisfy the orthonormal condition [44](19)Tr⁡{T̂(k,γ)T̂(k′,γ′)}=δkk′δγγ′,
where δkk′ denotes the Kronecker’s delta.


## 3.2.  Multipole Susceptibility

Now we move to the evaluation of multipole susceptibility. In order to determine the coefficientpk,γ(q) in (15), it is necessary to calculate the multipole susceptibility in the linear response theory. The multipole susceptibility is expressed as(20)χ(q,iνn)=∑kγ,k′γ′pk,γχkγ,k′γ′(q,iνn)pk′,γ′*,
where the susceptibility matrix is given by(21)χkγ,k′γ′(q,iνn)=∑m1~m4∑σ1~σ4Tm1σ1,m3σ3(k,γ)×χm1σ1m2σ2,m3σ3m4σ4(q,iνn)Tm2σ2,m4σ4(k,γ)*.
Then, χ and pk,γ are determined by the maximum eigenvalue and the corresponding normalized eigenstate of the susceptibility matrix equation (21).In order to calculate actually the multipole susceptibility, it is necessary to introduce an appropriate approximation. In this paper, we use a random phase approximation (RPA) for the evaluation of multipole susceptibility. For the purpose, we redivide the HamiltonianH into two parts as(22)H=H0+H1,
where H0 indicates the one-electron part given by H0=Hkin+Hso+HCEF and H1 is the interaction part, which is just equal to HC in the present case. Then, we consider the perturbation expansion in terms of the Coulomb interaction.The susceptibility diagrams are shown in Figure1 and they are expressed in a compact matrix form as
(23)χ̂=χ̂(0)[1̂-Ûχ̂(0)]-1+χ̂(0)[1̂+Ĵχ̂(0)]-1-χ̂(0),
where Û and Ĵ are, respectively, given by(24)Um1σ1m2σ2,m3σ3m4σ4=Im1m2,m3m4δσ1σ4δσ2σ3,Jm1σ1m2σ2,m3σ3m4σ4=Im1m2,m4m3δσ1σ3δσ2σ4,
and the dynamical susceptibility χ̂(0) is given by(25)χm1σ1m2σ2,m3σ3m4σ4(0)(q,iνn)=-T∑n′∑kGm1σ1,m4σ4(0)(k,iωn′)×Gm2σ2,m3σ3(0)(k+q,iωn′+iνn).
Here G(0) is the one-electron Green’s function defined by the noninteracting part H0.Figure 1
Feynman diagrams for multipole susceptibility in the RPA. The solid curve and broken line denote the noninteracting Green’s functionG(0) and Coulomb interaction, respectively.In order to determine the multipole ordering, it is necessary to detect the divergence ofχ at νn=0. We cannot evaluate the susceptibility just at a diverging point, but we find such a critical point by the extrapolation of 1/χmax as a function of U, where U indicates the energy scale of the Slater-Condon parameters and χmax denotes the maximum eigenvalue of susceptibility matrix equation (21) for νn=0. When we increase the magnitude of U, 1/χmax is gradually decreased from the value in the weak-coupling limit. In actual calculations, we terminate the calculation when 1/χmax arrives at a value in the order of unity. By using the calculated values of 1/χmax, we make an extrapolation of 1/χmax as a function of U. Then, we find a critical value of U at which 1/χmax becomes zero. As for the type of multipole and ordering vector in the ordered phase, we extract such information from the eignevectors of the susceptibility matrix corresponding to the maximum eigenvalue. By performing the above calculations, it is possible to find the multipole ordered phase from a microscopic viewpoint in principle.


## 4.  Results

In the previous sections, we have explained the model Hamiltonian and the procedure to determine the type of multipole ordering. We believe that the present procedure can be applied to actual materials, but there are so many kinds of materials and multipole phenomena. Here we show the calculated results for the case ofn=2 concerning Γ3 non-Kramers quadrupole ordering, in order to see how the present procedure works. The results for actual materials will be discussed elsewhere.


### 4.1.  CEF States

First we discuss the local CEF states in order to determine the CEF parameter. We consider the case ofn=2 corresponding to Pr3+ and U4+ ions. Since we discuss the local electron state, the energy unit is taken as F6. As for the spin-orbit coupling, here we take λ/F6=0.1. Concerning the value of W, it should be smaller than λ and we set W as W/F6=0.001.In Figure2, we show the CEF energies as functions of x. As easily understood from the discussion in the LS coupling scheme, the ground state multiplet for n=2 is characterized by J=4, where J is total angular momentum given by J=|L-S| with angular momentum L and spin momentum S. For n=2, we find L=5 and S=1 from the Hund’s rules and, thus, we obtain J=4. Due to the effect of cubic CEF, the nonet of J=4 is split into four groups as Γ1 singlet, Γ3 non-Kramers doublet, Γ4 triplet, and Γ5 triplet. In the present diagonalization of Hloc, we find such CEF states, as shown in Figure 2. When we compare this CEF energy diagram with that of the LS coupling scheme [38], we find that the shape of curves and the magnitude of excitation energy are different with each other. However, from the viewpoint of symmetry, the structure of the low-energy states is not changed between the LS and j-j coupling schemes [1]. Since we are interested in a possibility of Γ3 quadrupole ordering, we choose the value of x as x=0.0 in the following.Figure 2
CEF energy levels obtained by the diagonalization ofHloc for λ/F6=0.1 and W/F6=0.001 with F0=10F6, F2=5F6, and F4=3F6.


### 4.2.  Energy Bands

Next we consider the band structure obtained by the diagonalization ofH0=Hkin+HCEF+Hso. As for the Slater-Koster integrals, it is one way to determine them so as to reproduce the Fermi-surface sheets of actual materials, but here we determine them from a theoretical viewpoint as(26)-(ffσ)=(ffδ)=t,(ffπ)=(ffϕ)=-t/2,
where t indicates the magnitude of hopping amplitude. The size of t should be determined by the quasi-particle bandwidth, but here we simply treat it as an energy unit.In Figure3, we depict the eigen energies of H0 along the lines connecting some symmetric points in the first Brillouin zone. As for the spin-orbit coupling and CEF parameters, we set λ/t=0.1 and W/t=0.001. First we note that there exist seven bands and each band has double degeneracy due to time-reversal symmetry, which is distinguished by pseudospin. Since the magnitude of λ is not so large, we do not observe a clear splitting between j=7/2 octet and j=5/2 sextet bands. Around at Γ point, we find that j=5/2 sextet is split into two groups, Γ7 doublet and Γ8 quartet. Here we note that the energy of Γ8 quartet is lower than that of Γ7. Since the Γ8 has orbital degeneracy, it becomes an origin of the formation of Γ3 non-Kramers doublet, when we accommodate a couple of electrons per site.Figure 3
Energy band structure obtained by the diagonalization ofH0 for (ffσ)=-t, (ffπ)=-t/2, (ffδ)=t, (ffϕ)=-t/2, λ/t=0.1, and W/t=0.001. Note that we show the eigen energies along the lines of Γ→X→M→Γ→R→X.Note that the Fermi level is denoted by a horizontal line, which is determined by the condition ofn=2, where n is the average electron number per site. When we pay our attention to the band near the Fermi level, we find that the orbital degeneracy exists in the bands on the Fermi surface. For instance, we see the degenerate bands on the Fermi surface around the Γ point. Such orbital degeneracy in the momentum space is considered to be a possible source of Γ3 quadrupole ordering, which will be discussed in the next subsection. Finally, in the present case, we expect the appearance of the large-volume Fermi surface as well as the small-size pocket-like Fermi surface. Such mixture of the Fermi surface sheets with different topology may be an important issue for the appearance of higher-rank multipole ordering.


### 4.3.  Phase Diagram

Now we show the phase diagram of the multipole state. First it is necessary to calculate the susceptibility equation (25) at νn=0. As for the momentum q, we divide the first Brillouin zone into 16×16×16 meshes. Concerning the momentum integration in (25), we exploit the Gauss-Legendre quadrature with due care. At low temperatures such as T/t=0.01, it seems to be enough to divide the range between -π and π into 60 segments along each direction axis. As found in (25), χ(0) has 144 components in the spin-orbital space, but it is not necessary to calculate all the components due to the symmetry argument. We have checked that it is enough to evaluate 1586 components of χ(0).We set the parameters asλ/t=0.1, x=0.0, W/t=0.001, F2=0.5F0, F4=0.3F0, F6=0.1F0, -(ffσ)=(ffδ)=t, and (ffπ)=(ffϕ)=-t/2. Note that the ratio among the Slater-Condon parameters is the same as that in Figure 2. We also note that the hopping amplitude t is relatively large compared with local potential and interactions, since we consider the multipole ordering from the itinerant side. Here we emphasize that our framework actually works for the microscopic discussion on the multipole ordering. A way to determine more realistic parameters in the model will be discussed elsewhere.By changing the values of temperatureT/t, we depict the phase diagram in the plane of t/F0 and T/t. Note that t2/F0 corresponds to the typical magnitude of multipole-multipole interaction between nearest neighbor sites. As naively understood, when the temperature is increased, larger value of U is needed to obtain the ordered state. Then, the phase diagram is shown in Figure 4. We evaluate the maximum eigenvalue of the multipole susceptibility by increasing F0/t. One may think that the magnitude of t/F0 in Figure 4 is too small to obtain reasonable results in the RPA calculations. Here we note that the total bandwidth of the seven-orbital system is in the order of 10t, as shown in Figure 3. Namely, the critical value of the interaction Fc0 at low enough temperatures is considered to be in the order of the total bandwidth. In this sense, we consider that the value of t/F0 in Figure 4 is not small for the RPA calculations. Note also that when the temperature is increased, the magnitude of noninteracting susceptibility is totally suppressed, leading to the enhancement of Fc0. Thus, t/F0 is decreased when T is increased, as observed in Figure 4.Figure 4
Phase diagram of the multipole ordering forn=2 on the three-dimensional simple cubic lattice.At low temperatures asT/t<0.3, we obtain that the maximum eigen value of susceptibility matrix is characterized by the multipole with Γ3 symmetry and the ordering vector Q=(π,π,π). The component of the multipole depends on the temperature, but the 90% of the optimized multipole is rank 2 (quadrupole). Others are rank 4 (hexadecapole) and rank 6 (tetrahexacontapole) components, which are about 10%. Note again that the multipoles with the same symmetry are mixed in general, even if the rank of the multipole is different. Namely, quadrupole is the main component, while hexadecapole and tetrahexacontapole are included with significant amounts. Note also that the phase diagram is shown only in the region of T/t<1, but the boundary curve approaches the line of t/F0=0. Since the case with very large F0 is unrealistic, we do not pay our attention to the phase for T>t, although we can continue the calculation in such higher temperature region.When we increase the temperature, the magnetic phase is observed forT/t>0.3. The main component is Γ4 dipole and the ordering vector is Q=(0,0,0). Note that the susceptibility for Γ4 multipole moment does not mean magnetic susceptibility, which is evaluated by the response of magnetic moment L+2S, that is, J+S. At T/t=0.4, admixture of the multipole is as follows: rank 1 (dipole) 90.7%, rank 3 (octupole) 6.5%, rank 5 (dotriacontapole) 2.1%, and rank 7 (octacosahectapole) 0.7%. The amounts are changed by the temperature, but the main component is always dipole. We have found the low-temperature antiferroquadrupole state and the high-temperature ferromagnetic phase. Such a combination of nonmagnetic and magnetic phases can be observed in other parameter sets including quadrupole ordering.


## 4.1.  CEF States

First we discuss the local CEF states in order to determine the CEF parameter. We consider the case ofn=2 corresponding to Pr3+ and U4+ ions. Since we discuss the local electron state, the energy unit is taken as F6. As for the spin-orbit coupling, here we take λ/F6=0.1. Concerning the value of W, it should be smaller than λ and we set W as W/F6=0.001.In Figure2, we show the CEF energies as functions of x. As easily understood from the discussion in the LS coupling scheme, the ground state multiplet for n=2 is characterized by J=4, where J is total angular momentum given by J=|L-S| with angular momentum L and spin momentum S. For n=2, we find L=5 and S=1 from the Hund’s rules and, thus, we obtain J=4. Due to the effect of cubic CEF, the nonet of J=4 is split into four groups as Γ1 singlet, Γ3 non-Kramers doublet, Γ4 triplet, and Γ5 triplet. In the present diagonalization of Hloc, we find such CEF states, as shown in Figure 2. When we compare this CEF energy diagram with that of the LS coupling scheme [38], we find that the shape of curves and the magnitude of excitation energy are different with each other. However, from the viewpoint of symmetry, the structure of the low-energy states is not changed between the LS and j-j coupling schemes [1]. Since we are interested in a possibility of Γ3 quadrupole ordering, we choose the value of x as x=0.0 in the following.Figure 2
CEF energy levels obtained by the diagonalization ofHloc for λ/F6=0.1 and W/F6=0.001 with F0=10F6, F2=5F6, and F4=3F6.


## 4.2.  Energy Bands

Next we consider the band structure obtained by the diagonalization ofH0=Hkin+HCEF+Hso. As for the Slater-Koster integrals, it is one way to determine them so as to reproduce the Fermi-surface sheets of actual materials, but here we determine them from a theoretical viewpoint as(26)-(ffσ)=(ffδ)=t,(ffπ)=(ffϕ)=-t/2,
where t indicates the magnitude of hopping amplitude. The size of t should be determined by the quasi-particle bandwidth, but here we simply treat it as an energy unit.In Figure3, we depict the eigen energies of H0 along the lines connecting some symmetric points in the first Brillouin zone. As for the spin-orbit coupling and CEF parameters, we set λ/t=0.1 and W/t=0.001. First we note that there exist seven bands and each band has double degeneracy due to time-reversal symmetry, which is distinguished by pseudospin. Since the magnitude of λ is not so large, we do not observe a clear splitting between j=7/2 octet and j=5/2 sextet bands. Around at Γ point, we find that j=5/2 sextet is split into two groups, Γ7 doublet and Γ8 quartet. Here we note that the energy of Γ8 quartet is lower than that of Γ7. Since the Γ8 has orbital degeneracy, it becomes an origin of the formation of Γ3 non-Kramers doublet, when we accommodate a couple of electrons per site.Figure 3
Energy band structure obtained by the diagonalization ofH0 for (ffσ)=-t, (ffπ)=-t/2, (ffδ)=t, (ffϕ)=-t/2, λ/t=0.1, and W/t=0.001. Note that we show the eigen energies along the lines of Γ→X→M→Γ→R→X.Note that the Fermi level is denoted by a horizontal line, which is determined by the condition ofn=2, where n is the average electron number per site. When we pay our attention to the band near the Fermi level, we find that the orbital degeneracy exists in the bands on the Fermi surface. For instance, we see the degenerate bands on the Fermi surface around the Γ point. Such orbital degeneracy in the momentum space is considered to be a possible source of Γ3 quadrupole ordering, which will be discussed in the next subsection. Finally, in the present case, we expect the appearance of the large-volume Fermi surface as well as the small-size pocket-like Fermi surface. Such mixture of the Fermi surface sheets with different topology may be an important issue for the appearance of higher-rank multipole ordering.


## 4.3.  Phase Diagram

Now we show the phase diagram of the multipole state. First it is necessary to calculate the susceptibility equation (25) at νn=0. As for the momentum q, we divide the first Brillouin zone into 16×16×16 meshes. Concerning the momentum integration in (25), we exploit the Gauss-Legendre quadrature with due care. At low temperatures such as T/t=0.01, it seems to be enough to divide the range between -π and π into 60 segments along each direction axis. As found in (25), χ(0) has 144 components in the spin-orbital space, but it is not necessary to calculate all the components due to the symmetry argument. We have checked that it is enough to evaluate 1586 components of χ(0).We set the parameters asλ/t=0.1, x=0.0, W/t=0.001, F2=0.5F0, F4=0.3F0, F6=0.1F0, -(ffσ)=(ffδ)=t, and (ffπ)=(ffϕ)=-t/2. Note that the ratio among the Slater-Condon parameters is the same as that in Figure 2. We also note that the hopping amplitude t is relatively large compared with local potential and interactions, since we consider the multipole ordering from the itinerant side. Here we emphasize that our framework actually works for the microscopic discussion on the multipole ordering. A way to determine more realistic parameters in the model will be discussed elsewhere.By changing the values of temperatureT/t, we depict the phase diagram in the plane of t/F0 and T/t. Note that t2/F0 corresponds to the typical magnitude of multipole-multipole interaction between nearest neighbor sites. As naively understood, when the temperature is increased, larger value of U is needed to obtain the ordered state. Then, the phase diagram is shown in Figure 4. We evaluate the maximum eigenvalue of the multipole susceptibility by increasing F0/t. One may think that the magnitude of t/F0 in Figure 4 is too small to obtain reasonable results in the RPA calculations. Here we note that the total bandwidth of the seven-orbital system is in the order of 10t, as shown in Figure 3. Namely, the critical value of the interaction Fc0 at low enough temperatures is considered to be in the order of the total bandwidth. In this sense, we consider that the value of t/F0 in Figure 4 is not small for the RPA calculations. Note also that when the temperature is increased, the magnitude of noninteracting susceptibility is totally suppressed, leading to the enhancement of Fc0. Thus, t/F0 is decreased when T is increased, as observed in Figure 4.Figure 4
Phase diagram of the multipole ordering forn=2 on the three-dimensional simple cubic lattice.At low temperatures asT/t<0.3, we obtain that the maximum eigen value of susceptibility matrix is characterized by the multipole with Γ3 symmetry and the ordering vector Q=(π,π,π). The component of the multipole depends on the temperature, but the 90% of the optimized multipole is rank 2 (quadrupole). Others are rank 4 (hexadecapole) and rank 6 (tetrahexacontapole) components, which are about 10%. Note again that the multipoles with the same symmetry are mixed in general, even if the rank of the multipole is different. Namely, quadrupole is the main component, while hexadecapole and tetrahexacontapole are included with significant amounts. Note also that the phase diagram is shown only in the region of T/t<1, but the boundary curve approaches the line of t/F0=0. Since the case with very large F0 is unrealistic, we do not pay our attention to the phase for T>t, although we can continue the calculation in such higher temperature region.When we increase the temperature, the magnetic phase is observed forT/t>0.3. The main component is Γ4 dipole and the ordering vector is Q=(0,0,0). Note that the susceptibility for Γ4 multipole moment does not mean magnetic susceptibility, which is evaluated by the response of magnetic moment L+2S, that is, J+S. At T/t=0.4, admixture of the multipole is as follows: rank 1 (dipole) 90.7%, rank 3 (octupole) 6.5%, rank 5 (dotriacontapole) 2.1%, and rank 7 (octacosahectapole) 0.7%. The amounts are changed by the temperature, but the main component is always dipole. We have found the low-temperature antiferroquadrupole state and the high-temperature ferromagnetic phase. Such a combination of nonmagnetic and magnetic phases can be observed in other parameter sets including quadrupole ordering.


## 5.  Discussion and Summary

We have constructed the microscopic framework to discuss the multipole ordering due to the evaluation of multipole susceptibility inf-electron systems on the basis of the seven-orbital Hubbard model with strong spin-orbit coupling. For the evaluation of multipole susceptibility, we have used the RPA and found the critical point from 1/χmax. As an example of the present scheme, we have shown the results for the case of n=2 concerning quadrupole ordering on the three-dimensional simple cubic lattice. If we specify the lattice structure and determine the hopping parameter from the comparison with the experimental results on the Fermi-surface sheets, in principle, it is possible to determine the type of multipole ordering with the use of appropriate local CEF parameters and Coulomb interactions.Although the microscopic theory of multipole ordering has been proposed, it is necessary to elaborate the present scheme both from theoretical and experimental viewpoints. In order to enhance the effectiveness of the present procedure, we should increase the applicability of the theory. For instance, we have not considered at all the sublattice structure in this paper, but in actuality, the staggered-type multipole ordering has been observed. In order to reproduce the structure, it is necessary to maximize the multipole susceptibility by taking into account the sublattice structure. It is one of future problems from a theoretical viewpoint.It is also highly expected that the present scheme should be applied to actual materials in order to explain the origin of multipole ordering. For instance, it is interesting to seek for the origin of peculiar incommensurate quadrupole ordering observed in PrPb3 [45]. At the first glance, it seems to be quite difficult to explain the origin of the Γ3 quadrupole ordering with the ordering vector of Q=(π/2±δ,π/2±δ,0) with δ=π/8. However, if we use the present scheme, it may be possible to find a solution in a systematic way. Another issue is the revisit to octupole and higher-rank multipole ordering in NpO2. The significant amount of dotriacontapole component may be understood naturally in the present scheme.Another interesting future problem is the emergence of superconductivity near the multipole ordered phase. It has been widely accepted that anisotropicd-wave superconductivity appears in the vicinity of the antiferromagnetic phase, as observed in several kinds of strongly correlated electron materials. In general, near the quantum phase transition, anisotropic superconducting pairs are formed due to the effect of quantum critical fluctuations. Thus, also in the vicinity of multipole ordering, superconductivity is generally expected to occur. Even from purely theoretical interest, it is worthwhile to investigate superconductivity near the antiferroquadrupole phase in Figure 4. When we turn our attention to actual material, in PrIr2Zn20, superconductivity has been observed and quadrupole fluctuations have been considered to play some roles [46]. Within the RPA, it is possible to discuss the appearance of superconductivity in the vicinity of quadrupole ordering in the present scheme. It is another future problem.In summary, we have proposed the prescription to determine the type of multipole ordering from a microscopic viewpoint on the basis of the seven-orbital Hubbard model. The multipole susceptibility has been obtained in the RPA and the quadrupole ordering has been actually discussed in a way similar to that for the spin ordering in the single-orbital Hubbard model. The application to actualf-electron materials will be discussed elsewhere, but we believe that the present scheme is useful to consider the origin of multipole ordering. In addition, a possibility of superconductivity near the multipole ordering is an interesting future problem.



---

*Source: 762798-2012-05-13.xml*