高壓下鎂輝石彈性波速的第一性原理研究*

馬麥寧，張吉?jiǎng)P，韓林，張雨心，張璘翼

(中國(guó)科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院 中國(guó)科學(xué)院計(jì)算地球動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室， 北京 100049)

輝石(pyroxene)主要存在于地殼與上地幔中，是大部分鎂鐵質(zhì)到超鎂鐵質(zhì)火成巖和高級(jí)變質(zhì)巖中的重要組成礦物，同時(shí)也是上地幔中除橄欖石以外含量最多的礦物[1-3]。因此，對(duì)輝石物性的研究有助于更好地理解地球深部的物質(zhì)組成和地震波速異常及動(dòng)力學(xué)過(guò)程等。

目前，國(guó)際上對(duì)地球深部礦物物性的研究主要有兩種方式：一種是高溫高壓實(shí)驗(yàn)，就是利用各種靜態(tài)超高壓實(shí)驗(yàn)設(shè)備，包括大腔體高壓裝置和金剛石壓腔裝置，并與超聲、X射線衍射、紅外和拉曼等多種譜學(xué)測(cè)量方法相結(jié)合，進(jìn)行礦物物性原位測(cè)試[2,4-6]；另一種就是理論計(jì)算，主要是基于第一性原理(first principle)、分子動(dòng)力學(xué)(molecular dynamics)和密度泛函理論(density functional theory)等，對(duì)礦物在溫壓下的物性進(jìn)行模擬研究[3]。實(shí)驗(yàn)不可能完成所有材料在任意溫壓范圍的物性測(cè)試，而理論計(jì)算可以彌補(bǔ)這個(gè)不足。理論計(jì)算能夠提供從微觀到宏觀多尺度時(shí)空范圍內(nèi)的研究手段，是很多地球科學(xué)分支之間的紐帶和橋梁[7]。隨著計(jì)算機(jī)軟硬件的不斷發(fā)展，理論計(jì)算在礦物物性研究中的作用越來(lái)越大，越來(lái)越多的研究者開(kāi)始從事礦物物理計(jì)算研究。

輝石是上地幔中含量第二的礦物，前人的研究顯示輝石的高壓結(jié)構(gòu)非常復(fù)雜，其中斜方輝石在高壓下可能存在Pbca、P21/c、C2/c和P21ca等4種結(jié)構(gòu)相[8-18]。但前人的研究多注重于晶體結(jié)構(gòu)和相變，缺少對(duì)每種結(jié)構(gòu)相的彈性波速的理論研究。高溫高壓實(shí)驗(yàn)對(duì)輝石彈性波速的研究主要有Kung等[9]利用超聲波速測(cè)量并結(jié)合X射線衍射實(shí)驗(yàn)發(fā)現(xiàn)，斜方輝石多晶樣品在9～13/14 GPa出現(xiàn)速度軟化，這可能意味著樣品發(fā)生了Pbca→C2/c的相變；Li等[16]利用大腔體高壓裝置結(jié)合X射線衍射研究斜方鐵輝石的縱橫波速度，發(fā)現(xiàn)與鎂輝石不同；Li和Neuville[19]用超聲波干涉和XRD相結(jié)合的方法測(cè)量透輝石沿1 600 K絕熱地溫線的彈性波速和密度，表明在上地幔深度，透輝石的波速比斜方輝石大1%～3%。而輝石波速各向異性的研究，目前被廣泛引用的數(shù)據(jù)是Anderson[1]報(bào)道的常溫常壓下斜方輝石沿3個(gè)結(jié)晶軸方向的縱波速度。而高壓或高溫下的輝石單晶體不同方向的速度資料非常匱乏。因?yàn)樘烊划a(chǎn)出的斜方輝石以鎂輝石為主，所以本研究選擇鎂輝石為研究對(duì)象，利用第一性原理對(duì)其在高壓下4種結(jié)構(gòu)相的彈性波速進(jìn)行理論計(jì)算。

