馬超杰，吳 瀟，馬陽陽，何開華，姬廣富

（1.中國地質(zhì)大學(xué)（武漢）數(shù)學(xué)與物理學(xué)院，湖北 武漢 430074；2.中國工程物理研究院流體物理研究所沖擊波物理與爆轟物理重點實驗室，四川 綿陽 621999）

含碳固溶體的存在會影響地球內(nèi)部的物理和化學(xué)性質(zhì)[1-4]。近年來，菱鎂礦（MgCO3）被認為是碳進入地球深部的主要載體之一，因其在地球深部碳循環(huán)中的關(guān)鍵作用而引起廣泛關(guān)注[5-13]。Isshiki等[14]、Oganov等[15]的研究表明，菱鎂礦在下地幔的溫壓條件下能夠穩(wěn)定存在；Hazen等[13]、Oganov等[15]開展的高溫高壓實驗研究揭示，菱鎂礦和菱鐵礦的固溶體[(Mg,Fe)CO3]可以在低于100 GPa的壓力條件下穩(wěn)定存在。

鐵是多價態(tài)的過渡金屬，會對菱鎂礦和菱鐵礦固溶體的性質(zhì)產(chǎn)生非常重要的影響。此外，鐵的自旋在一定的壓力和溫度條件下可以發(fā)生轉(zhuǎn)變，引起物理性質(zhì)變化[16]。國內(nèi)外學(xué)者對鐵方鎂石(Mg,Fe)O的自旋轉(zhuǎn)變開展了大量研究，結(jié)果表明，鐵方鎂石中的鐵（Fe2+）在40～50 GPa范圍內(nèi)會從高自旋（High spin, HS）態(tài)向低自旋（Low spin, LS）態(tài)轉(zhuǎn)變，并伴隨著結(jié)構(gòu)、電子、光學(xué)、彈性和熱力學(xué)性質(zhì)的異常變化[16-28]。對于含鐵的菱鎂礦，已有的高壓穆斯堡爾光譜[29]、X射線發(fā)射光譜[30]、激光拉曼光譜[10]、X射線衍射[5,8,11]等實驗研究和第一性原理計算研究[31-32]都證實，(Mg,Fe)CO3中的鐵在約45 GPa時從HS態(tài)向LS態(tài)轉(zhuǎn)變，導(dǎo)致(Mg,Fe)CO3的體積比MgCO3減小6%～10%[11]。Liu等[11]、Hsu等[32]、Fu等[33]利用理論計算和實驗相結(jié)合的方法，詳細討論了彈性和地震波速在自旋轉(zhuǎn)變時的變化，得到了非常有意義的結(jié)果。其他熱力學(xué)參數(shù)如熱膨脹系數(shù)、格臨愛森常數(shù)和比熱容等在自旋轉(zhuǎn)變條件下的性質(zhì)尚缺少報道。

本研究利用第一性原理計算方法，開展高溫高壓下(Mg,Fe)CO3在HS、LS及HS和LS態(tài)共存時的混合自旋（Mixed spin, MS）態(tài)的熱力學(xué)性質(zhì)（熱膨脹系數(shù)、體變模量、體積、速度等熱力學(xué)參數(shù)）研究，并與不含鐵的MgCO3的相關(guān)性質(zhì)進行對比，分析引起(Mg,Fe)CO3熱力學(xué)性質(zhì)變化的機制。研究結(jié)果可為研究地幔深部碳的行為以及地幔在全球碳循環(huán)中的作用提供制約因素。

1 研究方法與理論

1.1 含鐵菱鎂礦的結(jié)構(gòu)和計算參數(shù)

菱鎂礦MgCO3屬于三角晶系，其空間群為，原胞中包含10個原子（2個Mg原子、2個C原子和6個O原子）。為了討論鐵及其在高溫高壓下的自旋轉(zhuǎn)變對菱鎂礦物理性質(zhì)的影響，本研究選擇用1個Fe原子替換1個Mg原子，得到Fe和Mg的摩爾比為1∶1，分子式為（Mg0.5Fe0.5）CO3。

