聚噻吩單鏈量子熱輸運的第一性原理研究?

吳宇蔡紹洪鄧明森3)4)孫光宇2)劉文江

1)(貴州大學大數據與信息工程學院,貴陽 550025)

2)(貴州師范學院物理與電子科學學院,應用物理研究所,貴陽 550018)

3)(貴州師范學院,貴州省納米材料模擬與計算重點實驗室,貴陽 550018)

4)(貴州財經大學,貴州省經濟系統仿真重點實驗室,貴陽 550025)

1 引 言

有機導電材料聚噻吩自問世以來,以奇異的電學和光學性質吸引了眾多研究人員的關注.2000年,諾貝爾化學獎授予艾倫·黑格、艾倫·麥克德爾米德和白川英樹,以表彰他們在有機導電材料研究領域中的開創性貢獻.2014年,世界上最小的發光二極管由聚噻吩單分子鏈實現,該材料的神奇用途得以展示[1],另外,在太陽能電池研究領域也能見到聚噻吩活躍的身影[2].

聚噻吩塊體通常被看作絕熱材料,熱導率低于1 W·m?1·K?1,其絕熱性由塊體內部高分子鏈的形變和大量界面狀態導致.事實上,沿聚噻吩分子鏈方向排列的無定形聚噻吩納米纖維,室溫下熱導率高達4.4 W·m?1·K?1,高于聚噻吩塊體的熱導率[3].這種奇特的微觀熱輸運機理值得探索.此外,聚噻吩作為熱電材料也受到人們的青睞[4].根據ZT=S2GT/(Kph+Ke)[5],其中ZT為熱電優值,表征材料的熱電轉換效率,S為材料的塞貝克系數,T為絕對溫度,G為電導,Kph為聲子熱導,Ke為電子熱導,由于本征聚噻吩單鏈的能隙約為2.0—2.3 eV[6],故Ke可忽略不計.探討同位素摻雜對聚噻吩單鏈聲子熱導Kph的抑制效果,提高聚噻吩材料的熱電轉換效率,也是一個頗具價值的問題.這些都與聚噻吩單鏈的熱輸運性質緊密相關.我們希望確定純聚噻吩單鏈熱導率的理論上限,這在實驗上難以實現,因為微觀上存在諸多復雜且難以控制的因素.其次,同位素雜質是影響材料熱輸運性質的重要因素.石墨烯的熱輸運及其同位素效應在理論[7?12]和實驗方面都有了較為深入的研究.溫度約為室溫時,12C石墨烯帶中摻入一半13C原子后,熱導率幾乎降低一半[13].室溫下,將自然豐度BN納米管(19.9%10B,80.1%11B)中的11B含量提升至99.5%后,熱導率提高了近一倍[14].

對于另一種有機高分子聚合物——聚乙烯,其出色的導熱性能在實驗和第一性原理角度已有研究[15,16].而聚噻吩單鏈量子熱輸運目前的理論研究僅限于分子動力學(MD)方法[3,17].為了體現該微觀輸運現象的量子特征,從量子力學出發,同時可有效避免MD中計算方法或經驗勢函數不同所導致的計算結果差異.鑒于聚噻吩在有機發光器件、太陽能電池、熱電材料等方面的廣泛用途,研究聚噻吩單鏈的量子熱輸運及其同位素效應,對認識聚噻吩的導熱潛力以及調控其熱性能十分必要.本文在密度泛函理論(DFT)計算基礎上,用中間插入延展(CIS)[18,19]方法結合非平衡格林函數(NEGF)[20,21]方法,考察純聚噻吩單鏈在室溫下熱導率的理論上限,以及一般摻雜比例下的熱導同位素效應.

2 計算模型與方法

計算模型中,散射區域包含32個聚噻吩單胞,共計448個原子,總長度L=25.107 nm,兩端分別為半無限長聚噻吩單鏈熱極.納米尺度下,聲子的彈性輸運占主導因素,若忽略聲子間以及電、聲子間的相互作用,且考慮特定摻雜比例下同位素雜質隨機構型不同引起的透射系數平均效應,Landauer公式[12,22,23]給出的平均熱導Kph為