1 鎂輝石的結(jié)構(gòu)特征

輝石屬于單鏈硅酸鹽礦物，根據(jù)晶體對(duì)稱性的不同，輝石分為斜方輝石(orthopyroxene，Opx)和單斜輝石(clinopyroxene，Cpx)兩個(gè)亞族，前者空間群為Pbca，后者為C2/c或P21/c[20-21]。

輝石主要由共角的TO4四面體和共邊的MO6八面體沿著[100]方向交互成層排列而成。輝石的化學(xué)式一般可寫(xiě)作M2M1T2O6，其中M位一般被Mg、Fe、Al、Ca、Na等金屬陽(yáng)離子占據(jù)，M1通常為6配位八面體，M2為6～8配位的畸變的八面體，當(dāng)M2被較小陽(yáng)離子占據(jù)，輝石就是斜方晶系，而如果M2被較大陽(yáng)離子占據(jù)，M2八面體就會(huì)比M1八面體略大，形成單斜晶系。O位根據(jù)成鍵特征可分為電價(jià)平衡的O1、未充分成鍵的O2，及充分成鍵的O3(又稱為橋氧)。常溫常壓下鎂輝石為斜方晶系，化學(xué)式可以簡(jiǎn)寫(xiě)為Mg2Si2O6，基本結(jié)構(gòu)單元是SiO4單鏈，鏈與鏈之間的八面體空隙由Mg2+占據(jù)，所以鎂輝石的結(jié)構(gòu)為MgO6八面體和SiO4四面體交互成層沿a軸方向排列(圖1)。

SiO4四面體為藍(lán)色，白色為另一層的SiO4四面體，MgO6八面體為黃色。圖1 鎂輝石在[100]方向的投影Fig.1 Projection of Mg-pyroxene onto the [100] direction

2 彈性波速計(jì)算方法

計(jì)算是利用Quantum Espresso軟件包[22]實(shí)現(xiàn)的，運(yùn)用基于密度泛函理論的第一性原理進(jìn)行計(jì)算，使用局域密度近似[23-24]作為交換關(guān)聯(lián)函數(shù)。離子勢(shì)能的近似基于平面波贗勢(shì)理論，其中鎂元素的贗勢(shì)采用Von Barth-Car方法產(chǎn)生的模守恒贗勢(shì)，硅和氧元素的贗勢(shì)采用Vanderbilt方法產(chǎn)生的超軟贗勢(shì)[25]。贗勢(shì)近似的價(jià)電子結(jié)構(gòu)分別是3s23p2(Si)、3s2(Mg)和2s22p4(O)。

2.1 計(jì)算參數(shù)選擇

斜方結(jié)構(gòu)的鎂輝石Pbca和P21ca的晶胞都含有80個(gè)原子，經(jīng)過(guò)能量收斂測(cè)試之后，最終選定平面波截?cái)嗄転?0 Ry，電子密度截?cái)嗄転槠矫娌ń財(cái)嗄艿?倍，即280 Ry；第一布里淵區(qū)k點(diǎn)取1×2×4[26]。

單斜結(jié)構(gòu)的鎂輝石P21/c和C2/c的晶胞里有40個(gè)原子，平面波截?cái)嗄茉O(shè)為80 Ry，電子密度截?cái)嗄芤彩瞧矫娌ń財(cái)嗄艿?倍，為320 Ry；k點(diǎn)取為4×4×4。

