陳 星， 劉 勛， 黃海軍

( 武漢理工大學 理學院 高壓物理和新材料研究中心 武漢 430070 )

1 引 言

地球形成過程會在地核及地幔中留下地球物理和地球化學方面的痕跡. 地震波觀測表明，地核是由鐵及少量輕元素組成. 確定地核中輕元素的種類和含量有助于構(gòu)建早期地球形成的模型，解釋地球物理現(xiàn)象以及預測地球未來的發(fā)展[1-5]. 由于碳元素在核-幔分離條件下對液態(tài)鐵具有高親和性，并且在太陽系和C1碳質(zhì)球粒隕石中含量較高，被認為是地核中重要的輕元素之一[6, 7]. 地球內(nèi)核具有異常高的泊松比，在眾多候選輕元素(碳，氫，氧，硫和硅)中，只有鐵碳化合物的泊松比接近內(nèi)地核[8, 9]. 因此，碳可能是內(nèi)地核中主要的輕元素，甚至一些研究表明內(nèi)地核幾乎完全是由鐵碳化合物組成[9]. 通常認為內(nèi)地核中鐵碳化合物的存在形式是Fe3C和Fe7C3，目前這兩種鐵碳化合物都已得到廣泛的研究[6, 10-16]. 雖然大部分研究認為在內(nèi)地核條件下Fe7C3比Fe3C更穩(wěn)定，但在最新的研究中Takahashi等人通過實驗發(fā)現(xiàn)Fe3C在內(nèi)地核條件下可以穩(wěn)定存在[17].

約束內(nèi)地核中碳含量最直接的方法是在地核條件下測量鐵碳化合物的密度和聲速，并將其與地震波模型(如初始地球參考模型-PREM)進行比較. 在內(nèi)地核條件(330-360 GPa，5000-7000 K)下直接測量物質(zhì)的性質(zhì)仍然非常困難. 目前對Fe3C狀態(tài)方程和聲速的測量局限于相對較低的壓力和溫度條件. 實驗測量的Fe3C聲速數(shù)據(jù)即使在室溫高壓條件下也顯示出很大的差異，將這些結(jié)果外推至內(nèi)地核條件會帶來很大的不確定性. 另一方面，基于第一性原理的模擬計算沒有上述限制. 本文通過第一性原理GGA-PBE方法計算了0-K下Fe3C的狀態(tài)方程和聲速，并將這些結(jié)果外推至內(nèi)地核溫度條件，通過與PREM模型比較以限定內(nèi)地核中的碳含量.

2 計算方法

常態(tài)下，F(xiàn)e3C為正交結(jié)構(gòu)(圖1)，空間群為Pnma，其晶格常數(shù)為a=5.077 ?，b=6.724 ?以及c=4.515 ?[18]. 本文采取GGA-PBE方法進行贗勢平面波計算[19]，計算選取的截斷能為500 eV，Monkhorst-Pack網(wǎng)格為6×4×7，幾何優(yōu)化算法為BFGS[21]，采用超軟贗勢來描述離子-電子間的相互作用[20]. 模擬晶胞包含12個Fe原子和4個C原子，計算包含F(xiàn)e 3d64s2和C 2s22p2電子. Fe3C在高壓下經(jīng)歷了從鐵磁性相(FM)到非磁性相(NM)的同構(gòu)型轉(zhuǎn)變，文獻中報道的相變壓力還存在爭議(4.3-55 GPa). 為了研究Fe3C在高壓下的磁性相變，我們對其進行了自旋限制和自旋非限制計算. 根據(jù)胡克定律計算了非磁性相Fe3C的彈性常數(shù)并進而得到其聲速. 我們也對不同壓力下的電子態(tài)密度進行了計算以獲得其對熱壓的貢獻.

圖1 Fe3C在常態(tài)下的晶體結(jié)構(gòu). 黑色和紅色的球分別代表C和Fe原子，每個C原子被六個Fe原子包圍形成一個八面體，此處僅顯示一個八面體. Fig. 1 Crystal structure of Fe3C at ambient conditions. Black and red balls represent C and Fe atoms, respectively. Avery C atom is surrounded by six Fe atoms to form an octahedron, only one octahedron is shown here.

3 計算結(jié)果

