有限溫度下一維Hubbard模型的化學勢泛函理論研究?

陸展鵬 魏興波 劉天帥 陳阿海 高先龍

(浙江師范大學物理系,金華 321004)

有限溫度下一維Hubbard模型的化學勢泛函理論研究?

陸展鵬 魏興波 劉天帥 陳阿海 高先龍?

(浙江師范大學物理系,金華 321004)

(2017年1月7日收到;2017年3月26日收到修改稿)

通過數值方法求解了有限溫度下一維均勻Hubbard模型的熱力學Bethe-ansatz方程組,得到了在給定溫度和相互作用強度情況下,比熱c、磁化率χ和壓縮比κ隨化學勢μ的變化圖像.基于有限溫度下一維均勻Hubbard模型的精確解,利用化學勢(μ)-泛函理論研究了一維諧振勢下的非均勻Hubbard模型,給出了金屬態和Mott絕緣態下不同溫度情況時局域粒子密度ni和局域壓縮比κi隨格點的變化情況.

熱力學Bethe-ansatz方程組,壓縮比,Hubbard模型,化學勢泛函理論

1 引言

20世紀初,玻色和愛因斯坦提出了著名的玻色-愛因斯坦凝聚(Bose-Einstein condensation,BEC)理論,他們預言,對于玻色氣體,當溫度低于某一個臨界值Tc時,所有的粒子將會占據能量最低的量子態[1].隨著1985年美國華裔科學家朱棣文等[2]利用激光冷卻技術獲得了接近絕對零度的溫度,冷原子物理開始成為研究的熱點.之后在1995年,Ketterle,Cornell和Wiemann首次在實驗中實現了Na和Rb的BEC凝聚體[3,4].而后人們開始設法冷卻費米氣體并得到費米凝聚體.由于費米氣體受到泡利不相容原理的限制,所以冷卻難度比玻色氣體更大.在1999年,Jin的研究組[5]實現了費米凝聚.隨著Feshbach共振技術的發展,人們可以通過磁場改變束縛態原子的散射長度,進而調節原子之間的相互作用,從而大大推動了冷原子實驗的發展[6,7].另外由于光晶格技術的發展,人們在冷原子體系中可以構造出非常干凈的一、二、三維晶格結構[8].

一維體系呈現出豐富的物理性質,如20世紀,Tomonaga[9]和Luttinger[10]提出的著名的Luttinger液體理論.他們指出一維費米體系之中,如果考慮相互作用,則朗道液體理論不再適用,體系中只有集體激發,沒有單粒子激發,其中最明顯的特征即自旋電荷分離[11].很多一維多粒子體系具有精確可解性,如1931年Bethe[12]提出的Bethe-ansatz理論,可求解一維Heisenberg模型;1963年,Lieb和Liniger[13]求解了一維玻色氣體;1967年,Yang[14]和Gaudin[15]求解了自旋為1/2的費米氣體;1968年,Lieb和Wu[16]求解了一維Hubbard模型.以上理論結果可以與一維冷原子實驗相互印證,并用于理解一維諧振勢受限體系中呼吸模的集體激發情況[17]和自旋電荷分離現象[18,19]等.

由于絕對零度無法達到,考慮到溫度對受限體系的影響具有十分現實的意義,所以本文研究一維有限溫度的晶格冷原子晶格體系.在1969年,楊振寧等[20]提出了研究一維可積模型的熱力學Bethe-ansatz(TBA)方法,隨后,日本物理學家Takahashi[21,22]研究了XXX鏈模型和具有排除相互作用的自旋?1/2模型的熱力學性質.此后,一維有限溫度的研究成為熱點,如利用TBA方程組結合Haldane統計討論低溫下一維強相互作用任意子簡并特征[23]和能態性質[24],研究有限溫度下吸引相互作用的費米氣體的配對特征和相變性質[25,26],有限溫下自旋為1的玻色子在一維諧振勢中的量子臨界現象[19,27],以及有限溫度下一維吸引相互作用的費米氣的Tomonaga-Luttinger液體相變[19,28]等.

由于冷原子實驗中的體系普遍是非均勻的,體系一般變得不可精確求解.利用一維均勻體系的精確解并結合局域密度近似的密度泛函理論求解一維非均勻體系,是一種可行的方法[29?32].利用此密度泛函理論可以研究有限溫度下一維晶格體系的交換關聯勢[29],有限溫度下一維諧振勢Hubbard模型的磁導率等熱力學量[33],以及用含時密度泛函理論研究一維費米體系自旋電荷分離[34,35]等性質.