截?cái)嗄芎筒祭餃Y區(qū)采樣網(wǎng)格點(diǎn)的選定，能夠保證每一次晶胞結(jié)構(gòu)優(yōu)化之后，原子的能量收斂標(biāo)準(zhǔn)在2.0×10-5Ry/atom，電子能量收斂標(biāo)準(zhǔn)在1.0×10-8Ry，原子受力不超過(guò)2.0×10-4Ry/bohr。為保證結(jié)果的可靠性，實(shí)際計(jì)算時(shí)先選定Pbca相的幾個(gè)壓力點(diǎn)進(jìn)行測(cè)試，并與Kung等[9]給出的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)照，驗(yàn)證計(jì)算數(shù)據(jù)的合理性及一致性后，完成4種結(jié)構(gòu)不同壓力下的彈性計(jì)算。這樣處理后選定的參數(shù)，能夠保證彈性計(jì)算的可靠性[3]。

2.2 彈性系數(shù)的計(jì)算

鎂輝石的初始晶胞結(jié)構(gòu)模型選擇Molin[27]及Thompson和Downs[28]的數(shù)據(jù)。利用Quantumn Espresso軟件包中的PWscf模塊對(duì)每個(gè)結(jié)構(gòu)進(jìn)行給定壓力下的結(jié)構(gòu)優(yōu)化，得到該壓力下的晶胞參數(shù)?？偣矊?duì)4種結(jié)構(gòu)相在0～30 GPa壓力范圍各自9個(gè)壓力點(diǎn)進(jìn)行結(jié)構(gòu)優(yōu)化，并得到相應(yīng)壓力下的晶胞參數(shù)。

彈性系數(shù)cij利用應(yīng)變-應(yīng)力的關(guān)系得到，即通過(guò)對(duì)平衡結(jié)構(gòu)的晶胞施加一個(gè)±1%的應(yīng)變，然后固定晶胞體積，只優(yōu)化原子位置，得到一個(gè)應(yīng)力矩陣之后，通過(guò)胡克定律得到彈性系數(shù)?？紤]輝石結(jié)構(gòu)的空間對(duì)稱性后，斜方結(jié)構(gòu)的鎂輝石含有9個(gè)彈性系數(shù)，分別為c11、c22、c33、c44、c55、c66、c12、c13和c23；而單斜結(jié)構(gòu)的鎂輝石含有13個(gè)彈性系數(shù)(另外4個(gè)分別為c15、c25、c35和c46)。

2.3 彈性波速及其各向異性的計(jì)算

得到不同壓力下的彈性系數(shù)后，利用彈性系數(shù)cijkl和Christoffel方程[29]可以得到沿著不同方向傳播的彈性波速度：

|cijklnjnl-ρV2δik|=0.

式中：n為波的傳播方向，δik為Kroncker函數(shù)，ρ為密度，V為波速。求解上面的方程可得到3個(gè)解，分別代表縱波速度和兩個(gè)橫波速度(快波和慢波)。

因?yàn)檩x石不僅在上地幔中含量豐富，而且各向異性也很強(qiáng)，因此特別計(jì)算了鎂輝石的波速各向異性。

得到縱波速度后，利用下式可以得到縱波的速度各向異性：

Ap=300%(VPmax-VPmin)/(VP(a)+VP(b)+VP(c))，

式中：VP(a)、VP(b)和VP(c)分別代表沿3個(gè)晶體軸方向傳播的縱波速度，VPmax和VPmin分別是3個(gè)值中的最大值和最小值。

而橫波的速度各向異性可通過(guò)下式獲得：

AS=200%(VS1-VS2)/(VS1+VS2)，

式中：VS1和VS2分別代表橫波快波和慢波速度。

3 計(jì)算結(jié)果

對(duì)于鎂輝石4種結(jié)構(gòu)相在高壓下的穩(wěn)定性和相變問(wèn)題，我們已經(jīng)在之前的文章中報(bào)道過(guò)[3]，這里只討論它們的彈性波速計(jì)算結(jié)果。

鎂輝石4種結(jié)構(gòu)相的軸向縱波速度隨壓力的變化如圖2所示：整體上，鎂輝石的軸向縱波速度是隨著壓力的增加而增大，但是當(dāng)壓力大于24 GPa時(shí)，Pbca的b軸方向的縱波速度與壓力呈負(fù)相關(guān)，即呈現(xiàn)出速度弱化趨勢(shì)。4種結(jié)構(gòu)中，除C2/c和個(gè)別低壓點(diǎn)外，其余3種基本上都是VP(a)>VP(c)>VP(b)，只有C2/c為VP(c)>VP(a)>VP(b)。