幾何結(jié)構(gòu)優(yōu)化和相關(guān)能量的計算[34]采用基于密度泛函理論（Density of functional theory, DFT）的第一性原理分子動力學(xué)計算軟件VASP完成。由于結(jié)構(gòu)中有鐵存在，需要考慮強關(guān)聯(lián)作用，因此計算中考慮了Hubbard參數(shù)U的影響，即LDA+U。已有的研究表明，U值會隨自旋轉(zhuǎn)變而改變，本研究選取Tsuchiya等[18]、Krukau等[35]通過線性響應(yīng)理論計算得到的U值，分別為ULS= 5.3 eV，UHS= 4.0 eV。K點采用Monkhorst-Pack方法生成以Γ點為中心的15 × 15 × 15網(wǎng)格，截斷能設(shè)置為1 000 eV[36]。總能量的收斂閾值設(shè)置為1 × 10?6eV/cell，原子力的收斂閾值為1 × 10?3eV/cell。采用基于準諧近似的PHONOPY軟件計算得到二階原子間力常數(shù)（IFCs）、聲子譜和其他熱力學(xué)參數(shù)[37]。構(gòu)造求解高對稱q點聲子頻率本征值的動力學(xué)矩陣，然后利用倒空間中動力學(xué)矩陣的傅里葉變換，計算其他一般q點的聲子頻率。所有計算中，以Γ為中心的q點網(wǎng)格選取為30 × 30 × 30。計算熱力學(xué)參數(shù)之前，用VASP軟件完成不同位移或不同體積下結(jié)構(gòu)的能量計算，計算過程中的參數(shù)設(shè)置與上述幾何優(yōu)化參數(shù)設(shè)置相同。

1.2 MS態(tài)熱力學(xué)參數(shù)

根據(jù)已有的研究結(jié)果，MS態(tài)時的吉布斯自由能可表示為

式中：n為LS態(tài)在MS態(tài)中所占的百分比，p為壓力，T為溫度，GHS、GLS分別為HS態(tài)和LS態(tài)的吉布斯自由能，Gmix為HS-LS混合態(tài)的吉布斯自由能。Gmix可以表示為

式中：XFe為鐵在含鐵菱鎂礦中的摩爾分數(shù)，kB為玻爾茲曼常數(shù)。在給定的壓力和溫度條件下，通過求最小化MS態(tài)吉布斯自由能得到LS態(tài)的占比n

式中：ΔGLS?HS為(Mg0.5Fe0.5)CO3在LS態(tài)和HS態(tài)的吉布斯自由能之差。對于Fe2+來說，自旋和軌道簡并量子數(shù)分別為s= 2，m= 3。得到n后可以推導(dǎo)出MS態(tài)的熱膨脹系數(shù)α(n)和等溫體積模量KT(n)分別為

式中：KS為絕熱體積模量，CV,m為定容比熱容，ρ為密度。根據(jù)上述公式可以計算出MS態(tài)下的體積V、速度v和格臨愛森常數(shù)γ。

2 結(jié)果與討論

2.1 Fe2+的自旋轉(zhuǎn)變

含鐵菱鎂礦(Mg0.5Fe0.5)CO3的LS態(tài)與HS態(tài)的焓差ΔH（HHS-HLS）隨壓力的變化趨勢如圖1所示。由圖1可知：在低壓時，ΔH＜ 0，HS態(tài)的焓值小，因此 HS 態(tài)更穩(wěn)定；在高壓時，ΔH＞ 0，LS 態(tài)的焓值變小，LS態(tài)變得更穩(wěn)定。本研究計算得到的自旋轉(zhuǎn)變壓力為44.5 GPa，與實驗觀察得到的結(jié)果（40～52 GPa）很好地符合[5,10-11,29-30]，與 Hsu 等[32]的計算結(jié)果（48 GPa）也符合較好，較小的差異源于鐵在含鐵菱鎂礦中的摩爾分數(shù)不同。