基于以上的研究,我們數值求解了有限溫度下一維均勻Hubbard模型的TBA方程組[36],主要分析了在不同相互作用和不同溫度下,磁化率χ、比熱c和壓縮率κ等熱力學量,之后基于這些一維均勻體系的性質,以化學勢為泛函結合密度泛函理論[37]方法研究了一維非均勻Hubbard模型,得到了不同溫度下金屬相和Mott絕緣相的粒子密度和局域壓縮比隨格點的變化情況.研究結果表明,在低溫情況下,一維均勻體系中的磁化率χ、比熱c和壓縮率κ能夠度量四個相變區域,分別為真空態、金屬態、Mott絕緣態和帶絕緣態.隨著溫度的升高,體系熱漲落逐漸占據主導地位,取代量子漲落成為主導,Mott絕緣體相會消失.在非均勻體系中,體系處于不同相的共存相.隨著溫度的升高,熱漲落會破壞勢阱中間的Mott絕緣相,相應的粒子密度分布ni和局域壓縮比κi的變化反映了這一過程.

2 理論模型

本文先討論一維均勻Hubbard模型,其哈密頓量的形式如下：

式中,t為近鄰格點的跳躍能;U為同格點粒子之間的相互作用,本文只研究U>0的排斥相互作用情況;μ為化學勢;h為磁場;為i格點自旋為σ的產生算符,表示為i+1格點自旋為σ的湮滅算符;為粒子數算符.

比熱c、磁化率χ和壓縮比κ的表達式如下：

熱力學勢g所滿足的TBA耦合方程的具體表達式及其相關的二階微分量的求解過程可以參考文獻[36].

我們首先通過數值方法求解一維均勻Hubbard模型的TBA耦合方程組,得到磁化率χ、比熱c、壓縮率κ、化學勢μ、粒子密度分布n等熱力學量,通過分析這些物理量了解系統的量子相變.基于對這些熱力學量的分析,我們研究一維受限諧振勢下的Hubbard模型,

式中,Vext為外勢強度,Na為總格點數.

對于一維非均勻Hubbard模型的求解,可以利用一維晶格密度泛函理論,即通過構造建立在Bethe-ansatz嚴格解基礎上的交換關聯勢,用基于局域密度近似來求體系的定態或動力學性質(BALDA方法[11,17,29,33,34])來研究.但是上述方法在處理具有Mott絕緣相的情況時,由于此時構造的交換關聯勢在n=1的半充滿情況下具有不連續性(即所謂的Mott能隙),此時基于BALDA的Kohn-Sham方程在不連續的交換關聯勢上下振蕩,很難收斂,從而產生不能自洽求解的問題[37].考慮到密度和化學勢之間具有一一對應的函數關系,基于密度的泛函理論等價于基于化學勢(μ)的泛函理論,相應的BALDA變成μ-BALDA,以化學勢為泛函的交換關聯勢是連續的,于是上述困難可以解決.

下面我們具體闡述μ-BALDA方法,討論如何通過一維均勻體系的化學勢μ(U,n)作為泛函,構造出非均勻體系的有效勢Veff,并利用Kohn-Sham方程求解出非均勻體系的局域粒子密度ni.

我們先構造出均勻體系的交換關聯勢Vhxc[n],

其中μ(U,n)代表均勻體系的化學勢(以下簡記為μ),μ(0,n)為均勻體系中U=0時的化學勢,即動能.μ(U,n)可以通過如下公式求出:

式中e(U,n)是均勻體系中單位格點的能量.考慮到(10)式中的粒子密度和化學勢具有一一對應的函數關系,n可以被反解出來,表示成以化學勢μ為自變量的函數,

于是方程(9)中的交換關聯勢Vhxc[n]就可以表示成以化學勢μ為泛函的函數Vhxc[μ],

基于密度泛函理論,多粒子體系的薛定諤方程等價于一個在有效勢場中的單粒子Kohn-Sham(KS)方程,

這里tij=t(當j=i±1)而其他情況下為0,其中系統的密度分布是由KS方程的軌道波函數自洽給出,

式中的α為能級指標,2來自于自旋簡并,f(ω)={1+exp[β(ω?μ)]}?1是費米分布,β =1/(kBT)是溫度的倒數,其中kB是玻爾茲曼常量.在實際計算中,非均勻體系的可由均勻系統的Vhxc[μ]在局域密度近似下給出：