式中?為約化普朗克常數,ω為聲子頻率,kB為玻爾茲曼常數,T為體系的絕對溫度,Tph(ω)為聲子透射系數,

式中I為單位矩陣,D為散射區域的動力學矩陣:

式中M為原子質量,力常數二階張量Kiα,jβ的定義為

式中i,j為原子編號,α,β={x,y,z}表示三個正交方向,E為微觀體系的總能量,R0為體系能量最低狀態下的穩定幾何構型,Qiα是第i個原子α方向的小位移量,F是原子間的相互作用力.

圖1 含隨機同位素雜質13C的聚噻吩單鏈示意圖Fig.1.Illustration of a polythiophene chain containing random con fi guration of13C isotopes.

為了體現同位素隨機摻雜后透射系數的統計平均效應,對每一種確定的同位素摻雜比例(100%摻雜的情形除外,此時物質為純聚噻吩單鏈,沒有隨機因素,體系的平均透射系數與透射系數相同,平均熱導與熱導相同)計算10個隨機同位素雜質構型,逐一計算各隨機構型下的透射系數,再求其算術平均值,進而獲得該摻雜比例下的平均透射系數〈Tph(ω)〉以及平均熱導Kph.

聚噻吩單鏈的優化和力常數二階張量的計算采用量子化學軟件Gaussian09.計算選用雜化密度泛函b3lyp及6-31g(d)基組,優化后沿聚噻吩單鏈的方向單胞的晶格常數為7.846 ?(1 ?=0.1 nm),與X射線衍射(XRD)實驗[24]測得的聚噻吩薄膜的晶格常數7.56 ?相近.構建的超胞包含3個聚噻吩單胞,在周期性邊界條件下進行優化后,可得到聚噻吩單鏈的穩定幾何結構,此時各原子的最大受力數值都低于10?3eV/?,在此基礎上,計算相應的力常數張量,并使最終得到的力常數張量與原物理體系的對稱性一致[25,26].

在DFT計算獲取超胞力常數張量的前提下,通過CIS方法進一步構造了25.107 nm長的聚噻吩單鏈的動力學矩陣.鑒于DFT計算與實驗結果的差異,最終的聲子透射譜經過計算化學比較和基準數據庫(CCCBDB)的頻率校正,校正因子取為0.960[27].CIS方法是一種電子結構線性標度算法,在計算大型準周期電子體系結構和量子輸運中已得到嚴格驗證和廣泛應用[28?30].利用CIS方法擴展體系的力常數張量,極大地提高了計算效率,可以處理長達25.107 nm的聚噻吩單鏈及其同位素摻雜的聲子散射問題.采用CIS方法得到的純12C,32S聚噻吩單鏈的聲子透射譜和熱導,與直接用超胞力常數張量計算的結果完全相同,這印證了CIS方法計算結果的正確性.

聲子在傳播過程中與同位素雜質發生散射,使體系熱導減小.納米體系中,若直接用第一性原理方法計算包含數百原子體系的力常數張量,計算量極為巨大.作為簡化處理,討論聲子與雜質散射常見的方法有級聯散射模型[31,32]和標度理論方法[33,34].二者假定在低摻雜濃度下,雜質對聲子的散射總效應可看作一系列相繼單道散射效果的疊加,但這一假設對同位素雜質濃度(原子百分數,下同)高于10%的情況并不成立[12,31].用CIS方法在保持第一性原理精度的條件下大大提高了計算效率,因而可以采用直接模擬同位素雜質隨機分布的辦法,對包含448個原子的聚噻吩單鏈在各同位素摻雜比例下的平均熱導進行研究.該方法可體現單、多道聲子散射效應同時存在的一般情形,理論上能比較準確完整地描述一般同位素摻雜比例情況下的平均熱導變化.

C的穩定同位素有12C,13C兩種,其中12C的豐度為98.93%;S的穩定同位素有32S,33S,34S,36S,其中32S的豐度為94.99%.根據聚噻吩C,S原子的自然豐度,研究沒有同位素雜質時純聚噻吩體系的熱導,即先考查純凈聚噻吩單鏈在以12C,32S作為基本組分時的熱輸運規律,再研究單獨摻入13C,36S雜質時熱導的同位素效應,最后討論聚噻吩單鏈熱導的原子質量效應.