圖1 (Mg0.5Fe0.5)CO3的LS態(tài)與HS態(tài)的焓差ΔH（HHS-HLS）隨壓力的變化Fig.1 Pressure dependence of the enthalpy difference (ΔH)between LS state and HS state (HHS-HLS） for (Mg0.5Fe0.5)CO3

2.2 含鐵菱鎂礦HS態(tài)與LS態(tài)的熱力學(xué)性質(zhì)

圖2(a)給出了菱鎂礦含鐵前后的體積對比。由圖2(a)可以看出，菱鎂礦含鐵后HS態(tài)對應(yīng)的體積比LS態(tài)對應(yīng)的體積大，在30 GPa、300 K的溫壓條件下，(Mg0.5Fe0.5)CO3在HS和LS態(tài)下的體積分別為77.138 10 ?3和73.366 07 ?3，自旋轉(zhuǎn)變后體積減小5.0%左右。對比含鐵菱鎂礦與不含鐵菱鎂礦的體積可知，含鐵菱鎂礦LS態(tài)的體積減小，HS態(tài)的體積則在低溫端增大，在高溫端減小，說明含鐵菱鎂礦的HS態(tài)對應(yīng)的熱膨脹系數(shù)小于不含鐵菱鎂礦。本研究中自旋轉(zhuǎn)變導(dǎo)致的體積變化幅度比Liu等[11]實驗測量得到的變化幅度稍小，這是由于實驗所用樣品為(Mg0.35Fe0.65)CO3，體積變化幅度的差異源于測量樣品中鐵的摩爾分數(shù)不同。

圖2 (Mg0.5Fe0.5)CO3的HS態(tài)與LS態(tài)的體積V（a）和熱膨脹系數(shù)α（b）隨溫度的變化關(guān)系Fig.2 Temperature dependence of the (a) volume V, (b) thermal expansion coefficient α for both the HS and LS state of (Mg0.5Fe0.5)CO3

熱膨脹系數(shù)是研究礦物熱力學(xué)性質(zhì)的重要參數(shù)，本研究計算了含鐵對菱鎂礦熱膨脹系數(shù)的影響。圖2(b)分別給出了含鐵菱鎂礦的HS、LS態(tài)在30 GPa和60 GPa時的熱膨脹系數(shù)，為了便于比較，圖2(b)中給出了對應(yīng)壓力下MgCO3的熱膨脹系數(shù)以及部分前人的結(jié)果[11,32]。由圖2(b)可知，在相同的溫壓條件下，(Mg0.5Fe0.5)CO3的HS與LS兩種自旋態(tài)的熱膨脹系數(shù)均小于MgCO3，即含鐵會導(dǎo)致菱鎂礦的熱膨脹系數(shù)急劇減小。在30 GPa、300 K的溫壓條件下，計算得到MgCO3與(Mg0.5Fe0.5)CO3的熱膨脹系數(shù)分別為 8.05 × 10?6K?1和 3.13 × 10?6K?1，減小幅度為 61.0%；(Mg0.5Fe0.5)CO3的 HS 和 LS 態(tài)的熱膨脹系數(shù)分別為 3.13 × 10?6K?1和 4.22 × 10?6K?1，鐵的自旋轉(zhuǎn)變導(dǎo)致的熱膨脹系數(shù)增加幅度為 34.8%。此外，熱膨脹系數(shù)會隨著壓力的增加而減小，并且對溫度的依賴性減弱。這與含鐵方鎂石在低壓下HS態(tài)的熱膨脹系數(shù)比LS態(tài)大，而在高壓下HS態(tài)的熱膨脹系數(shù)反而比LS態(tài)小的變化趨勢有所差別[28]。根據(jù)格臨愛森定律

由于影響熱膨脹系數(shù)的主要物理量為格臨愛森常數(shù)γ、定容比熱容CV,m、體積模量ΚT及體積V，因此可以從上述熱力學(xué)參數(shù)的變化來解釋含鐵菱鎂礦熱膨脹系數(shù)的變化機制。

圖3給出了(Mg0.5Fe0.5)CO3的熱力學(xué)參數(shù)隨溫度和壓力的變化關(guān)系。由圖3可知，定容比熱容在低溫段先隨溫度升高而急劇增大，隨后趨于飽和，MgCO3的定容比熱容為 231.75 J/(mol·K)，(Mg0.5Fe0.5)CO3在HS、LS態(tài)時的定容比熱容分別為 232.99 J/(mol·K)和 231.58 J/(mol·K)，彼此之間的差異很小，因此定容比熱容對熱膨脹系數(shù)的影響可以忽略不計。由1.2節(jié)的研究結(jié)果可知，菱鎂礦含鐵及鐵的自旋轉(zhuǎn)變引起的體積變化在5.0%左右，無法成為熱膨脹系數(shù)劇烈變化的決定性因素。由圖3(b)和圖3(c)可知，菱鎂礦含鐵后格臨愛森常數(shù)明顯減小，由HS態(tài)轉(zhuǎn)變?yōu)長S態(tài)后明顯增大，與熱膨脹系數(shù)的變化趨勢很好地符合。需要注意的是，含鐵菱鎂礦LS態(tài)的格臨愛森常數(shù)γ比不含鐵時大，而熱膨脹系數(shù)卻比不含鐵時小，由此可以推斷，體變模量KT對熱膨脹系數(shù)的變化也起到至關(guān)重要的作用。因此，菱鎂礦含鐵及鐵的自旋轉(zhuǎn)變導(dǎo)致的熱膨脹系數(shù)變化是由格臨愛森常數(shù)γ與體變模量KT共同決定的。

