Drude Weight of the Two-Dimensional Hubbard Model -- Reexamination of Finite-Size Effect in Exact Diagonalization Study --

0.7 MB·English

by  Hiroki Nakano

Typeset with jpsj2.cls <ver.1.2> Full Paper
Drude Weight of the Two-Dimensional Hubbard Model
– Reexamination of Finite-Size Effect in Exact Diagonalization Study –
Hiroki Nakano∗, Yoshinori Takahashi and Masatoshi Imada1
7 Graduate School of Material Science, University of Hyogo, Kouto 3-2-1, Kamiori, Ako-gun, 678-1297
0 1Department of Applied Physics, University of Tokyo, Hongo 7-3-1, Tokyo 113-8656
0
(ReceivedFebruary6,2008)
2
TheDrudeweightoftheHubbardmodelonthetwo-dimensionalsquarelatticeisstudiedby
n
the exact diagonalizations applied to clusters up to 20 sites. We carefully examine finite-size
a
effects by consideration of the appropriate shapes of clusters and the appropriate boundary
J
condition beyond the limitation of employing only the simple periodic boundary condition.
0
We successfully capture the behavior of the Drude weight that is proportional to the squared
3
holedopingconcentration.Ourpresentresultgivesaconsistentunderstandingofthetransition

betweentheMottinsulatoranddopedmetals.Wealsofind,inthefrequencydependenceofthe
]
l opticalconductivity,thatthemid-gapincoherentpartemergesmorequicklythanthecoherent
e
- part and rather insensitive to the doping concentration in accordance with the scaling of the
r Drudeweight.
t
s
t. KEYWORDS: Hubbardmodel,Lanczosmethod,metal-insulatortransition,Drudeweight,opticalconduc-
a tivity
m
-
d
1. Introduction ical quantity as the Drude weight of the strongly corre-
n
o Metal-insulatortransitionisafundamentalissueinthe lated system without approximations. One possible way
c condensed-matter physics. Particular interest has been is to apply the exact diagonalization method to finite-
[ focused on properties of metals near the transition be- size clusters. The first work along this line concerning
1 tween the Mott insulator and the metallic phase. Many with the two-dimensional Hubbard model was reported
v studies have already been made from both theoretical byDagottoetal.2 Basedonthetreatmentontheclusters
5 and experimental points of view.1 upto16sites,theysuccessfullycapturedaroughfeature
3
Oneofthequantitiesbywhichonecandistinguishbe- oftheDrudeweightasafunctionoftheholedopingcon-
7
1 tween metallic or insulating phases is the Drude weight centration δ. They tentatively concluded that D ∝ δ,
0 D,namelythecoherentcomponentoftheopticalconduc- although the Drude weight sometimes becomes nega-
7 tivityσ(ω).The systemismetallicorinsulatingdepend- tivewithlargeabsolutevaluesespeciallynearthe metal-
0 ing on D 6= 0 or D = 0, respectively. In the Mott insu- insulator transition point. Such negative Drude weights
t/ lator,the Drude weightthereforealwaysvanishes.Ifone occur as a consequence of the finite-size effect. In order
a
m dopestheMottinsulatorwithcarriers,itbecomesmetal- to get rid of such a finite-size effect, Nakano and Imada
lic and the Drude weight increases with further doping. examined the properties under various boundary condi-
d- TheHubbardmodelisbyfarthesimplestmicroscopic tions in the system ofthe same size with 16sites.3 They
n model which can show the Mott-insulator-metal transi- consequently found the concave behavior of the Drude
o tion. It is the most appropriate model for intensive the- weight as a function of δ in the small-δ region, which
c oretical studies. Precise determination of its properties implies the dependence D ∝δ2. Recently, Tohyama and
:
v by some methods with no approximations contributes his co-workers treated larger clusters of 18 and 20 sites
Xi much to the fundamental understanding of the metal- under the periodic boundary condition.4 Their result of
insulatortransition.Inspiteofthesimplicity,variousas- the Drude weight seems to suggest D ∝ δ again. Thus,
r
a pectsoftypicalmany-bodyproblemsareinvolvedinthis the issue of the δ-dependence of D is still controversial
model due to the presence of the Coulomb interaction. at present.
Many of them are still controversial even now. Among Under these circumstances, we examine the finite-size
themis the hole-dopingdependence ofthe Drude weight effectsoftheDrudeweightofclusterswithvariousshapes
near the transition between the Mott insulator and the and sizes up to 20 sites under various boundary condi-
metal. Primary interest of this paper is the dependence tions very carefully and do our best in minimizing the
of the two-dimensional Hubbard model for small doping finite-size effects when we calculate the Drude weight of
concentration; this is because the Mott insulator of this the finite-size clusters. The purpose of the study in this
model at half filling is considered to become metallic in paper is to capture thereby a reliable exponent for the
anontrivialwaywhenthe doping is switchedonandbe- δ-dependence of the Drude weight near the transition
causethismodelhasbeenextensivelystudiedasamodel between the metal and the Mott insulator. Our present
of high-T cuprates. results suggest that D ∝δ2.
c
Generally,it is notso easy to calculate sucha dynam- Thispaperisorganizedasfollows.Inthenextsection,
themodelHamiltonianandthemethodofnumericalcal-
culation are introduced. Section 3 is devoted to the ex-
∗E-mailaddress:[email protected]
2 J.Phys.Soc.Jpn. FullPaper H.Nakano,Y.Takahashi andM.Imada
amination of the finite-size effect in the non-interacting
case. In Section 4, results in the interacting case are re-
portedanddiscussed.Summaryandconclusionaregiven
in the final section.
2. Model and Method
We consider the two-dimensional Hubbard model. Its
Hamiltonian is given by
H = H +H , (1)
hop int
H = −t c† c ,
hop i,σ j,σ
X
hi,ji,σ
Hint = U ni↑ni↓, Fig. 1. Finite size clusters for Ns = 16, 18 and 20. The blue
Xi and green vectors represent two translational vectors for each
cluster. The atomic sites are denoted by the red circles inside
where the creation (annihilation) of an electron at site
thequadrangleformedbythesetwovectorsandtheotheredges
i with spin σ is denoted by c†i,σ (ci,σ) and ni,σ is the denotedbyblacksolidlines.
corresponding number operator, namely n = c† c .
i,σ i,σ i,σ
In this paper, we treat only the nearest-neighbor hop-
ping on the two-dimensional square lattice. Shapes for deviation from the value in the thermodynamic limit.
finite-size clusters will be illustrated later. Energies are In ref. 2, a remarkable finite-size effect was reported. In
measured in units of t; therefore, we set t=1 hereafter. view of this situation, effects of other types of boundary
The system consists of N electrons on N atomic sites. conditions, namely the antiperiodic boundary condition
e s
The hole doping concentration is a controllable param- and the mixed boundary condition, were studied in ref.
eter given by δ = 1−n, where the electron density is 3 to see whether one can reduce the finite-size effects
defined as n=Ne/Ns. for the fixed Ns or not. Consequently, the wrong nega-
We perform numerical diagonalization of finite-size tivelargeDrudeweightsweresuccessfullyresolved.Note
clusters of the model (1) using Lanczos algorithm. This here that when one imposes the mixed boundary condi-
method is non-biased against the effects of Coulomb tion, the two lattice directions are no longer equivalent
interaction. The method is also useful to obtain the and the average of the two directions gives a reasonable
ground-statewavefunctionwiththeground-stateenergy result. This suggests that adhering to the equivalence of
E very precisely. On the other hand, the disadvantage the two directions is not necessary. The two directions
g
of this method is that available system sizes are lim- becomeinequivalentwhenonerelaxestheconditionthat
ited to be small. This is because the dimension of the these quadrangles of the clusters have the right angle at
Hilbert space grows exponentially with respect to the their vertices. Therefore we rather additionally examine
systemsize.5 Inorderto overcomethis disadvantage,we the parallelogramclusters depicted in Fig. 1 (d) and (e)
perform parallel calculations in the computer systems for Ns = 16 and 18 in place of squared-shaped clusters,
of the distributed-memory type. Thereby we can treat respectively.Animportantpoint is thata setof possible
arbitrary shapes of clusters up to N = 20 if we use ap- wave number vectors corresponding to the clusters (d)
s
propriate supercomputer systems. We also employ the and (e) is different from the one corresponding to the
continued-fractionexpansionmethod,6whichprovidesus clusters (a) and (b), respectively. The system is never-
with dynamical quantities including the optical conduc- theless still bipartite. For Ns = 20, on the other hand,
tivityandtheDrudeweight.Thedetailedprocedurewill it is impossible to preserve the bipartite nature when
beexplainedaftertheexaminationsofclustershapesand one distorts the cluster (c) so that it realizes a parallel-
boundary conditions. ogram cluster with a different set of possible wave num-
ber vectors. Then we study only one type of cluster for
3. Finite-Size Effect for U = 0 and Boundary N =20. Therefore we consider these clusters (a)-(e) up
s
Condition to N =20.
s
In this section, we examine how finite-size effects oc- Next, we show how to select the appropriate cluster
cur in various types of boundary conditions of clusters under the appropriate boundary condition for a given
for a given Ns. When one considers the system of two- Ns.Forthis purpose,letusfirstsee the ground-stateen-
dimensional square lattice by using a finite-size cluster, ergy Eg of the clusters (a)-(e) at half filling for U = 0,
regular squares under the periodic boundary condition presented in Table I. Note that this quantity is propor-
are usually adopted as in refs. 2 and 4. Such regular tionalto theDrude weightforU =0.Anoticeablepoint
squaresforN =16,18and20aredepictedinFig.1(a), of the results of clusters (a), (b) and (c) under the pe-
s
(b) and (c), respectively. The merit of adopting these riodic boundary condition is that the relative difference
squares is that the x and y directions are equivalent. (E∞ −Eg)/E∞ does not monotonically converge as Ns
When one treats such a quantity as the Drude weight, is increased, where E∞ denotes the ground-state energy
however,calculationsbasedontheregularsquaresunder in the limit of Ns → ∞. This suggests that the system
the periodic boundaryconditionreveala largefinite-size size Ns up to 20 is not free from the finite-size effect as
J.Phys.Soc.Jpn. FullPaper H.Nakano,Y.TakahashiandM.Imada 3
far as we stick to the periodic boundary condition. We
have thus to be careful in choosing an appropriate clus-
ter under anappropriateboundary conditionfor a given
N .Letusnowcomparethedifference(E −E )/E of
s ∞ g ∞
alltheclusters(a)-(e)withvariousboundaryconditions.
One can find the best cluster and boundary condition
fromamongthe casesshowninTable I sothatthe abso-
lute value of the difference is the smallest. For N = 16
s
and 20, two boundary conditions give the smallest one.
Sets of possible wave number vectors of the two cases
for eachN are equivalent.Thus, the same result will be
s
obtained if one selects either of the two.
Let us compare the δ-dependence of the non-
interacting E between two cases, the one of the most
g
appropriate cluster under the most appropriate bound-
ary condition with the best E at half filling and the
g
other of the regular square under the periodic boundary
condition; results are depicted in Fig. 2. One can easily
observe the results of the former case fall closer to the
curve of the thermodynamic limit in the whole regionof
δ,irrespectiveofthe sizeN .Onthe otherhand,thelat-
s
terregular-squarecaseundertheperiodicboundarycon-
ditiongivessignificantoverestimatesparticularlyat(N ,
s
N ) = (18, 18) and (20, 18). This strongly suggests that
e
theDrudeweightsofthesecaseswillbepossiblyoveresti-
matedevenwhenU isconsiderablylarge.Italsosuggests
that if the intrinsic Drude weight shows the concave be-
havioras a function of δ for large U, this overestimation
will possibly cancelthis behavior.Figure 2 alsoprovides
another important finite-size effect that finite-size data
Table I. List of ground-state energy per site in the non-
interactingcaseU =0athalffillingforvariousshapesofclusters
andvariousboundaryconditionstogetherwithrelativedifference
from the energy in the thermodynamic limit. Head BC1 (BC2)
denotes the boundary condition along the green (blue) transla-
tionalvectorinFig.1.PandArepresenttheperiodicboundary
conditionandtheantiperiodicboundarycondition,respectively.
Noteherethattheboundaryconditions givingthesameenergy
for each cluster correspond to the identical set of possible wave
numbervectors.MarkedcasesinlastcolumnwithheadPWshow Fig. 2. Ground-state energies as functionofδ innon-interacting
the selected clusters under the selected boundary condition in caseforgivenNspresentedintheformof−Eg/(4Ns).Thisquan-
whichweactuallyperformcalculations inthepresentwork. titycorresponds toD forU =0.Greentrianglesrepresentdata
cluster BC1 BC2 −Eg/Ns (E∞−Eg)/E∞ PW for the regular square under the periodic boundary condition.
(a) P P 1.5 −0.0747 Redcirclesrepresentdata fortheappropriatecluster under the
A P 1.70711 +0.0530 appropriate boundary condition based on the analysis in Table
P A 1.70711 +0.0530 I. The solid curve denotes the corresponding quantity for the
A A 1.41421 −0.1276 thermodynamiclimit.
(b) P P 1.77778 +0.0966
A P 1.53960 −0.0503
P A 1.53960 −0.0503
show nonzero curvature in their δ-dependence only at δ
A A 1.33333 −0.1775
correspondingtothe closedelectronicshellstructure.At
(c) P P 1.69443 +0.0452
A P 1.65065 +0.0182 • δ for the open electronic shell structure, on the other
P A 1.65065 +0.0182 hand, Eg follows linear dependence. The linear depen-
A A 1.52169 −0.0613 dence from this origin is likely to survive more or less
(d) P P 1.70711 +0.0530 evenforU >0.Thisis,infact,confirmedinthequantum
A P 1.63099 +0.0061 •
Monte Carlostudy atnonzeroU.7 Thus,it is difficult to
P A 1.5 −0.0747
judgewhetherthelinearδ-dependencerevealedfromnu-
A A 1.63099 +0.0061
(e) P P 1.72417 +0.0636 mericalstudiesintheinteractingcaseintheregionofthe
A P 1.64519 +0.0148 • openshellstructureisintrinsicorissimplyattributedto
P A 1.51621 −0.0647 this particular electronic structure. In view of this situ-
A A 1.44675 −0.1076
ation, it is desirable to confine our analysis to the data
4 J.Phys.Soc.Jpn. FullPaper H.Nakano,Y.Takahashi andM.Imada
at δ for the closed shell structure in order to extract the
intrinsic behavior. In this paper, therefore, we take the
appropriateclustersunderthetheappropriateboundary
conditionforeachN (showninTable I)andadoptdata
s
for the filling forming the closed shell structure.
In the selected clusters and boundary conditions, the
equivalence of the two directions denoted by α (= x,y)
doesnothold,aswehavealreadymentionedintheabove.
Then,wetakethedirectionaverageaccordingtothepro-
cedure used in ref. 3. The optical conductivity is given
by σ(ω) = 2πDδ(ω)+σreg(ω). Here, the second term is
theregularpartofthe opticalconductivity,whichisalso
called the incoherent part. The regular part σreg(ω) is
obtainedby the averageσreg(ω)=[σreg(ω)+σreg(ω)]/2;
x y
each of σreg(ω) given by
α
Fig. 3. Ground-stateenergiesasfunctionofδininteractingcases
σreg(ω)= π |hℓ|jα|0i|2δ(ω−E +E ), (2) presented in the form of −Eg/(4Ns). The solid curve denotes
α Ns X Eℓ−E0 ℓ 0 −Eg/(4Ns) for U =0 for the thermodynamic limit.Closed cir-
ℓ(6=0) cles and closed triangles denote the present results from exact
is calculated by the continued-fraction expansion diagonalization calculations for U = 4 and U = 35, respec-
tively. Open symbols (triangles, squares, reversed triangles and
method. Here, j is a current operator along the α-
α diamonds) denote results of Ns = 6×6, 8×8, 10×10 and
direction defined as 12×12 for U =4 from the quantum Monte Carlo simulations,
j =−i t(c† c −c† c ), (3) respectively.7
α i,σ i+δα,σ i+δα,σ i,σ
X
i,σ
whereδ istheunitvectoralongtheα-direction,|ℓirep- the finite-size effects in D and N based on the present
α eff
resentsaneigenstatewiththeenergyeigenvalueE .Note procedure are also small if the smooth δ-dependence of
ℓ
that the groundstate is describedby ℓ=0.One can ob- D and N is confirmed irrespective of the size N .
eff s
tain the Drude weight D from the combination of σ(ω) In this paper, we take three cases of U = 10, 20 and
and the sum rule 35.Inthe systemwithanonzeroU weakerthanU =10,
∞ itisdifficulttodistinguishtheupper-Hubbardbandand
σ(ω)dω =πK, (4)
Z the lower-Hubbard band in the responses of σ(ω); the
0
estimation of N is not available. On the other hand,
eff
where−4Krepresentsthekineticenergypersite,namely
in the system for U > 35, the ground state is possibly
−4K = hH i/N . For U = 0, the incoherent part is
hop s ferromagnetic near half filling. Once the ferromagnetism
absent;thereforeonehasD =K.Itiswidelyknownthat
appears, it is too difficult to study the scaling of D near
for U >0, the incoherentpart appears.When U is large
half filling.
enough, one can easily recognize in σ(ω) the responses
towardtheupper-Hubbardband.Insuchacase,onecan 4. Results for Interacting Case and Discussions
definitely determine a frequency just below the upper-
Weshownumericalresultsfortheδ-dependenceofN
eff
Hubbardbandtobeω ;wealsocalculatedefiniteintegral
c connected by lines for interacting cases in Fig. 4. The δ-
1 ωc dependenceoftheDrudeweightislatershown.Although
N = σ(ω)dω, (5)
eff π Z the data consistof those of various system sizes for each
0
case of U, δ-dependence of them is fairly smooth, inde-
which is called an effective carrier density.
pendent of the value U, as shown in Fig. 4(a). It clearly
As for the validity of extending the present procedure
indicates that deviations of our numerical results due to
to cases with finite U, there is no convincing theoretical
the finite-size effect are very small for every U that we
argument on the condition of reducing the finite-size ef-
treat. One can also observe that as δ is decreased to-
fect.Thoughnotreliableinastrictandrigoroussense,it
wardthe half-filling caseδ =0,N goesto a verysmall
isatleastbettertostartfromsituationswheretheeffects eff
value for each of non-zero U. The values of N at half
are reduced in the case of U = 0 from the continuity of eff
filling are very close to the corresponding D. This is be-
behavior as the function of U. Another evidence for the
cause no significant responses within 0 < ω < ω are
validityofthepresentprocedureisgivenfromourexam- c
observed in the optical conductivity; this behavior has
inationofE forU >0;theresultsaredepictedinFig.3.
g already been reported in all the previous works2–4 and
Althoughallthe possible data ofN =16,18and 20are
s willbeconfirmedlaterinthepresentstudy.Then,wewill
shown there with the same closed symbol for each case
later come back to these values of N at half filling in
of U >0, δ-dependence of E is fairly smooth. This fact eff
g the discussion of the Drude weight. In order to examine
implies that finite-size deviations ofthe presentdata are
the δ-dependence of N near half filling, Fig. 4(b) de-
small. One can actually confirm that the present results eff
pictstheδ-dependenceofN /δinsteadofrawN .One
for U = 4 agree well with the results from the quan- eff eff
can observe the linear behavior irrespective of the value
tum Monte Carlo simulations7 on larger systems up to
of U in the wide region of δ. They extrapolate to finite
N =12×12.Therefore,itisreasonabletopresumethat
s
J.Phys.Soc.Jpn. FullPaper H.Nakano,Y.TakahashiandM.Imada 5
Fig. 6. Log-log plot of D and Neff as function of δ for U = 35.
Circlesandtrianglesdenotes D andNeff,respectively.
inthe thermodynamic limit.Allthe values forvariousU
are so small, less than 0.003, that we cannot distinguish
their symbols from the origin as shown in Fig. 4(a) and
5, where our data of the half-filling case for N =20 are
s
plotted. For N = 16, we obtain D ∼ 0.009 for U = 10
s
at most; since its difference from the data of N = 20 is
s
onlywithinthesymbolsize,itisnotshowninthefigure.
On the other hand, ref. 4 reported the Drude weight of
about 0.02, which is substantially larger than ours, for
thecluster(c)withN =20undertheperiodicboundary
Fig. 4. δ-dependence ofNeff in(a)andNeff/δin(b).Closedtri- s
condition.Ourpresentresultisclosertothecorrectvalue
angles, open squares and closed circles denote data for U =10,
20 and 35, respectively. In (a), total weight for U = 0 is also ofzeroinspiteofthesmallersizeNs,whichindicatessu-
depicted bythesolidcurveforcomparison. periority of our present calculations in reproducing well
the insulating state at half filling.
In the metallic region for δ >0, all the data obtained
from studies on various system sizes are presented in
Figs. 4(a) and 5 according to the discussion in section
3.Recallhere thatrefs.2and4 reportedanomalousand
discontinuous deviations of D at some particular hole
concentrations, resulting from a finite-size effect.11 On
the other hand, the present data for all the values of δ
in Fig. 5 show smooth variation on the whole, though
they consist of mixture for various N . It indicates that
s
the finite-size effect is well eliminated in our estimates.
Asfortheδ-dependenceoftheDrudeweight,theconvex
behaviorisapparentinalltheregionofδ forU =10.For
U =20 and35,on the other hand, the concavebehavior
occurs near half filling. The change in this δ-dependence
with increasing U will be discussed below.
InordertoconfirmtheconcavebehaviorofD,weplot
Fig. 5. δ-dependence of Drude weight for U = 10, 20 and 35.
the data of D on the logarithmicscale in Fig.6 together
Notations for closed triangles, open squares and closed circles
are the same as those in Fig. 4. For comparison, solid curve with Neff for U = 35. Our data of D near half filling
represents D forU =0,whichisthesameasthetotal weight. indicate the dependence
D ∝δ2. (7)
intercepts of the N /δ axis, resulting in the following This dependence is in contrast to the linear dependence
eff
dependence ofN ofeq.(6).Formoredetailedanalysisofthedepen-
eff
dence of D in the region of small δ, let us next examine
N ∝δ, (6)
eff the δ-dependence of the ratio D/N . The correspond-
eff
near half filling δ = 0. Note that the above dependence ing quantity wasalsostudied inthe two-dimensionalt-J
(6) is in agreement with all the previous studies.2–4 model.12 The present results are shown in Fig. 7 in the
Wedepictδ-dependenceoftheDrudeweightinFig.5. formof(D/N )/δ.ForU =10and20,asδ isdecreased
eff
Note that the Drude weight at δ =0 vanishes for U >0 toward half filling, the quantity (D/N )/δ shows steep
eff
6 J.Phys.Soc.Jpn. FullPaper H.Nakano,Y.Takahashi andM.Imada
Fig. 8. δ-dependenceofD/δ.Notationsforclosedtriangles,open
Fig. 7. δ-dependence ofD/(Neffδ) andits inverse. Notations for
squaresandclosedcirclesarethesameasthoseinFig.4.Broken
closedtriangles,opensquaresandclosedcirclesarethesameas
linesarethefittinglinesin0.25≤δ≤0.5.Dottedlinesareguides
thoseinFig.4.
foreyes intheassumptionofD∝δ2.
increase.ForU =35,ontheotherhand,nosuchbehavior
tem sizes not yet treated by exact-diagonalization stud-
isobserved.TheincreasingbehaviorforU =10and20is
ies. The QMC calculation15 also demonstrated that the
examined in the plot of the inverse of (D/N )/δ shown
eff localizationlengthsuggeststhedivergenceξ ∝|µ−µ |−ν
in the inset of Fig.7. Our data for U =10 and20 reveal l c
with ν = 1/4 as the chemical potential approaches the
linear dependence in a wide region of δ. For U = 35, a
charge gap µ from the insulating side. The scaling the-
similar linear dependence is observed in 0.2 <∼ δ <∼ 0.5 c
ory indicates zν = 1, which also leads to the exponent
togetherwithadeviationfromitinthesmaller-δ region.
z = 4. Therefore, our present result is consistent with
Each of the linear extrapolation of plots intersects at a
these QMC results in views of the scaling theory. Re-
non-zerointerceptwith the axis of ordinatesin the limit
cent study of the metal-insulator transition16–18 clari-
of δ →0. If one takes its weak δ-dependence for U =35
fied that the exponent z = 4 appears at the marginal
intoaccount,itisnaturaltoconsiderthat[(D/N )/δ]−1
eff
quantum critical point between the continuous metal-
forvariousU staysfiniteinthelimitofδ →0ascommon
insulator transition with T = 0 and the discontinuous
behavior; there is no sign of vanishing [(D/N )/δ]−1 c
eff
transition with nonzero T . The appearance of the same
from our data of finite size clusters. The above behav- c
marginal exponent z = 4 was also found in the one-
ior, together with eq. (6), thus leads us to the scaling
dimensional Hubbard model with next-nearest-neighbor
(7) valid not only for U = 35 but also for U = 10 and
hopping.10,19 The unified understandingofthe marginal
20. This analysis is particularly important since it sug-
point occurring in the one- and two-dimensionalcases is
gests the validity of the scaling (7) even for U = 10 in
the subject for future studies. Although it is difficult to
spite of the fact that the δ-dependence of D for this U
determine the marginalquantumcriticalpointprecisely,
in Fig. 5 does not show concave behavior which appears
the presentresults suggestthat the critical regionof the
for U =20 and 35.
marginal point is wide extending to all the values of U
Let us now discuss the present result from the view-
we studied. This allows us to capture successfully the
pointofthescalingtheorydevelopedbyImada,ageneral
distinct signature characteristic to the marginal point.
theory of the metal-insulator transition.13 This helps us
Next, let us consider the width of the critical region
to justify the dependence (7) in relation to other behav-
of δ in which the behavior (7) is observed. For this
iorofthetwo-dimensionalHubbardmodel.Accordingto
purpose, we depict the δ-dependence of D/δ in Fig. 8.
the scalingtheory,the Drude weightshowsthe following
critical behavior near the metal-insulator transition, In 0.25 <∼ δ <∼ 0.5, one can observe U-irrelevant lin-
ear behavior. In Fig. 8, we fit our numerical data in
D ∝δ1+(z−2)/d, (8) 0.25 ≤ δ ≤ 0.5 with broken lines. For U = 20 and 35,
our data near half filling digress from the broken lines.
where z denotes the dynamical exponent and d is the
In particular, for U = 35, two points shows clear de-
spatial dimension. Since d = 2 in the present problem,
viation toward the origin in Fig. 8, while one point for
one obtains z =4 from (7). The scaling theory also pre-
U = 20. Here, let us draw the dotted guide lines for
dictsthecriticalbehaviorofthecompressibilitynearthe
these digressingdataby assumingthe presentresult(7).
transition,
Thedottedandthebrokenlinescrosswitheachotherfor
κ∝δ1−z/d. (9) eachU.OnecanregardthecrossingpointforeachU asa
rough estimation of the boundary of the critical region.
Substitutionofz =4andd=2intheabovedependence
Then, the estimated boundaries are δ ∼ 0.14 and 0.22
of κ leads to κ ∝ δ−1. The diverging behavior of κ was
for U =20 and 35, respectively. Though crude are these
reported by the quantum Monte Carlo (QMC) simula-
estimations, this is the first study to report the bound-
tions7,14 performed on the same model with larger sys-
J.Phys.Soc.Jpn. FullPaper H.Nakano,Y.TakahashiandM.Imada 7
ary. Calculations of systems with larger N with proper
s
treatment of finite-size effects will make the precise es-
timations possible in future studies. For U = 10, there
arenodatapointsshowingdigressingtowardstheorigin.
If the δ2 dependence of D, i.e. eq. (7), is true, then the
above fact suggests that the critical region for U =10 is
narrowerthan δ <1/9.This is consistent with the rapid
decrease of the estimated boundary as U decreasesfrom
35 to 20. The doping concentration at the extrapolated
boundary for U =10 may be very small.
Let us discuss the frequency dependence of σreg(ω)
obtained by eq. (2) for the insulating and metallic cases
wherethe systemisinthe criticalregion;we presentour
resultsinFig.9.Figures9(a)and9(b)depictthe results
for insulating cases. The structures of σreg(ω) between
N = 16 and N = 20 are similar for each U. This fact
s s
suggests that finite-size effects are rather irrelevant for
the structure of σreg(ω). For each U, broad responses
corresponding to the transitions to the upper-Hubbard
band appear around ω = U. For U = 10, a sharp peak
emerges at the lower edge of the band. According to ref.
4, this sharppeak is ascribed to the spin-polaronforma-
tioninthephotoexcitedstate,andtheweightofthepeak
becomes smaller for larger N . On the other hand, our
s
sharp peak for U =10 does not show any significant N
s
dependence,whichimpliesthatitwillsurviveinthether-
modynamiclimit,althoughthe peaksgetrapidlyweaker
for largerU =20 and 35. The broad responses are unaf-
fectedbythechangeinU betweenU =20and35except
that the center of responses moves with U. Figures 9(c),
9(d) and 9(e) depict the metallic cases; one can observe
not only the high-frequency spectrum due to the transi-
tions to the upper-Hubbardband but also a pronounced
structure in the low-energy region below the Mott gap.
First, let us see the transitions to the upper-Hubbard
band.Byswitchingδ on,these transitionsgraduallylose
their weights by transferring them to the lower-energy
structure, as already been reported in refs. 2, 3 and 4.
Next, let us see the structure in the low-frequency re-
gion. We find strong responses growing rapidly towards
ω =0 with decreasing ω. It is noticeable that the shape
seems to be a common property irrespective of U and δ
within the present parameter ranges. The same shapes
are produced in the present results of both Ns =18 and Fig. 9. Frequency dependence ofregularpartoftheoptical con-
20.Theω-dependenceoftheshaperemindsusofthe1/ω ductivity in (a) for Ns = 16 and Ne = 16, namely δ = 0, (b)
tail and the mid-IR peak observedin the experiments of for Ns = 20 and Ne = 20, namely δ = 0, (c) for Ns = 18 and
high-T compound.20 Ne=16,namelyδ=1/9,(d)forNs=18andNe=16,namely
c
Finally,letusmakecomparisonswithexperimentalre- δ=1/9and(e)forNs=20andNe=16,namelyδ=0.2.Blue,
green and red lines denote the cases for U = 10, U = 20 and
sultsofopticalconductivities.Recallherethatthecoher- U =35, respectively. Delta functions arebroadened withwidth
ent and incoherent parts are not definitely separated in ofη=0.05.
the experimentally observed optical conductivity. Thus,
it is not so easy to make direct comparison between ex-
perimental and numerical σreg(ω). Theoretically, on the the system is in the insulating state. This comes from
otherhand,onehastointroducethebroadeningofdelta the tiny D due to the finite-size effect. For δ > 0, on
functions with some widths in dealing with finite sys- the other hand, the Drude weight is considerably large
tems as the present ones. Keeping these circumstances in comparison with this tiny D at half filling; there oc-
in mind, we consider the frequency dependence of σ(ω) curs a huge peak at ω = 0 as the coherent component.
including the Drude part; Figure 10 depicts results of ItstailinFig.10appearsduetothebroadeningofwidth
U = 10, typical value of the Coulomb interaction for η =0.2.Noteherethateventhetailissizableincompar-
most of the high-T materials. For δ = 0, the present ison with the intrinsic responses of σreg(ω) in the inner-
c
spectrum reveals a clear increase as ω → 0, even when gap region. This means that in order to understand op-
8 J.Phys.Soc.Jpn. FullPaper H.Nakano,Y.Takahashi andM.Imada
lowdopingsisindeedcommontoallthehigh-T cuprates.
c
Itisevenmoreuniversalinthetransitionmetaloxidesin
general.Ourcalculatedresultforthefurtherdopingpro-
cessfrommoderatetooverdopingisalsoconsistentwith
the universaltrend of the experimentalresults in transi-
tion metal oxides including the high-T cuprates, where
c
the mid-IR structure is absorbed to the Drude weight.
The δ2 scaling of the Drude weight is tightly related to
this profound and universal trend in the doping process
of the whole spectral weight in the optical conductivity.
5. Conclusion and Discussion
We have investigated the Drude weight of the Hub-
bardmodel onthe squarelattice by means of large-scale
parallelized exact diagonalizations. The hole-doping de-
Fig. 10. Frequency dependence ofoptical conductivities forU = pendence of the Drude weight near the transition point
10. The insulating case of Ns = 20 and Ne = 20 is denoted by between the metal and the Mott insulator has been a
the sky-blue dotted line. Red, deep-blue, green and black solid
controversialproblem. We have very carefully examined
lines correspond to the metallic cases of Ns =18 and Ne =16,
finite-size effects and we have successfully captured the
Ns =20 and Ne =16, Ns =16 and Ne =12 and Ns =18 and
Ne=12,respectively.Deltafunctionsarebroadenedwithwidth intrinsic behavior of the Drude weight. Our results sug-
of η = 0.2. For comparison, the coherent part in σ(ω) is also gest that the Drude weight is scaled by D ∝ δ2, which
presentedforthetwocasesofNs=18bythebrokenlinesofthe is consistent with results from various kinds of analy-
samecolors.
sis. In particular,it is important that this scaling agrees
with the quantum Monte Carloresults throughthe scal-
ing theory of the metal-insulator transition. The behav-
tical responses of the system, it is important to capture
ior D ∝ δ2 is a characteristic feature near the marginal
boththebehaviorofthecoherentcomponentandthatof
quantumcriticalpointofthetransitionbetweentheMott
the incoherentone included in the total σ(ω) simultane-
insulator and the doped metal in two-dimensional sys-
ously.InFig.10,qualitativelythesamechangebydoping
tems. The restructuring of the electronic structure upon
in the transitions to the upper-Hubbard band is seen as
doping into the Mott insulator first causes the spectral
that mentioned above for U =20 and 35. The structure
weight transfer from the weight above the gap mainly
in this ω region should be associated with the charge
to a mid-IR incoherent component below the gap. The
transfer from the copper d to the oxygen p states in the
spectral weight is further progressively transferred from
case of the cuprates, hence it is, in narrow sense, dif-
the mid-IR component to the Drude weight. This large-
ferent from the present case of the Mott-Hubbard type.
scale restructuring seen in the frequency and doping de-
Aside from the character of the gap, i.e., either it is the
pendences of σ(ω) in our result is consistent with the
Mott-Hubbardgaporthechargetransfergap,theevolu-
experimental results on transition metal oxides includ-
tion of doping observedin the experiments are correctly
inghigh-T cuprates.The δ2 scalingoftheDrude weight
captured in our numerical results, the spectral weight c
comes from this underlying physics involving the high
transfer from the higher to lower energies as well as the
energy scale.
collapse of the sharp peak at the low-energy edge.
AtU =0,theHubbardmodelwestudiedhasaperfect
Theevolutionofthelow-energystructurewithdoping
nesting just at half filling and the Fermi level coincides
contains both the growth of the Drude weight and the
with the van Hove singularity. The effect of deviation
developmentofthe incoherentmid-IRstructure.Theto-
from this perfect nesting by introducing the hopping to
talsum ofthese two contributionsis the effective carrier
further sites (such as the next-nearest-neighbor hopping
density N , whose weight is proportional to the dop-
eff
t′) isleftfor separatestudy.The QMCstudy14 indicates
ing concentration. In the doping process of the evolu-
that the scalingbehavior of the compressibility does not
tion from the insulator to low doping, the mid-IR in-
sensitivelydependonthepresenceorabsenceoft′,which
coherent part below the Mott gap grows quickly in our
implies thatthe presentscalingofthe Drude weightalso
calculated results and is already prominent even in un-
holds. Next, let us recall the marginal point with z = 4
derdoped systems (see the result at δ ∼ 0.111), with a
in the one-dimensional Hubbard model with the next-
rather insensitive dependence on the further doping to
nearest-neighbor hopping10,19 mentioned above. In this
the overdoped systems (see δ = 0.333). On the contrary
case, a singularity in the density of states is absent at
to the mid-IR part, the Drude part grows slowly at low
half filling at U = 0 because the next-nearest-neighbor
doping, while it starts growing quickly from moderate
hopping is small. This fact also suggests that the metal-
to over doping, which is the origin of the scaling given
by D ∝ δ2 rather than D ∝ δ. The total sum resulting insulatortransitionwithz =4occursirrespectiveofsuch
a singularity in the density of states.
in N evolves smoothly and linearly with δ because of
eff
The finite-size effects in such a quantity as the Drude
the compensation of these two contributions. The fact
weight are considerably large. Thus, the careful exam-
thattheDrudepartremainsonlyasmallfractionofN
eff
ination of the finite-size effects is necessary in possible
whereasthemid-IRstructurerapidlygrowsfromzeroto
ways as much as we can. If one calculates the quantity
J.Phys.Soc.Jpn. FullPaper H.Nakano,Y.TakahashiandM.Imada 9
without considering the effects, obtained results could 2) E.Dagotto, A.Moreo,F.Ortolani,D.PoilblancandJ.Riera:
severelysufferfromsuperficialfinite-sizeeffects.Bysuch Phys.Rev.B45(1992)10741.
results,unfortunately, the incorrectconclusionwould be 3) H.NakanoandM.Imada:J.Phys.Soc.Jpn.68(1999) 1458.
4) T.Tohyama,Y.Inoue,K.TsutsuiandS.Maekawa:Phys.Rev.
derived. In this meaning, the present problem of the
B72(2005)045113.
Drude weight is an instructive example. As a way to re-
5) The largest dimension of the subspace of the Hubbard model
duce the finite-size effect, the present work employs ap- with a given Ns is realized at half filling with the same num-
propriatelychosenclustersundertheappropriatebound- berofup-spinelectronsanddown-spinones.Forexamples,the
ary conditions. Thereby, we have captured the critical dimensions are 165 636 900 for Ns = 16, 2 363 904 400 for
behavior of D ∝ δ2 near the half-filling insulator. Av- Ns = 18 and 34 134 779 536 for Ns = 20. It is possible to
reducethesedimensionsbyusingthesymmetriesoftheHamil-
eraging over twisted boundary conditions21 may be an-
tonian;however,theprogramforthereducedcaseisnotneces-
other possibility for the reduction. Quite recently such sarilyapplicable to any Hamiltonians having arbitrary shapes
an investigationhas been done in the N =20 cluster of ofclustersincontrasttoourprogram.Theversatilefeatureof
s
thet-J modelplusthree-siteterms.22Althoughthiswork our code makes it possible to examine the various cases very
easily.Althoughthedimensionisveryhuge,ourcodecosts29.7
givesaresultoftheδ-dependenceofD,itissosubtlethat
hoursforthe260-timeiterationstoobtaintheground-stateen-
one cannot conclude either of D ∝ δ2 or D ∝ δ. In the
ergy of the 20-site half-filling Hamiltonian for U = 35, when
not-too-distantfuture,computerswillbeevenmorepow- we use 16 nodes in the supercomputer SR11000 in Hokkaido
erful; numerical diagonalizations of the Hubbard model University.
withlargersystemsizeswouldbe possible.Suchcalcula- 6) E.GaglianoandC.Balseiro:Phys.Rev.Lett.59(1987)2999.
7) N.FurukawaandM.Imada:J.Phys.Soc.Jpn.60(1991)3604;
tions could clarify the validity of the present procedure
ibid.61(1992)3331.
for reducing the finite-size effects in D and provide us 8) H.NakanoandY.Takahashi:J.Phys.Soc.Jpn.72(2003)1191.
withamoredefinite resultforthe criticalbehaviorofD. 9) H.NakanoandY.Takahashi:J.Mag.Mag.Mat.272-276(2004)
We hope that the knowledge in the present study will
487.
serve for such future studies. 10) H.NakanoandY.Takahashi:J.Phys.Soc.Jpn.73(2004)983.
11) The negative Drude weights with large absolute values were
Acknowledgements reportedatδ=0andδ=1/8forNs=16inref.2.Inref.4,the
discontinuousdependencewasreportedatδ=1/3forNs=18.
These works use only the periodic boundary condition. Note
ThisworkwassupportedbyaGrant-in-AidforYoung
thatsuchbehavioroccursonlyattheholeconcentrationwhen
Scientists (B) from the Ministry of Education, Culture, theelectronicshellstructureisopen.
Sports, Science and Technology of Japan (15740221). A 12) H. Tsunetsugu and M. Imada: J. Phys. Soc. Jpn. 67 (1998)
part of the computations was performed using facilities 1864.
13) M.Imada: J.Phys.Soc.Jpn.64(1995) 2954.
of the Information Initiative Center, Hokkaido Univer-
14) N.FurukawaandM.Imada:J.Phys.Soc.Jpn.62(1993)2557.
sity and the Supercomputer Center, Institute for Solid
15) F.AssaadandM.Imada:Phys.Rev.Lett.76(1996) 3176.
State Physics (ISSP), University of Tokyo. The authors 16) M.Imada: J.Phys.Soc.Jpn.74(2005) 859.
thank Prof. Hajime Takayama and the staffs of the Su- 17) M.Imada: Phys.Rev.B72(2005) 075113.
percomputer Center, ISSP for their support to perform 18) T. Misawa, Y. Yamaji and M. Imada: J. Phys. Soc. Jpn. 75
(2006) 083705.
our large-scale parallelized calculations in a special job
19) H.Nakano,Y.TakahashiandM.Imada:PhysicaB359-361C
class of the supercomputer of ISSP.
(2005) 657.
20) S. Uchida, T. Ido, H. Takagi, T. Arima, Y. Tokura and S.
Tajima:Phys.Rev.B43(1991) 7942.
1) M. Imada, A. Fujimori and Y. Tokura: Rev. Mod. Phys. 70 21) D.Poilblanc:Phys.Rev.B44(1990) 9562.
(1998)1039andreferencestherein. 22) T.Tohyama:J.Phys.Chem.Solids67(2006)2210.