馬麥寧，張吉凱，韓林，張雨心，張璘翼

(中國科學院大學地球與行星科學學院 中國科學院計算地球動力學重點實驗室， 北京 100049)

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

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

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

1 鎂輝石的結構特征

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

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

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

2 彈性波速計算方法

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

2.1 計算參數(shù)選擇

斜方結構的鎂輝石Pbca和P21ca的晶胞都含有80個原子，經(jīng)過能量收斂測試之后，最終選定平面波截斷能為70 Ry，電子密度截斷能為平面波截斷能的4倍，即280 Ry；第一布里淵區(qū)k點取1×2×4[26]。

單斜結構的鎂輝石P21/c和C2/c的晶胞里有40個原子，平面波截斷能設為80 Ry，電子密度截斷能也是平面波截斷能的4倍，為320 Ry；k點取為4×4×4。

截斷能和布里淵區(qū)采樣網(wǎng)格點的選定，能夠保證每一次晶胞結構優(yōu)化之后，原子的能量收斂標準在2.0×10-5Ry/atom，電子能量收斂標準在1.0×10-8Ry，原子受力不超過2.0×10-4Ry/bohr。為保證結果的可靠性，實際計算時先選定Pbca相的幾個壓力點進行測試，并與Kung等[9]給出的實驗數(shù)據(jù)進行對照，驗證計算數(shù)據(jù)的合理性及一致性后，完成4種結構不同壓力下的彈性計算。這樣處理后選定的參數(shù)，能夠保證彈性計算的可靠性[3]。

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

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

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

2.3 彈性波速及其各向異性的計算

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

|cijklnjnl-ρV2δik|=0.

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

因為輝石不僅在上地幔中含量豐富，而且各向異性也很強，因此特別計算了鎂輝石的波速各向異性。

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

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

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

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

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

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

3 計算結果

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

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

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

圖3～圖6分別為4種結構相的軸向橫波速度隨壓力的變化。結果顯示，Pbca和P21/c的3個結晶軸方向快波和慢波的速度差都隨壓力升高逐漸最大；而C2/c和P21ca的橫波速度變化比較復雜：隨壓力升高，有些方向速度差變化不大(如a軸方向)，有些逐漸增大(如P21ca的b軸方向)，有些逐漸減小(如C2/c的b軸方向和P21ca的c軸方向)，而個別是先減小后增大(如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種結構相的縱波速度各向異性如圖7所示，在壓力約低于12 GPa時，C2/c的各向異性最大；在更高的壓力下時，Pbca的各向異性最強，且隨著壓力急劇上升，C2/c相的各向異性越來越弱；壓力高于14 GPa后P21ca相的各向異性最小。

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

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

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

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

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

5 結論與展望

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

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