從方程(13)和(14)中計算密度分布等價于尋找以化學勢為泛函的不動點問題：

如前所述,對這些方程的普通迭代求解辦法不能解決由于交換關聯勢不連續帶來的不收斂問題.而更精妙的迭代方法如Newton-Rhapson方法也不能解決該問題,因為它用到了非線性方程對于未知變量的一級導數,其在不連續處會趨于無窮,同樣會導致耦合方程的不收斂問題.這里我們擬用多維二分法來解決上述以化學勢為泛函的不動點問題,即通過構造

來找方程左邊的函數與橫軸的交點,即其零解,詳細的數值解法參見文獻[29].

3 結果分析與討論

3.1 一維均勻Hubbard模型

本節先給出一維均勻Hubbard模型中比熱c、磁化率χ和壓縮比κ在不同溫度和相互作用強度下的結果,分析這些物理量在度量量子相變中的作用.圖1給出了低溫(T=0.1)情況、不同相互作用強度U下的比熱c、磁化率χ和壓縮比κ隨化學勢μ變化的性質.

圖1中比熱c和壓縮比χ的圖像中均出現了四個峰值點或拐點,這四個峰值點分別對應著系統的四次相變:真空態向金屬態轉變,金屬態向Mott絕緣體態轉變,Mott絕緣體態向金屬態轉變,金屬態向帶絕緣體態轉變.從壓縮比κ的圖像看出,當相互作用U不斷增加時,壓縮比κ中間出現了為零的區域,這部分區域對應Mott絕緣體相(絕緣相不可壓縮),區域的大小也反映了Mott能隙的大小,隨著相互作用U的增大,Mott絕緣體區域的范圍也在增大.

圖1 (網刊彩色)溫度T=0.1,不同相互作用強度(U=1,4,8,15)下的壓縮比κ、比熱c和磁化率χ隨化學勢μ的變化情況,其中對化學勢進行了平移μ?U/2,使系統具有中心對稱(數值計算方便取t=kB=1)Fig.1.(color online)Speci fi c heat c,susceptibility χ,and compressibility κ as a function of the chemical potentialμfor temperature T=0.1 and di ff erent interaction strength(U=1,4,8,15).In order to have the system of the central symmetry,we shift the chemical potentialμby U/2(for convenience,we let t=kB=1).

圖2 (網刊彩色)T=0.1,U=1時壓縮比κ和粒子密度n隨μ的變化(圖中的虛線從左到右把系統劃分為真空區、金屬相和帶絕緣體相)Fig.2.(color online)Compressibility κ and density n as a function of the chemical potentialμfor temperature T=0.1 and interaction strength U=1.The dotted lines from left to right divide the system into the vacuum phase,the metal phase,and the bandinsulating phase.

圖3 (網刊彩色)T=0.1,U=4壓縮比時κ和粒子密度n隨μ的變化(圖中的虛線從左到右把系統劃分為真空區、金屬相、Mott絕緣體相、金屬相和帶絕緣體相)Fig.3.(color online)For the temperature T=0.1 and interaction strength U=4,compressibility κ,and density n as a function of the chemical potentialμ.The dotted line from left to right divide the system into the vacuum phase,the metal phase,the Mott-insulating phase,the metal phase,and the bandinsulating phase.

為了清晰地看出Mott絕緣相的區域,我們將不同相互作用U下壓縮比κ和粒子密度n的結果放在一起比較,結果見圖2、圖3和圖4.從圖2可以看出,在低溫(T=0.1)下,相互作用U較弱時,此時并非Mott絕緣體相(這與零溫時有限相互作用在半充滿情況下即為Mott絕緣相不同).此時壓縮比κ的圖像中沒有反映出不可壓縮的區域,而粒子密度n的圖像上也沒有出現n=1的半充滿的Mott平臺.從圖3可以看出,由于粒子之間的相互作用U增強到4,壓縮比κ的圖像中間出現為零的不可壓縮的Mott絕緣區域,而對應的粒子密度n也出現了n=1的半充滿的Mott平臺(圖中的中間兩條虛線部分對應的是Mott絕緣體系區域的范圍).從圖4可以看出,隨著粒子之間的相互作用U增強到8,壓縮比κ的圖像中Mott絕緣體系區域的范圍擴大了,而對應的粒子密度n的圖像中的Mott平臺也相應地增大,兩者的變化是一致的.由此說明壓縮比κ可以分辨各個量子相.如果將圖1中的磁化率χ和比熱c的結果和粒子密度n相比較也能夠得出上述類似的結論.