3 結果與討論

圖2中純凈12C,32S聚噻吩單鏈聲子的透射譜呈明顯的量子化特征,熱導量子化的基本單元為[35](h為普朗克常數),從第一個顯著的量子化臺階可以看出,體系在低頻極限時透射系數Tph=3.圖2中最高聲子頻率在1500 cm?1左右,該頻段對應C—C鍵伸縮振動,這與聚噻吩薄膜晶體拉曼光譜的測量結果一致[36].

為了考察不同頻率的聲子對熱導的貢獻,以圖2中純凈12C,32S聚噻吩單鏈聲子的透射譜為例,將(1)式的被積函數在T=300 K時作無量綱化處理,除以其在積分區間內的最大值,得到不同頻率聲子對熱導貢獻的相對強度r(ω).由圖3可知,低頻段聲子的透射能力對體系的熱導有決定性的貢獻.

圖2 100%12C和25%13C同位素摻雜的平均透射譜Fig.2.Average transmittance of pure12C sample and sample with13C isotope concentration of 25%.

圖3 T=300 K時不同頻率聲子對熱導貢獻的相對強度Fig.3.Relative intensity of a phonon contributing to the thermal conductance at different frequencies at 300 K.

如圖4所示,繼續提高13C原子的摻雜濃度,當C原子全為13C時,體系重新成為純凈的聚噻吩單鏈,此時透射譜再次呈現整齊的臺階特征,這說明同位素雜質概念的相對性.該透射譜形狀與純凈12C聚噻吩的透射譜非常相似,只是整個透射譜向低頻方向移動,截止頻率從1550 cm?1減少到1500 cm?1,相應的C—C鍵振動頻率降低.

從圖5可以看出,隨著溫度的升高,體系的平均熱導單調增加.對于12C聚噻吩單鏈,隨著13C摻雜濃度由低到高,體系的平均熱導先減小后變大,其中13C摻雜濃度為25%與75%時兩條曲線基本重合.純12C聚噻吩單鏈的熱導比純13C的聚噻吩單鏈略大.其他條件一定,聚噻吩單鏈中12C,13C原子等比例摻雜時,平均熱導總是最小.

圖4 50%13C同位素摻雜和100%13C的平均透射譜Fig.4.Average transmittance of samples with13C isotope concentrations of 50%and 100%.

圖5 聚噻吩單鏈各同位素摻雜比例下的平均熱導溫度曲線Fig.5.Average thermal conductance versus temperature for different isotopic compositions.