圖2 不同壓力下鎂輝石4種結(jié)構(gòu)相的軸向縱波速度Fig.2 Axial compressional wave velocities of the four structural phases of Mg-pyroxene as a function of pressure

圖3～圖6分別為4種結(jié)構(gòu)相的軸向橫波速度隨壓力的變化。結(jié)果顯示，Pbca和P21/c的3個(gè)結(jié)晶軸方向快波和慢波的速度差都隨壓力升高逐漸最大；而C2/c和P21ca的橫波速度變化比較復(fù)雜：隨壓力升高，有些方向速度差變化不大(如a軸方向)，有些逐漸增大(如P21ca的b軸方向)，有些逐漸減小(如C2/c的b軸方向和P21ca的c軸方向)，而個(gè)別是先減小后增大(如C2/c的c軸方向)。

圖3 不同壓力下鎂輝石Pbca相的軸向橫波速度Fig.3 Axial shear wave velocity of the Pbca phase of Mg-pyroxene as a function of pressure

圖4 不同壓力下鎂輝石P21/c相的軸向橫波速度Fig.4 Axial shear wave velocity of the P21/c phase of Mg-pyroxene as a function of pressure

圖5 不同壓力下鎂輝石C2/c相的軸向橫波速度Fig.5 Axial shear wave velocity of the C2/c phase of Mg-pyroxene as a function of pressure

圖6 不同壓力下鎂輝石P21ca相的軸向橫波速度Fig.6 Axial shear wave velocity of the P21ca phase of Mg-pyroxene as a function of pressure

4 波速各向異性

鎂輝石4種結(jié)構(gòu)相的縱波速度各向異性如圖7所示，在壓力約低于12 GPa時(shí)，C2/c的各向異性最大；在更高的壓力下時(shí)，Pbca的各向異性最強(qiáng)，且隨著壓力急劇上升，C2/c相的各向異性越來(lái)越弱；壓力高于14 GPa后P21ca相的各向異性最小。

圖7 鎂輝石4種結(jié)構(gòu)相的縱波速度各向異性Fig.7 Compressional wave velocity anisotropy of the four structural phases of Mg-pyroxene

圖8為鎂輝石4種結(jié)構(gòu)相的軸向橫波速度各向異性，可以看出：Pbca相，a軸橫波速度各向異性最明顯；3個(gè)方向的橫波速度各向異性隨著壓力的增大而增大(P>4 GPa)。P21/c相，a軸橫波速度各向異性最明顯；3個(gè)方向的橫波速度各向異性隨著壓力的增大而增大(c軸方向在16～20 GPa之間有一個(gè)減小，20 GPa之后繼續(xù)單調(diào)遞增)。C2/c相，壓力小于24 GPa時(shí)b軸橫波速度各向異性最明顯，更高的壓力下c軸橫波速度各向異性最大；b軸方向的橫波速度各向異性隨著壓力的增大而減?。籧軸先減小，壓力達(dá)到12 GPa后開(kāi)始增大；a軸的變化不明顯。P21ca相，a軸橫波速度各向異性最大；隨著壓力的增大，a軸橫波速度各向異性變化不明顯；b軸橫波速度各向異性隨著壓力的增大而增大，c軸橫波速度各向異性隨著壓力的增大而減小，兩者的交叉點(diǎn)在16 GPa附近，此壓力點(diǎn)之前b軸各向異性最小，之后c軸各向異性最小。

圖8 鎂輝石4種結(jié)構(gòu)相的軸向橫波速度各向異性Fig.8 Shear wave velocity anisotropy of the four structural phases of Mg-pyroxene