圖5給出了高溫情況(T=1)、不同相互作用(U=1,4,8,15)下的壓縮比κ、比熱c和磁化率χ隨化學勢μ變化的結果.

從圖5中我們發現熱力學漲落破壞了原來的相變,四個相變點已經基本消失.在非常強的相互作用(U=15)下,Mott絕緣體區域仍然穩健,這可以從壓縮比κ的圖像中看出,其不可壓縮的κ為零的區域依然存在,此時熱力學漲落仍然沒有破壞Mott絕緣體相.從上面的結果中我們發現低溫下熱力學量壓縮比κ、比熱c和磁化率χ可以清晰地反映出一維均勻Hubbard模型的相變情況,拐點和峰值點可以反映體系的相變.隨著溫度升高,熱漲落會產生影響并很快破壞Mott絕緣相,但如果相互作用足夠大,Mott絕緣體區域在高溫(T=1.0)時仍然穩健.更詳細的相圖可由Bethe-ansatz解析解分析得到,具體可見文獻[19,26,27].

圖4 (網刊彩色)T=0.1,U=8壓縮比時κ和粒子密度n隨μ的變化(圖中的虛線從左到右把系統劃分為真空區、金屬相、Mott絕緣體相、金屬相和帶絕緣體相)Fig.4.(color line)Compressibility κ and density n as a function of the chemical potentialμfor temperature T=0.1 and interaction strength U=8.The dotted lines from left to right divide the system into the vacuum phase,the metal phase,the Mott-insulating phase,the metal phase,and the band-insulating phase.

圖5 (網刊彩色)溫度T=1.0,不同相互作用強度(U=1,4,8,15)下的壓縮比κ、比熱c和磁化率χ隨化學勢μ的變化Fig.5.(color online)Speci fi c heat c,susceptibility χ,and compressibility κ as a function of the chemical potentialμfor temperature T=0.1 and di ff erent interaction strength(U=1,4,8,15).

3.2 一維束縛勢中的Hubbard模型

本節分析光晶格中的一維費米冷原子體系,由于束縛勢的存在,體系變成非均勻的.在單帶近似和緊束縛模型下,體系可由非均勻的Hubbard模型[方程(8)]描寫.下面我們討論金屬相和Mott絕緣相中的局域粒子密度ni和局域壓縮比κi在不同溫度下的分布.

圖6是相互作用U=8,粒子數N=30時,不同溫度(T=0.01,0.1,0.6,1.0)下金屬態局域粒子密度ni隨格點i的變化.其他物理參數為:晶格長度Na=200,諧振勢強度Vext=2.5×10?3,這些參數也用在了以下其他圖的計算中.

圖6 (網刊彩色)相互作用U=8,粒子數N=30時,不同溫度T=0.01,0.1,0.6,1.0下金屬態的局域粒子密度ni隨著格點i的變化Fig.6.(color online)The local density niof the metal phase as a function of the lattice site i for interaction strength U=8,the number of particles N=30,and di ff erent temperature T=0.01,0.1,0.6,1.0.

從圖6可以看出,溫度較低(T=0.01)時,量子漲落引起的Friedel振蕩仍然可見,隨著溫度的升高,熱力學漲落占據主導作用,Friedel振蕩不復存在,密度分布呈高斯型.

圖7為相互作用U=8,粒子數N=70時,不同溫度(T=0.01,0.1,0.6,1.0)下Mott絕緣態局域粒子密度ni的分布情況,從圖7可以看出,在溫度比較低的情況下,整個體系是金屬相和Mott絕緣相共存,Mott平臺(ni=1)存在于勢阱中部,兩邊屬于金屬相.溫度T從0.01升高至0.1時,Mott絕緣相整體保持穩定,逐漸從Mott平臺的兩邊開始受到破壞.分析表明,Mott平臺穩健存在于T<0.3,破壞Mott相的溫度和均勻體系相近.當溫度T升高到1.0時,熱力學漲落占主導作用,體系的密度分布接近于高斯型.