考慮由12C,32S組成的純凈聚噻吩單鏈的熱導率.在聲子平均自由程以內彈性輸運條件下,聚噻吩單鏈的熱導率σ與散射區域的長度L成正比,σ=KphL/A[16],其中A為聚噻吩單鏈的橫截面積.根據聚噻吩晶體薄膜X射線衍射實驗的結果[24],用聚噻吩晶體薄膜的晶格常數a=10.80 ?,b=4.74 ?的乘積來定義橫截面積.實驗測定聚噻吩單鏈中的聲子平均自由程非常困難[3].從實驗結果來看,沿聚噻吩分子鏈方向排列但處于無定形狀態的聚噻吩纖維熱導率能夠達到4.4 W·m?1·K?1[3],可以推測純凈的一維聚噻吩單鏈晶體當有更為可觀的量子熱輸運能力,12C,32S組成的純聚噻吩單鏈中的聲子平均自由程大于32 nm.為了便于比較,取L=32 nm,得到T=300 K時32 nm長的純聚噻吩單鏈的熱導率為30.2 W·m?1·K?1(圖6),約為室溫下相同長度聚乙烯單鏈熱導率的1/3[16],而與鉛的熱導率35 W·m?1·K?1相近. 文獻[3]運用平衡分子動力學(EMD),在周期性邊界條件下,選用ReaxFF勢函數計算得到室溫下32 nm長聚噻吩單鏈的熱導率為43.3 W·m?1·K?1.另外,文獻[17]在MD框架下,同樣選取ReaxFF勢函數,用格林久保模態分析(GKMA)方法,研究了聚噻吩單鏈熱導率與鏈長的關系,結果表明室溫下32 nm長的聚噻吩單鏈熱導率約7.5 W·m?1·K?1,且在該輸運尺度以內,熱導率與鏈長成正比;120 nm長的聚噻吩單鏈熱導率超過120 W·m?1·K?1.該結論有力支持了純聚噻吩單鏈中聲子平均自由程大于32 nm推測的合理性.從量子力學方法出發可得,室溫下32 nm長聚噻吩單鏈熱導率為30.2 W·m?1·K?1,該結果恰好處于EMD和GKMA方法給出的熱導率之間(43.3 W·m?1·K?1＞30.2 W·m?1·K?1＞7.5 W·m?1·K?1).我們認為EMD的計算結果偏高是采用經典力學近似和選用的經驗勢函數ReaxFF導致的.GKMA方法在ReaxFF勢模擬基礎上考慮了聲子之間的相互作用,體現出聲子散射的非諧效應,故而得出的熱導率7.5 W·m?1·K?1小于我們在彈性輸運條件下用量子力學方法計算得到的結果.上述三個獨立的理論研究結果一致表明,絕熱的聚噻吩材料的熱輸運潛力十分可觀,仍有相當大的提升空間.

T=300 K時12C聚噻吩鏈中由低到高摻入不同濃度的13C同位素雜質時,體系平均熱導變化規律呈U型曲線.這與石墨烯帶和硅納米線熱輸運的同位素效應規律相似[13,37,38].作為對比,我們計算了純32S聚噻吩單鏈中摻入同樣比例36S同位素雜質的情形,體系平均熱導變化的規律相近,但變化量相對較小.

從圖7還可以發現,聚噻吩12C,13C熱導的同位素效應在等比例摻雜時最為突出,平均熱導比純凈的12C聚噻吩單鏈降低了31%;32S,36S的情況也是如此,只是其平均熱導比未摻雜時降低了9%.12C,13C同位素原子各以50%比例隨機摻雜時對聲子的散射最強烈,此時平均熱導最小.聲子與同位素雜質作用時單、多道散射同時存在,且多道散射的影響具有決定性作用,體現出整體的散射行為.誤差棒的長度特征表明,同樣摻雜比例下,C原子同位素效應引起聚噻吩單鏈的熱導分散性比S原子更加突出.

圖6 T=300 K時32 nm長純12C,32S聚噻吩單鏈熱導率Fig.6.Thermal conductivity of pure12C and32S polythiophene chain with a length of 32 nm at 300 K.

圖7 T=300 K時不同13C,36S雜質摻雜比例的聚噻吩單鏈平均熱導Fig.7.Average thermal conductance of polythiophene doped with various13C and36S concentrations at 300 K.

從結構化學的角度來看,C=C鍵的鍵能(約600 kJ/mol)比C—S鍵的鍵能(約270 kJ/mol)大,即相鄰C,C原子間的耦合比C,S原子間的耦合更為緊密,而S的原子質量約為C原子質量的3倍,故S原子的存在將抑制聲子沿聚噻吩單鏈中C—S—C結構的傳播,因而聚噻吩單鏈中的聲子主要從相鄰C原子為骨架的主鏈中傳輸;另外,一個聚噻吩單胞中C原子數與S原子數的比例為4:1,這就意味著同一摻雜濃度下,聚噻吩單鏈中聲子與C原子同位素雜質散射的次數將是S原子同位素雜質數目的4倍.因此,對于同樣的隨機摻雜濃度,聚噻吩單鏈中C原子同位素效應導致的體系熱導削弱比S原子顯著得多.可以推測,同樣摻雜比例下,C=C主鏈上出現C同位素雜質對體系熱導的擾動幅度比非C=C主鏈上出現S原子同位素雜質時更為明顯,事實上,圖7中,在相同的摻雜比例下,C,S同位素雜質引起體系熱導波動的誤差棒長度特征有力地印證了這一點.