圖3 MgCO3與（Mg0.5Fe0.5）CO3的定容比熱容CV,m(a)、格臨愛森常數(shù)γ(b)、體變模量KT(c)與溫度的關(guān)系Fig.3 Temperature dependence of the (a) specific heat capacity of constant volume CV,m, (b) Grüneisen parameter γ,(c) bulk modulus KT of MgCO3 and (Mg0.5Fe0.5)CO3

2.3 含鐵菱鎂礦MS態(tài)的熱力學(xué)性質(zhì)

Liu 等[11]、Hsu 等[32]、Fu 等[33]的研究表明，含鐵菱鎂礦的MS態(tài)在高溫高壓下會引起彈性及聲速的異常變化。研究MS態(tài)熱力學(xué)參數(shù)的前提是計算出自旋共存時LS態(tài)的占比n。根據(jù)式（1）～式（8）計算出MS態(tài)下LS的占比n，n對p和T的一階導(dǎo)數(shù)以及α、V、v等熱力學(xué)參數(shù)。

圖4給出了不同溫度下n隨壓力的變化關(guān)系。由圖4可知：在低溫時，n急劇增大到1，即HS態(tài)與LS態(tài)共存的壓力區(qū)間很窄；隨著溫度升高，HS態(tài)與LS態(tài)共存的壓力區(qū)間明顯增大；當溫度為300 K時，HS態(tài)與LS態(tài)共存的壓力區(qū)間約為5 GPa；溫度為1 500 K時，兩者共存的壓力區(qū)間增大到約15 GPa，與Liu等[11]、Hsu等[32]通過實驗和理論計算得到的趨勢一致。同時，通過計算發(fā)現(xiàn)：自旋轉(zhuǎn)變的壓力隨溫度的升高而增加；溫度為300 K時，自旋轉(zhuǎn)變發(fā)生在約40 GPa；而溫度為3 000 K時，自旋轉(zhuǎn)變發(fā)生在約60 GPa。該自旋轉(zhuǎn)變行為在鐵方鎂石的理論計算和實驗研究中也有報道[18,20,22,38-40]。HS態(tài)與LS態(tài)共存區(qū)間和自旋轉(zhuǎn)變隨壓力的變化關(guān)系決定了n對溫度和壓力的一階導(dǎo)數(shù)，如圖5所示。計算表明：當壓力為50 GPa時異常變化的起始溫度為500 K左右；70 GPa時異常變化的起始溫度為1 500 K左右。同時峰的寬度也由50 GPa時的1 000 K增大到70 GPa時的1 700 K，由此可知，自旋轉(zhuǎn)變引起的峰隨著壓強的增加向高溫方向移動，同時峰的寬度也隨之增加。圖5(b)給出了隨壓強的變化關(guān)系，其峰高和峰寬隨壓力的變化趨勢與類似。