為了更好地刻畫溫度和相互作用對Mott絕緣相的影響,我們計算了局域壓縮比,定義為κi= κ(n)|n→ni.圖8為溫度T=0.1,相互作用U=8時,局域粒子密度分布ni和局域壓縮比κi的對比圖.從圖8可以看出,體系的局域壓縮比κi為零的不可壓縮區域對應著Mott絕緣平臺(虛線分割的中間部分正好對應Mott絕緣區)和真空態,分別對應著密度分布ni=1和ni=0.圖9為高溫情況T=1.0的局域壓縮比κi和局域粒子密度ni的對比,可以看出局域壓縮比在中間區域有下降的趨勢但不為零,這正說明了隨著溫度的升高,由于相互作用和束縛勢產生的Mott絕緣平臺被破壞,體系重新變成了可壓縮的金屬相,相應的密度為ni<1.

圖7 (網刊彩色)相互作用U=8,粒子數N=70時,不同溫度T=0.01,0.1,0.6,1.0下金屬相和Mott絕緣相共存時的局域粒子密度ni隨著格點i的變化Fig.7.(color online)The local density niof the metal-Mott insulating mixed phase as a function of the lattice site i for interaction strength U=8,the number of particles N=70,and di ff erent temperature T=0.01,0.1,0.6,1.0.

圖8 (網刊彩色)溫度T=0.1,相互作用U=8,粒子數N=70時,局域粒子密度ni和局域壓縮比κi隨著格點i的變化 (圖中的虛線從左到右把系統劃分為真空態、金屬相、Mott絕緣體相、金屬相、真空態)Fig.8.(color online)The local density niand local compressibility κias a function of the lattice site i for temperature T=0.1,interaction strength U=8,and the number of particles N=70.The dotted line from left to right divide the system into the vacuum phase,the metal phase,the Mott-insulating phase,the metal phase,and the vacuum phase.

圖9 (網刊彩色)溫度T=1.0,相互作用U=8,粒子數N=70時,局域粒子密度ni和局域壓縮比κi隨著格點i的變化Fig.9.(color online)The local density niand the local compressibility κias a function of the lattice site i for temperature T=1.0,interaction strength U=8,and the number of particles N=70.

為了更好地理解有限溫度下金屬相向Mott絕緣相的轉變過程,我們給出了一維束縛勢中的Hubbard模型在給定粒子數N=70、晶格長度Na=200和諧振勢強度Vext=2.5×10?3下的T-U相圖,如圖10所示.

圖10 (網刊彩色)一維束縛勢中的Hubbard模型在給定粒子數N=70、晶格長度Na=200和諧振勢強度Vext=2.5× 10?3下的溫度-相互作用(T-U)相圖(圖中橫坐標表示相互作用U,縱坐標表示溫度T,圖中的點是兩個相的分界點)Fig.10.(color online)T-U phase Diagram of a harmonically trapped Hubbard model in one dimension with the number of particles N=70,the length of the lattice Na=200 and the harmonic potential strength Vext=2.5×10?3.The points in the divide the region into two quantum phases.

從圖10中我們可以看出,Mott絕緣相只存在于溫度T 6 0.3和相互作用U>3.2的區域,隨著溫度的升高,體系的Mott絕緣相會被破壞掉.

4 總 結

本文通過數值方法求解了一維均勻Hubbard模型的TBA方程組,得到了不同溫度和相互作用下壓縮比κ、比熱c和磁化率χ等熱力學量.發現在低溫下,這三個熱力學量(特別是壓縮比κ)能夠清晰地反映出體系的四個相變,分別對應著體系從真空態向金屬態轉變、金屬態向Mott絕緣體態轉變、Mott絕緣體態向金屬態轉變、金屬態向帶絕緣體態轉變.隨著溫度的升高,熱漲落逐漸取代量子漲落成為主導,溫度會破壞掉原有體系的Mott絕緣體相,明顯的相變點消失.而后利用一維均勻Hubbard模型的Bethe-ansatz嚴格解的數值結果,我們構造了可以用于研究一維非均勻Hubbard模型的化學勢泛函方法:μ-BALDA,該方法從理論上講是嚴格的密度泛函理論的變種.化學勢泛函理論的構造解決了密度泛函理論中Kohn-Sham耦合方程在交換關聯勢不連續時的不收斂問題,并把問題轉變成用多維二分法來解決以化學勢為泛函的不動點問題,從而收斂性可以得到保障.