圖2顯示了有和沒有自旋限制的情況下，計算所得Fe3C的0-K等溫線. 結(jié)果顯示沒有自旋限制的等溫線在40 GPa處出現(xiàn)不連續(xù)，表明存在相變，相變處的密度增加了約2%，壓力超過40 GPa后，這兩者趨于一致. 計算的 Fe3C隨壓力變化的磁矩繪制在圖3中，當壓力小于40 GPa時，磁矩隨著壓力的增加而緩慢減小，當壓力大于40 GPa時，磁矩突然下降至接近0，表明在40 GPa時Fe3C由鐵磁相轉(zhuǎn)變?yōu)榉谴判韵? 在靜高壓實驗中，F(xiàn)e3C的自旋相變已獲得廣泛的研究. 使用X射線發(fā)射譜，Lin等人報道Fe3C在25 GPa存在磁坍塌[22]. 基于同步輻射穆斯保爾譜測量，Gao等人報道Fe3C在4.3-6.5 GPa出現(xiàn)磁性相變[23]. 根據(jù)X射線磁圓二色性測量，Acet等人報道Fe3C在9 GPa附近存在磁性轉(zhuǎn)變[24]. Ono等人通過X射線衍射技術(shù)測量了0-70 GPa壓力范圍內(nèi)Fe3C的狀態(tài)方程，他們發(fā)現(xiàn)在55 GPa處存在1.7%的體積塌縮并將其解釋為磁性相變[25]. 最近，Chen等人通過X射線發(fā)射譜重新檢測了Fe3C在高壓下的磁性性質(zhì)，他們的研究結(jié)果表明Fe3C在10-50 GPa之間會逐步發(fā)生鐵磁性到順磁性，再到非磁性的轉(zhuǎn)變[8]，我們的計算結(jié)果與這個結(jié)論吻合. 計算得到的0-K等溫線與靜高壓實驗獲得的300-K等溫線數(shù)據(jù)也基本一致[16, 25-28]. 對鐵磁性和非磁性相，使用三階Birch-Murnaghan狀態(tài)方程對計算結(jié)果進行擬合可得K0= 235 GPa，K0′= 4.0(鐵磁相)和K0= 332 GPa，K0′= 4.6 (非磁性相)，其中K0為體積模量，K0′ 為其壓力偏導數(shù). 在表I中，我們將擬合獲得的狀態(tài)方程參數(shù)與實驗結(jié)果進行了比較，對于鐵磁性相，計算的體積模量略高于實驗結(jié)果. 對于非磁性相，所得結(jié)果與Takahashi等人的最新結(jié)果非常吻合[17].

圖2 Fe3C的靜態(tài)等溫線具有(藍色實心圓)和沒有(紅色實心圓)自旋限制Fig. 2 Static isothermals of Fe3C calculated with (blue solid circles) and without (red solid circles) spin restriction

圖3 高壓下Fe3C的磁矩Fig. 3 Magnetic moments of Fe3C at high pressure

計算所得的彈性常數(shù)如圖4(a)所示. 因為我們關(guān)注的是Fe3C在內(nèi)地核壓力下的性質(zhì)，因而只計算了非磁性相Fe3C的彈性常數(shù). 所有的彈性常數(shù)都隨壓力升高而線性增大. 發(fā)生彈性應變的晶體只有在其內(nèi)能增加的情況下才保持穩(wěn)定. 對于正交晶體如Fe3C，該穩(wěn)定性判據(jù)表現(xiàn)為[29]：

表1 Fe3C的狀態(tài)方程參數(shù)

C22+C33-2C23>0,

(1)

C11+C22+C33+2C12+2C13+2C23>0,

(2)

C11>0,C22>0,C33>0,

C44>0,C55>0,C66>0.

(3)

計算所得的C44在零壓下為負值，表明非磁性相Fe3C在零壓下不穩(wěn)定. 但隨著壓力的增加，非磁性相穩(wěn)定性增加，C44也變?yōu)檎? (C22+C33- 2C23)為正值，并且隨著壓力升高而增加，如圖4(b)所示，表明非磁性相Fe3C在內(nèi)地核壓力下保持機械穩(wěn)定.

圖4 (a)NM相Fe3C在高壓下的彈性常數(shù)；(b)(C22 + C33-2C23)的值與壓力的關(guān)系Fig. 4 (a) Elastic constants of NM phase Fe3C at high pressure; (b) The value of C22+C33-2C23 versus pressure

為了獲得電子的熱壓貢獻，我們在0 ～ 360 GPa壓力范圍內(nèi)計算了非磁性相Fe3C的電子態(tài)密度，壓力間隔為40 GPa，計算結(jié)果如圖5(a)所示. 不同壓力下的電子總態(tài)密度相似，費米能級處的電子態(tài)密度不為零，表明Fe3C在整個壓力范圍內(nèi)為金屬，因而其在極端條件下的電子熱壓貢獻不可忽略. 圖5(b)所示為P= 120 GPa時的部分電子態(tài)密度，可見主要貢獻來自于鐵的3d電子.

圖5 (a)最高360 GPa，間隔為40 GPa的Fe3C的總電子態(tài)密度；(b)P = 120 GPa時的典型局部電子態(tài)密度，黑色虛線表示費米能級Fig. 5 (a) Total electronic DOSs of Fe3C up to 360 GPa with an interval of 40 GPa; (b) Typical partial DOS at P=120 GPa. Black dashed lines indicate Fermi level.

4 熱力學計算和地球物理意義

地球內(nèi)核處于高壓高溫狀態(tài)，為了約束地核中的碳含量，必須將計算所得的0-K等溫線外推至高溫條件. 我們將熱壓寫為：

(4)

其中V=1/r為比容，g為格林愛森參數(shù). 晶格熱容Cvib=3R/m，其中R為理想氣體常數(shù)，m為摩爾質(zhì)量. 電子熱容Ce=β0(V/V0)γeT，其中b0為電子比熱系數(shù)，ge為電子格林愛森參數(shù)[30]. 有效格林愛森參數(shù)可表述為[31]：

