馬超杰，吳 瀟，馬陽(yáng)陽(yáng)，何開(kāi)華，姬廣富

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

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

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

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

1 研究方法與理論

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

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

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

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

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

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

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

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

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

2 結(jié)果與討論

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

含鐵菱鎂礦(Mg0.5Fe0.5)CO3的LS態(tài)與HS態(tài)的焓差ΔH（HHS-HLS）隨壓力的變化趨勢(shì)如圖1所示。由圖1可知：在低壓時(shí)，ΔH＜ 0，HS態(tài)的焓值小，因此 HS 態(tài)更穩(wěn)定；在高壓時(shí)，ΔH＞ 0，LS 態(tài)的焓值變小，LS態(tài)變得更穩(wěn)定。本研究計(jì)算得到的自旋轉(zhuǎn)變壓力為44.5 GPa，與實(shí)驗(yàn)觀察得到的結(jié)果（40～52 GPa）很好地符合[5,10-11,29-30]，與 Hsu 等[32]的計(jì)算結(jié)果（48 GPa）也符合較好，較小的差異源于鐵在含鐵菱鎂礦中的摩爾分?jǐn)?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)給出了菱鎂礦含鐵前后的體積對(duì)比。由圖2(a)可以看出，菱鎂礦含鐵后HS態(tài)對(duì)應(yīng)的體積比LS態(tài)對(duì)應(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%左右。對(duì)比含鐵菱鎂礦與不含鐵菱鎂礦的體積可知，含鐵菱鎂礦LS態(tài)的體積減小，HS態(tài)的體積則在低溫端增大，在高溫端減小，說(shuō)明含鐵菱鎂礦的HS態(tài)對(duì)應(yīng)的熱膨脹系數(shù)小于不含鐵菱鎂礦。本研究中自旋轉(zhuǎn)變導(dǎo)致的體積變化幅度比Liu等[11]實(shí)驗(yàn)測(cè)量得到的變化幅度稍小，這是由于實(shí)驗(yàn)所用樣品為(Mg0.35Fe0.65)CO3，體積變化幅度的差異源于測(cè)量樣品中鐵的摩爾分?jǐn)?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ù)，本研究計(jì)算了含鐵對(duì)菱鎂礦熱膨脹系數(shù)的影響。圖2(b)分別給出了含鐵菱鎂礦的HS、LS態(tài)在30 GPa和60 GPa時(shí)的熱膨脹系數(shù)，為了便于比較，圖2(b)中給出了對(duì)應(yīng)壓力下MgCO3的熱膨脹系數(shù)以及部分前人的結(jié)果[11,32]。由圖2(b)可知，在相同的溫壓條件下，(Mg0.5Fe0.5)CO3的HS與LS兩種自旋態(tài)的熱膨脹系數(shù)均小于MgCO3，即含鐵會(huì)導(dǎo)致菱鎂礦的熱膨脹系數(shù)急劇減小。在30 GPa、300 K的溫壓條件下，計(jì)算得到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ù)會(huì)隨著壓力的增加而減小，并且對(duì)溫度的依賴(lài)性減弱。這與含鐵方鎂石在低壓下HS態(tài)的熱膨脹系數(shù)比LS態(tài)大，而在高壓下HS態(tài)的熱膨脹系數(shù)反而比LS態(tài)小的變化趨勢(shì)有所差別[28]。根據(jù)格臨愛(ài)森定律

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

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

圖3 MgCO3與（Mg0.5Fe0.5）CO3的定容比熱容CV,m(a)、格臨愛(ài)森常數(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)在高溫高壓下會(huì)引起彈性及聲速的異常變化。研究MS態(tài)熱力學(xué)參數(shù)的前提是計(jì)算出自旋共存時(shí)LS態(tài)的占比n。根據(jù)式（1）～式（8）計(jì)算出MS態(tài)下LS的占比n，n對(duì)p和T的一階導(dǎo)數(shù)以及α、V、v等熱力學(xué)參數(shù)。

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

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

3 結(jié) 論

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

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

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

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