聚噻吩熱導的同位素效應可以作為熱導調制的一個有效途徑.從散熱角度考慮,希望盡量提高聚噻吩的熱導率,這對設計諸如聚噻吩單分子發光二極管之類的納米電子器件,尤其是防止局域熱導致的熔斷,至關重要;而熱電材料的開發卻需要抑制聚噻吩熱導以提高熱電轉換效率,這可以通過調節C同位素雜質的最佳摻雜濃度實現.理論預測可知C元素的熱導同位素效應比S元素更顯著,因而調節效果更佳.

圖8所示為T=300 K時,聚噻吩單鏈中保持32S不變,碳元素分別為12C,13C,14C的質量熱導效應曲線;以及保持12C不變,S元素分別為32S,33S,34S,36S的質量熱導效應曲線.從圖8可以看出,隨著C的相對原子質量變大,熱導基本呈反比例減小;精確的計算表明,隨著S的相對原子質量變大,熱導反而緩慢增加,后逐漸趨于恒定.這說明相對原子質量熱導效應曲線隨著微觀體系結構的不同可以有復雜的表現形式,由整個體系的質量分布和力常數張量共同決定.

圖8 T=300 K時純聚噻吩單鏈熱導隨相對原子質量的變化Fig.8.Thermal conductance of pure polythiophene chain versus relative atom mass at 300 K.

4 結 論

應用DFT,CIS方法結合NEGF方法,預測了室溫下32 nm長純聚噻吩單鏈熱導率的理論上限為30.2 W·m?1·K?1,并與MD方法的結果進行了比較.另外,對長度為25.107 nm、包含448個原子的聚噻吩單鏈在一般摻雜比例下的熱導同位素效應進行了研究,克服了級聯散射模型及標度理論僅適用于低摻雜濃度(≤10%)的不足.計算發現,同等摻雜比例下,聚噻吩中C元素熱導的同位素效應比S元素更為顯著.室溫下,12C,13C隨機同位素摻雜的聚噻吩單鏈,摻雜比例為1:1時熱導同位素效應最為顯著,此時的平均熱導比未摻雜前至少降低了30%.這對于認識聚噻吩單鏈熱導的調控因素,挖掘聚噻吩材料在熱輸運方面的潛力,優化其在有機發光器件、太陽能電池和熱電材料方面的性能具有積極的意義.同時,探索準一維高分子鏈體系熱導同位素效應的特點,對于調制納米線的熱導[39,40]及提高納米材料的熱電優值[41,42],也有一定的參考借鑒價值.

[1]Reecht G,Scheurer F,Speisser V,Dappe Y J,Mathevet F,Schull G 2014Phys.Rev.Lett.112 047403

[2]Bulumulla C,Du J,Washington K E,Kularatne R N,Nguyen H Q,Michael C B,Stefan M C 2017J.Mater.Chem.A5 2473

[3]Singh V,Bougher T L,Weathers A,Singh V,Bougher T L,Weathers A,Cai Y,Bi K,Pettes M T,McMenamin S A,Lv W,Resler D P,Gattuso T R,Altman D H,Sandhage K H,Shi L,Henry A,Cola B A 2014Nature Nanotech.9 384

[4]Cowen L M,Atoyo J,Carnie M J,Baran D,Schroeder B C 2017ECS J.Solid State Sci.Technol.6 3080

[5]Chen X B,Duan W H 2015Acta Phys.Sin.64 186302(in Chinese)[陳曉彬，段文暉 2015物理學報 64 186302]

[6]Bouzzine S M,Salgado-Morán G,Hamidi M,Bouachrine M,Pacheco A G,Glossman-Mitnik D 2015J.Chem.2015 296386

[7]Tan Z W,Wang J S,Chee K G 2011Nano Lett.11 214

[8]Xu Y,Chen X B,Gu B L,Duan W H 2009Appl.Phys.Lett.95 233116

[9]Xie Z X,Tang L M,Pan C N,Li K M,Chen K Q,Duan W H 2012Appl.Phys.Lett.100 073105