我們用μ-BALDA方法我們數值求解了有限溫度下一維諧振勢中的Hubbard模型,得到了不同溫度下粒子的密度分布情況.從結果中可以看出在給定相互作用U的情況下,逐漸升高溫度,熱漲落會破壞金屬態中的Friedel振蕩和Mott絕緣體平臺,最終變成平滑的高斯分布.而Mott絕緣體在一定溫度范圍內保持穩定,平臺的破壞是從邊界開始到中心地帶,從而使整個Mott平臺被破壞.最后通過對比非均勻體系的局域壓縮比κi,并與粒子密度分布ni進行了對比,驗證了溫度對Mott絕緣態的破壞.在給定粒子數、晶格長度和諧振勢強度下,我們給出了一個簡單的T-U相圖.由相圖可知在此情況下Mott絕緣相只存在于溫度T 6 0.3和相互作用U>3.2的區域內.

本文所運用的μ-BALDA方法可以進一步推廣到非均勻的具有短程相互作用的Gaudin-Yang模型和吸引Hubbard模型中去,其中自旋激發有能隙,相應的交換關聯勢在自旋為零處不連續.另外也可以推廣到一維XXZ模型 (可用Jordan-Wigner變換變成一維近鄰相互作用的極化費米子模型,在Tomonaga-Luttinger液體相和公度的電荷密度波相的臨界點處n=0.5,交換關聯勢不連續)[38]中去,或推廣到有限溫度(或多組分)的上述模型[26]中去.

感謝浙江師范大學物理系王沛副教授在論文構思階段部分物理內容的討論.

參考文獻

[1]Wang Z C 2003 Thermodynamics Statistical Physics(Beijing:Higher Education Press)p300(in Chinese)[汪志誠1993 熱力學和統計物理學(北京:高等教育出版社)第300頁]

[2]Chu S,Hollberg L,Bjorkholm J E,Cable A,Ashkin A 1985 Phys.Rev.Lett.55 48

[3]Davis K B,Mewes M O,Andrews M R,van Druten N J,Durfee D S,Kurn D M,Ketterle W 1995 Phys.Rev.Lett.75 3969

[4]Anderson M H,Ensher J R,Matthews M R,Wieman C E,Cornell E A 1995 Science 269 198

[5]DeMarco B,Jin D S 1999 Science 285 1703

[6]Feshbach H 1958 Ann.Phys.5 357

[7]Batchelor M T,Bortz M,Guan X W,Oelkers N 2005 Phys.Rev.A 72 061603

[8]Pachos J K,Knight P L 2003 Phys.Rev.Lett.91 107902

[9]Tomonaga S 1950 Prog.Theo.Phys.5 544

[10]Luttinger J M 1963 J.Math.Phys.4 1154

[11]Gao X L 2010 Phys.Rev.B 81 104306

[12]Bethe H 1931 Z.Phys.71 205

[13]Lieb E H,Liniger W 1963 Phys.Rev.130 1605

[14]Yang C N 1967 Phys.Rev.Lett.19 1312

[15]Gaudin M 1967 Phys.Lett.A 24 55

[16]Lieb E H,Wu F Y 1968 Phys.Rev.Lett.20 1445

[17]Hu H,Gao X L,Liu X J 2014 Phys.Rev.A 90 013622

[18]Lee J Y,Guan X W,Sakai K,Batchelor M T 2012 Phys.Rev.B 85 085414

[19]Guan X W,Batchelor M T,Lee C 2013 Rev.Mod.Phys.85 1633

[20]Yang C N,Yang C P 1969 J.Math.Phys.10 1115

[21]Takahashi M 1969 Prog.Theo.Phys.42 1098

[22]Takahashi M 1972 Prog.Theo.Phys.47 69

[23]Batchelor M T,Guan X W 2006 Phys.Rev.B 74 195121

[24]Batchelor M T,Guan X W,Oelkers N 2006 Phys.Rev.Lett.96 210402

[25]Guan X W,Batchelor M T,Lee C,Bortz M 2007 Phys.Rev.B 76 085120

[26]Jiang Y Z,Chen Y Y,Guan X W 2015 Chin.Phys.B 24 050311

[27]Kuhn C C N,Guan X W,Foerster A,Batchelor M T 2012 Phys.Rev.A 86 011605

[28]Guan X W,Lee J Y,Batchelor M T,Yin X G,Chen S 2010 Phys.Rev.A 82 021606

[29]Gao X L,Chen A H,Tokatly I V,Kurth S 2012 Phys.Rev.B 86 235139

[30]Gao X L 2012 J.Phys.B 45 225304