γ=(γvibCvib+γeCe)/CV,

(5)

其中g(shù)vib=gvib,0(V/V0)q為晶格對格林愛森參數(shù)的貢獻，gvib,0是零壓和零溫下的晶格格林愛森參數(shù)，q為常數(shù). 由上述表達式可以得到：

(6)

上式中有四個參數(shù)gvib,0，q，b0和ge需要確定. 一般需要通過計算聲子譜來確定gvib，這里我們使用德拜格林愛森參數(shù)gD來代替，它與德拜聲速vD的關(guān)系如下[32]：

(7)

圖6 計算所得高壓下非磁性相Fe3C的德拜Grüneison參數(shù)gDFig.6 Calculated Debye Grüneison parameter gD of NM phase Fe3C at high pressure

電子格林愛森參數(shù)ge和比熱系數(shù)b0可以通過計算出的電子態(tài)密度來確定. 電子內(nèi)能為：

(8)

其中D(e) 為電子態(tài)密度，f(e) 為費米分布函數(shù)[34]. 如果電子態(tài)密度在費米能級附近沒有快速變化，將Sommerfeld展開應用于(8)式，則可得電子比熱為：

(9)

現(xiàn)根據(jù)(6)式和上面確定的熱力學參數(shù)，將Fe3C的0-K等溫線外推至內(nèi)地核條件. 內(nèi)地核的溫度剖面由Tg=TICB(rI/rICB)1.5確定，其中rI是地核密度，TICB和rICB分別是內(nèi)-外地核邊界處的溫度和密度. 根據(jù)Hirose的討論，取TICB=5400 K，計算結(jié)果如圖7所示，相同壓力下Fe3C的密度比PREM值低約0.67 g/cm3[37]. 我們也根據(jù)Fei等人報道的狀態(tài)方程參數(shù)計算了內(nèi)地核條件下純鐵的密度[38]. 如果Fe和Fe3C為理想混合，計算結(jié)果表明需要47 wt%的Fe3C來解釋內(nèi)地核的密度虧損，這意味著如果碳是內(nèi)核中唯一的輕元素且以Fe3C的形式存在，則內(nèi)地核含有最多3.1wt%的碳. Murphy等人通過實驗測量發(fā)現(xiàn)鐵的gvib比gD大約10%，我們通過計算發(fā)現(xiàn)晶格格林愛森參數(shù)的微小變化不會影響上述討論結(jié)果.

圖7 內(nèi)地核條件下Fe3C(紅線)和純鐵(黑線)的密度. 藍線為根據(jù)理想混合計算的Fe+47wt%Fe3C的密度，黑色十字為PREM模型. Fig.7 Densities of Fe3C (red line) and pure iron (black line) at the inner core conditions. The blue lines are sound velocities of pure iron and black crosses are from PREM.

地核中適當?shù)妮p元素應能夠同時解釋純鐵和地震波觀測值之間的密度和聲速差異. 實驗測量的Fe3C在高壓下的聲速數(shù)據(jù)存在很大差異，如圖8所示[8, 23, 39, 40]，將這些結(jié)果外推到內(nèi)地核條件將得出截然不同的結(jié)論. 壓力較低時，我們計算的結(jié)果與實驗結(jié)果非常吻合，但是由于非磁性相Fe3C的不穩(wěn)定性導致零壓下的Vs明顯較小. 在較高的壓力下，Vp和Vs均高于Chen等人的實驗結(jié)果[8]. 在內(nèi)地核壓力下，F(xiàn)e3C的Vp高于純鐵，且Vs與純鐵相當[41]. 地核中的輕元素應降低純鐵的Vp和Vs使其更接近PREM值，碳因而不是一個好的選擇. 從目前的研究結(jié)果看，溫度對Fe3C聲速的影響還尚未可知. Gao等人在低壓下對Fe3C進行的實驗測量表明溫度對Vp的影響可忽略不計，但對Vs的影響較為明顯，他們推斷在內(nèi)地核的溫度條件下，Vs可能最多會下降30%[42]. 我們計算的Vs比PREM值高2.14 km / s，減小30%后仍然比PREM值高0.43 km / s，因此我們的計算結(jié)果表明碳不是內(nèi)地核中主要的輕元素.

圖8 計算得到的Fe3C聲速與實驗結(jié)果的比較(藍線是純鐵的聲速，黑色十字為PREM值)Fig.8 Comparison of the calculated Fe3C sound velocities with available experimental results (The blue lines are sound velocities of pure iron and black crosses are from PREM)

5 結(jié) 論

總結(jié)，我們通過第一性原理GGA-PBE方法計算了高壓下Fe3C的熱狀態(tài)方程和聲速. 計算結(jié)果表明Fe3C在40 GPa由鐵磁性相轉(zhuǎn)變?yōu)榉谴判韵? 將Fe3C的0-K等溫線外推至內(nèi)地核溫壓條件發(fā)現(xiàn)內(nèi)地核中碳含量最多為3.1wt%. 另一方面，計算得出的Fe3C聲速數(shù)據(jù)表明碳不應是內(nèi)地核中主要的輕元素.