[10]Ouyang T,Chen Y P,Xie Y,Wei X L,Yang K K,Yang P,Zhong J X 2010Phys.Rev.B82 245403

[11]Zhang H J,Lee G,Fonseca A F,Borders T L,Cho K 2010J.Nanomater.7 537657

[12]Sevin?li H,Sevik C, ?a?n T,Cuniberti G 2013Nature.Sci.Rep.3 1228

[13]Chen S S,Wu Q Z,Mishra C,Kang J Y,Zhang H J,Cho K,Cai W W,Balandin A A,RuoffR S 2012Nature Mater.11 203

[14]Chang C W,Fennimore A M,Afanasiev A,Okawa D,Ikuno T,Garcia H,Li D Y,Majumdar A,Zettl A 2006Phys.Rev.Lett.97 085901

[15]Shen S,Henry A,Tong J,Zheng R T,Chen G 2010Nature Nanotech.5 251

[16]Jiang J W,Zhao J H,Zhou K,Rabczuk T 2012J.Appl.Phys.111 124304

[17]Lv W,Winters M,Deangelis F,Weinberg G,Henry A 2017J.Phys.Chem.A121 5586

[18]Gao B,Jiang J,Liu K,Wu Z Y,Lu W,Luo Y 2007J.Comput.Chem.29 434

[19]Jiang J,Liu K,Lu W,Luo Y 2006J.Chem.Phys.124 214711

[20]Taylor J,Guo H,Wang J 2001Phys.Rev.B63 245407

[21]Wang J S,Wang J,Lü J T 2008Eur.Phys.J.B62 381

[22]Yamamoto T,Watanabe S,Watanabe K 2004Phys.Rev.Lett.92 075502

[23]Mingo N,Yang L 2003Phys.Rev.B68 245406

[24]Satoh M,Yamasaki H,Aoki S,Yoshino K 1988Mol.Cryst.Liq.Cryst.Inc.Nonlinear Opt.159 289

[25]Mingo N,Stewart D A,Broido D A,Srivastava D 2008Phys.Rev.B77 033418

[26]Nikoli? B K,Saha K K,Markussen T,Thygesen K S 2012J.Comput.Electron.11 78

[27]NIST ComputationalChemistry Comparison and Benchmark Database http://cccbdb.nist.gov/vibscalejust.asp

[28]Hu W P,Jiang J,Nakashima H,Luo Y,Kashimura Y,Chen K Q,Shuai Z,Furukawa K,Lu W,Liu Y Q,Zhu D B,Torimitsu K 2006Phys.Rev.Lett.96 027801

[29]Jiang J,Gao B,Han T T,Fu Y 2009Appl.Phys.Lett.94 092110

[30]Jiang J,Sun L,Gao B,Wu Z Y,Lu W,Yang J L,Luo Y 2010J.Appl.Phys.108 094303

[31]Savic I,Mingo N,Stewart D A 2008Phys.Rev.Lett.101 165502

[32]Stewart D A,Savic I,Mingo N 2009Nano Lett.9 81

[33]Markussen T,Jauho A P,Brandbyge M 2009Phys.Rev.B79 035415

[34]Markussen T,Rurali R,Jauho A P,Brandbyge M 2007Phys.Rev.Lett.99 076803

[35]Rego L G C,Kirczenow G 1998Phys.Rev.Lett.81 232

[36]Fu M X,Shi G Q,Chen F G,Hong X Y 2002Phys.Chem.Chem.Phys.4 2685

[37]Jiang J W,Lan J H,Wang J S,Li B W 2010J.Appl.Phys.107 054314

[38]Yang N,Zhang G,Li B W 2008Nano Lett.8 276

[39]Hu M,Giapis K P,Goicochea J V,Zhang X,Poulikakos D 2011Nano Lett.11 618

[40]Liu Y Y,Zhou W X,Tang L M,Chen K Q 2014Appl.Phys.Lett.105 203111

[41]Zhou W X,Chen K Q 2014Nature.Sci.Rep.4 7150

[42]Zhou W X,Chen K Q 2015Carbon85 24