[31]Gao X L,Asgari R 2008 Phys.Rev.A 77 033604

[32]Hu J H,Wang J J,Gao X L,Okumura M,Igarashi R,Yamada S,Machida M 2010 Phys.Rev.B 82 014202

[33]Campo V L 2015 Phys.Rev.A 92 013614

[34]Gao X L,Polini M,Rainis D,Tosi M P,Vignale G 2008 Phys.Rev.Lett.101 206402

[35]Li W,Gao X L,Kollath C,Polini M 2008 Phys.Rev.B 78 195109

[36]Takahashi M,Shiroishi M 2002 Phys.Rev.B 65 165104

[37]Ying Z J,Brosco V,Lorenzana J 2014 Phys.Rev.B 89 205130

[38]Wang C J,Chen A H,Gao X L 2012 Acta Phys.Sin.61 127501(in Chinese)[王嬋娟,陳阿海,高先龍 2012物理學報61 127501]

PACS:67.25.bd,71.10.–wDOI:10.7498/aps.66.126701

Chemical potential-functional-theory about the properties of one-dimensional Hubbard model at fi nite temperature?

Lu Zhan-Peng Wei Xing-Bo Liu Tian-Shuai Chen A-Hai Gao Xian-Long?

(Department of Physics,Zhejiang Normal University,Jinhua 321004,China)

7 January 2017;revised manuscript

26 March 2017)

In this paper,we numerically solve the thermodynamic Bethe-ansatz coupled equations for a one-dimensional Hubbard model at fi nite temperature and obtain the second order thermodynamics properties,such as the speci fi c heat,compressibility,and susceptibility.We fi nd that these three quantities could embody the phase transitions of the system,from the vacuum state to the metallic state,from the metallic state to the Mott-insulating phase,from the Mott-insulating phase to the metallic state,and from the metallic state to the band-insulating phase.With the increase of temperature,the thermal fl uctuation overwhelms the quantum fl uctuations and the phase transition points disappear due to the destruction of the Mott-insulating phase.But in the case of the strong interaction strength,the Mott-insulating phase is robust,embodying the compressibility.Furthermore,we study the thermodynamic properties of the inhomogeneous Hubbard model with trapping potential.Making use of the Bethe-ansatz results from the homogeneous Hubbard model,we construct the chemical potential-functional theory(μ-BALDA)for the inhomogeneous Hubbard model instead of the commonly used density-functional theory,in order to solve the in-convergence problem of the Kohn-Sham equation in the case of the divergence appearing in the exchange-correlation potential.We further point out a multi-dimensional bisection method which changes the Kohn-Shan equation into a problem of fi nding the fi xed points.Throughμ-BALDA we numerically solve the one-dimensional homogeneous Hubbard model of trapping potential.The density pro fi le and the local compressibility are obtained.We fi nd that at a given interaction strength,the metallic phase and the Mottinsulating phase are destroyed and the density pro fi le becomes a Guassian distribution with increasing temperature.To the metallic phase,Friedel oscillation caused by quantum fl uctuations is still visible at low temperature.With increasing temperature,Friedel oscillation will disappear.This situation re fl ects the fact that the thermal fl uctuation overwhelms the quantum fl uctuations.For the Mott-insulating phase,the Mott-insulating plateau is robust at a certain temperature and only the boundary of the Mott-insulating plateau is destroyed.With increasing temperature,the Mott insulating plateau will be destroyed.And the change of the local compressibility provides the information about such a change.So we conclude that the thermal fl uctuation destroys the original quantum phase.Through our analysis,we fi nd that theμ-BALDA can be used to study the fi nite temperature properties for the system of the exchange-correlation potential divergence with high efficiency.

thermodynamic Bethe-ansatz equations,compressibility,Hubbard model,chemical potentialfunctional theory

10.7498/aps.66.126701

?國家自然科學基金(批準號:11374266)、新世紀優秀人才支持計劃和浙江省自然科學基金(批準號:Z15A050001)資助的課題.

?通信作者.E-mail:gaoxl@zjnu.edu.cn

?2017中國物理學會Chinese Physical Society

http://wulixb.iphy.ac.cn

*Project supported by the National Natural Science Foundation of China(Grant No.11374266),the Program for New Century Excellent Talents in University,China,and the Natural Science Foundation of Zhejiang Province,China(Grant No.Z15A050001).

?Corresponding author.E-mail:gaoxl@zjnu.edu.cn