圖4 n隨溫度和壓力的變化關(guān)系Fig.4 Ratio n as the function of pressure at different temperature

圖6 MS態(tài)的體積V(a)、熱膨脹系數(shù)α(b)、速度v(c)隨壓力的變化關(guān)系Fig.6 Volume V (a), thermal expansion coefficient α (b), velocity v (c) of MS state as the function of pressure

圖6(b)和圖6(c)分別給出了熱膨脹系數(shù)α和速度v隨壓力的變化關(guān)系。由MS態(tài)時熱膨脹系數(shù)的計算公式(4)可知，α主要受到LS態(tài)的占比n對溫度的一階導(dǎo)數(shù)的影響，其中α與成反比，由圖5(a)可知，隨著HS態(tài)向LS態(tài)發(fā)生自旋轉(zhuǎn)變，在自旋共存區(qū)間內(nèi)的變化趨勢為先突然減小然后增大，因此熱膨脹系數(shù)α的變化趨勢為先突然增大然后減小，這與體積的變化趨勢不同；由式（5）、式（6）、式（8）可知，MS態(tài)的速度主要由LS態(tài)的占比n對壓強的一階導(dǎo)數(shù)決定，且v與成反比，由圖5(b)可知，在自旋共存區(qū)間內(nèi)，的變化趨勢為先突然增大后減小，因此速度v在自旋共存區(qū)間的變化趨勢為先突然減小然后增大。在300 K時，MS態(tài)(Mg0.5Fe0.5)CO3在約40 GPa時產(chǎn)生異常變化，而在1 500 K時，這類異常變化出現(xiàn)在約50 GPa，該變化趨勢與的變化一致，表明熱力學(xué)參數(shù)的變化受n的影響。α、V、v由于自旋轉(zhuǎn)變引起的異常隨溫度的增大向高壓方向移動，并且變得更加平滑，但在2 500 K左右時，α、V、v的異常仍然比較明顯。本研究計算結(jié)果與前人的實驗研究均符合較好，在相同溫度下，僅峰的位置存在較小的偏差，這是由于實驗中Fe2+的濃度較高，其自旋轉(zhuǎn)變的壓力也隨之相應(yīng)地提高。(Mg0.5Fe0.5)CO3的MS態(tài)體積和速度的異常變化也會影響其在地球深部條件下的熱傳導(dǎo)特征。

3 結(jié) 論

采用第一性原理計算方法研究了含鐵菱鎂礦(Mg0.5Fe0.5)CO3中鐵的摻入以及自旋轉(zhuǎn)變對其熱力學(xué)性質(zhì)的影響，得出以下主要結(jié)論。

（1）含鐵后菱鎂礦晶體的體積發(fā)生明顯變化，HS態(tài)在低溫端體積增大，而在高溫端體積比不含鐵時減小，LS態(tài)的體積明顯減小，表明鐵的自旋轉(zhuǎn)變會導(dǎo)致含鐵菱鎂礦晶體體積減小。

（2）含鐵后菱鎂礦的熱膨脹系數(shù)顯著減小，但是鐵的自旋轉(zhuǎn)變導(dǎo)致熱膨脹系數(shù)有所增大，但仍然比不含鐵時要??；格臨愛森常數(shù)與體變模量的變化是熱膨脹系數(shù)變化的主要決定因素。

若要進一步明確含鐵菱鎂礦的熱力學(xué)性質(zhì)對地幔在碳循環(huán)中的作用，還需全面考慮地幔中其他礦物對菱鎂礦的影響，如鐵方鎂石和鈣鈦礦等在地球深部同時存在時所表現(xiàn)出的物理化學(xué)性質(zhì)。