以上計(jì)算結(jié)果顯示，不同結(jié)構(gòu)相鎂輝石的波速各向異性特征復(fù)雜且多變。在地幔條件下，根據(jù)Kung等[9]、Nestola等[11]、Jahn[13]、Zhang等[15]、Li等[16]、Finkelstein等[17]、Yu等[30]和Xu等[31]的最新研究結(jié)果，隨著壓力的增加，鎂輝石很有可能相變?yōu)椴煌母邏合嘟Y(jié)構(gòu)。這種相變完全可能引起鎂輝石各向異性特征的改變，而以鎂輝石為代表的斜方輝石是上地幔中的主要組成礦物，因此高壓下鎂輝石的相變引起的波速各向異性的變化非?？赡茉斐缮系蒯Ｖ杏^測(cè)到的地震波速異常。反過(guò)來(lái)，根據(jù)4種結(jié)構(gòu)相的波速計(jì)算結(jié)果，無(wú)論鎂輝石在高壓下相變?yōu)楹畏N結(jié)構(gòu)，其各向異性都有明顯的差別，所以如果上地幔局部地區(qū)的各向異性是由鎂輝石引起的，也可以通過(guò)將地震波觀測(cè)到的各向異性與鎂輝石的4種結(jié)構(gòu)的各向異性特征進(jìn)行對(duì)比分析，就可能判斷出實(shí)際的鎂輝石在上地幔中的真實(shí)相結(jié)構(gòu)。最新的Zhang和Bass通過(guò)對(duì)天然含F(xiàn)e的Pbca相單晶的實(shí)驗(yàn)研究，認(rèn)為Pbca→P21/c的相變不可能是250～325 km深度X不連續(xù)面的成因，反而可能是上地幔更深處地震波各向異性的成因[32]。這體現(xiàn)了對(duì)鎂輝石波速各向異性研究的意義。

鎂輝石在高壓下結(jié)構(gòu)相變和不同結(jié)構(gòu)相各向異性復(fù)雜的原因，本質(zhì)是由其晶體結(jié)構(gòu)導(dǎo)致的。由于鎂輝石的結(jié)構(gòu)為MgO6八面體和SiO4四面體交互成層沿a軸方向排列(圖1)，所以在高壓下a軸方向最抗壓，波速較大；而b軸和c軸的分別與四面體、八面體鏈的滑移以及鏈的旋轉(zhuǎn)、扭結(jié)有關(guān)，所以b軸方向容易變形，波速較小；c軸方向中等[3, 16, 18]。

5 結(jié)論與展望

本文利用第一性原理對(duì)鎂輝石高壓下的4種結(jié)構(gòu)相的彈性波速及其各向異性進(jìn)行研究，發(fā)現(xiàn)：4種結(jié)構(gòu)的鎂輝石彈性波速度均表現(xiàn)出明顯的各向異性；壓力小于12 GPa時(shí)，高壓?jiǎn)涡毕郈2/c的縱波波速各向異性最大，而在更高的壓力下，Pbca的各向異性最大，且隨著壓力急劇上升；除C2/c結(jié)構(gòu)外，其他結(jié)構(gòu)相都是a軸方向的橫波波速各向異性最明顯。研究結(jié)果為上地幔的地震波速異常的解釋提供了新的依據(jù)。

目前的計(jì)算只考慮壓力，并未考慮溫度的影響。而實(shí)際地幔中是高溫高壓環(huán)境，溫度的加入將會(huì)使輝石的結(jié)構(gòu)、彈性波速及其各向異性變得更為復(fù)雜。但是，僅僅壓力的改變，就會(huì)引起輝石結(jié)構(gòu)和彈性波速及其各向異性的明顯變化，所以，未來(lái)對(duì)于輝石礦物物理的研究，不僅要繼續(xù)關(guān)注其相變問(wèn)題，更要重視其彈性各向異性及其相變前后的變化。