铁的自旋转变对菱镁矿热力学性质的影响

马超杰，吴 潇，马阳阳，何开华，姬广富

（1.中国地质大学（武汉）数学与物理学院，湖北 武汉 430074；2.中国工程物理研究院流体物理研究所冲击波物理与爆轰物理重点实验室，四川 绵阳 621999）

含碳固溶体的存在会影响地球内部的物理和化学性质[1-4]。近年来，菱镁矿（MgCO3）被认为是碳进入地球深部的主要载体之一，因其在地球深部碳循环中的关键作用而引起广泛关注[5-13]。Isshiki等[14]、Oganov等[15]的研究表明，菱镁矿在下地幔的温压条件下能够稳定存在；Hazen等[13]、Oganov等[15]开展的高温高压实验研究揭示，菱镁矿和菱铁矿的固溶体[(Mg,Fe)CO3]可以在低于100 GPa的压力条件下稳定存在。

铁是多价态的过渡金属，会对菱镁矿和菱铁矿固溶体的性质产生非常重要的影响。此外，铁的自旋在一定的压力和温度条件下可以发生转变，引起物理性质变化[16]。国内外学者对铁方镁石(Mg,Fe)O的自旋转变开展了大量研究，结果表明，铁方镁石中的铁（Fe2+）在40～50 GPa范围内会从高自旋（High spin, HS）态向低自旋（Low spin, LS）态转变，并伴随着结构、电子、光学、弹性和热力学性质的异常变化[16-28]。对于含铁的菱镁矿，已有的高压穆斯堡尔光谱[29]、X射线发射光谱[30]、激光拉曼光谱[10]、X射线衍射[5,8,11]等实验研究和第一性原理计算研究[31-32]都证实，(Mg,Fe)CO3中的铁在约45 GPa时从HS态向LS态转变，导致(Mg,Fe)CO3的体积比MgCO3减小6%～10%[11]。Liu等[11]、Hsu等[32]、Fu等[33]利用理论计算和实验相结合的方法，详细讨论了弹性和地震波速在自旋转变时的变化，得到了非常有意义的结果。其他热力学参数如热膨胀系数、格临爱森常数和比热容等在自旋转变条件下的性质尚缺少报道。

本研究利用第一性原理计算方法，开展高温高压下(Mg,Fe)CO3在HS、LS及HS和LS态共存时的混合自旋（Mixed spin, MS）态的热力学性质（热膨胀系数、体变模量、体积、速度等热力学参数）研究，并与不含铁的MgCO3的相关性质进行对比，分析引起(Mg,Fe)CO3热力学性质变化的机制。研究结果可为研究地幔深部碳的行为以及地幔在全球碳循环中的作用提供制约因素。

1 研究方法与理论

1.1 含铁菱镁矿的结构和计算参数

菱镁矿MgCO3属于三角晶系，其空间群为，原胞中包含10个原子（2个Mg原子、2个C原子和6个O原子）。为了讨论铁及其在高温高压下的自旋转变对菱镁矿物理性质的影响，本研究选择用1个Fe原子替换1个Mg原子，得到Fe和Mg的摩尔比为1∶1，分子式为（Mg0.5Fe0.5）CO3。

几何结构优化和相关能量的计算[34]采用基于密度泛函理论（Density of functional theory, DFT）的第一性原理分子动力学计算软件VASP完成。由于结构中有铁存在，需要考虑强关联作用，因此计算中考虑了Hubbard参数U的影响，即LDA+U。已有的研究表明，U值会随自旋转变而改变，本研究选取Tsuchiya等[18]、Krukau等[35]通过线性响应理论计算得到的U值，分别为ULS= 5.3 eV，UHS= 4.0 eV。K点采用Monkhorst-Pack方法生成以Γ点为中心的15 × 15 × 15网格，截断能设置为1 000 eV[36]。总能量的收敛阈值设置为1 × 10−6eV/cell，原子力的收敛阈值为1 × 10−3eV/cell。采用基于准谐近似的PHONOPY软件计算得到二阶原子间力常数（IFCs）、声子谱和其他热力学参数[37]。构造求解高对称q点声子频率本征值的动力学矩阵，然后利用倒空间中动力学矩阵的傅里叶变换，计算其他一般q点的声子频率。所有计算中，以Γ为中心的q点网格选取为30 × 30 × 30。计算热力学参数之前，用VASP软件完成不同位移或不同体积下结构的能量计算，计算过程中的参数设置与上述几何优化参数设置相同。

1.2 MS态热力学参数

根据已有的研究结果，MS态时的吉布斯自由能可表示为

式中：n为LS态在MS态中所占的百分比，p为压力，T为温度，GHS、GLS分别为HS态和LS态的吉布斯自由能，Gmix为HS-LS混合态的吉布斯自由能。Gmix可以表示为

式中：XFe为铁在含铁菱镁矿中的摩尔分数，kB为玻尔兹曼常数。在给定的压力和温度条件下，通过求最小化MS态吉布斯自由能得到LS态的占比n

式中：ΔGLS−HS为(Mg0.5Fe0.5)CO3在LS态和HS态的吉布斯自由能之差。对于Fe2+来说，自旋和轨道简并量子数分别为s= 2，m= 3。得到n后可以推导出MS态的热膨胀系数α(n)和等温体积模量KT(n)分别为

式中：KS为绝热体积模量，CV,m为定容比热容，ρ为密度。根据上述公式可以计算出MS态下的体积V、速度v和格临爱森常数γ。

2 结果与讨论

2.1 Fe2+的自旋转变

含铁菱镁矿(Mg0.5Fe0.5)CO3的LS态与HS态的焓差ΔH（HHS-HLS）随压力的变化趋势如图1所示。由图1可知：在低压时，ΔH＜ 0，HS态的焓值小，因此 HS 态更稳定；在高压时，ΔH＞ 0，LS 态的焓值变小，LS态变得更稳定。本研究计算得到的自旋转变压力为44.5 GPa，与实验观察得到的结果（40～52 GPa）很好地符合[5,10-11,29-30]，与 Hsu 等[32]的计算结果（48 GPa）也符合较好，较小的差异源于铁在含铁菱镁矿中的摩尔分数不同。

图1 (Mg0.5Fe0.5)CO3的LS态与HS态的焓差Δ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态与LS态的热力学性质

图2(a)给出了菱镁矿含铁前后的体积对比。由图2(a)可以看出，菱镁矿含铁后HS态对应的体积比LS态对应的体积大，在30 GPa、300 K的温压条件下，(Mg0.5Fe0.5)CO3在HS和LS态下的体积分别为77.138 10 Å3和73.366 07 Å3，自旋转变后体积减小5.0%左右。对比含铁菱镁矿与不含铁菱镁矿的体积可知，含铁菱镁矿LS态的体积减小，HS态的体积则在低温端增大，在高温端减小，说明含铁菱镁矿的HS态对应的热膨胀系数小于不含铁菱镁矿。本研究中自旋转变导致的体积变化幅度比Liu等[11]实验测量得到的变化幅度稍小，这是由于实验所用样品为(Mg0.35Fe0.65)CO3，体积变化幅度的差异源于测量样品中铁的摩尔分数不同。

图2 (Mg0.5Fe0.5)CO3的HS态与LS态的体积V（a）和热膨胀系数α（b）随温度的变化关系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

热膨胀系数是研究矿物热力学性质的重要参数，本研究计算了含铁对菱镁矿热膨胀系数的影响。图2(b)分别给出了含铁菱镁矿的HS、LS态在30 GPa和60 GPa时的热膨胀系数，为了便于比较，图2(b)中给出了对应压力下MgCO3的热膨胀系数以及部分前人的结果[11,32]。由图2(b)可知，在相同的温压条件下，(Mg0.5Fe0.5)CO3的HS与LS两种自旋态的热膨胀系数均小于MgCO3，即含铁会导致菱镁矿的热膨胀系数急剧减小。在30 GPa、300 K的温压条件下，计算得到MgCO3与(Mg0.5Fe0.5)CO3的热膨胀系数分别为 8.05 × 10−6K−1和 3.13 × 10−6K−1，减小幅度为 61.0%；(Mg0.5Fe0.5)CO3的 HS 和 LS 态的热膨胀系数分别为 3.13 × 10−6K−1和 4.22 × 10−6K−1，铁的自旋转变导致的热膨胀系数增加幅度为 34.8%。此外，热膨胀系数会随着压力的增加而减小，并且对温度的依赖性减弱。这与含铁方镁石在低压下HS态的热膨胀系数比LS态大，而在高压下HS态的热膨胀系数反而比LS态小的变化趋势有所差别[28]。根据格临爱森定律

由于影响热膨胀系数的主要物理量为格临爱森常数γ、定容比热容CV,m、体积模量ΚT及体积V，因此可以从上述热力学参数的变化来解释含铁菱镁矿热膨胀系数的变化机制。

图3给出了(Mg0.5Fe0.5)CO3的热力学参数随温度和压力的变化关系。由图3可知，定容比热容在低温段先随温度升高而急剧增大，随后趋于饱和，MgCO3的定容比热容为 231.75 J/(mol·K)，(Mg0.5Fe0.5)CO3在HS、LS态时的定容比热容分别为 232.99 J/(mol·K)和 231.58 J/(mol·K)，彼此之间的差异很小，因此定容比热容对热膨胀系数的影响可以忽略不计。由1.2节的研究结果可知，菱镁矿含铁及铁的自旋转变引起的体积变化在5.0%左右，无法成为热膨胀系数剧烈变化的决定性因素。由图3(b)和图3(c)可知，菱镁矿含铁后格临爱森常数明显减小，由HS态转变为LS态后明显增大，与热膨胀系数的变化趋势很好地符合。需要注意的是，含铁菱镁矿LS态的格临爱森常数γ比不含铁时大，而热膨胀系数却比不含铁时小，由此可以推断，体变模量KT对热膨胀系数的变化也起到至关重要的作用。因此，菱镁矿含铁及铁的自旋转变导致的热膨胀系数变化是由格临爱森常数γ与体变模量KT共同决定的。

图3 MgCO3与（Mg0.5Fe0.5）CO3的定容比热容CV,m(a)、格临爱森常数γ(b)、体变模量KT(c)与温度的关系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态的热力学性质

Liu 等[11]、Hsu 等[32]、Fu 等[33]的研究表明，含铁菱镁矿的MS态在高温高压下会引起弹性及声速的异常变化。研究MS态热力学参数的前提是计算出自旋共存时LS态的占比n。根据式（1）～式（8）计算出MS态下LS的占比n，n对p和T的一阶导数以及α、V、v等热力学参数。

图4给出了不同温度下n随压力的变化关系。由图4可知：在低温时，n急剧增大到1，即HS态与LS态共存的压力区间很窄；随着温度升高，HS态与LS态共存的压力区间明显增大；当温度为300 K时，HS态与LS态共存的压力区间约为5 GPa；温度为1 500 K时，两者共存的压力区间增大到约15 GPa，与Liu等[11]、Hsu等[32]通过实验和理论计算得到的趋势一致。同时，通过计算发现：自旋转变的压力随温度的升高而增加；温度为300 K时，自旋转变发生在约40 GPa；而温度为3 000 K时，自旋转变发生在约60 GPa。该自旋转变行为在铁方镁石的理论计算和实验研究中也有报道[18,20,22,38-40]。HS态与LS态共存区间和自旋转变随压力的变化关系决定了n对温度和压力的一阶导数，如图5所示。计算表明：当压力为50 GPa时异常变化的起始温度为500 K左右；70 GPa时异常变化的起始温度为1 500 K左右。同时峰的宽度也由50 GPa时的1 000 K增大到70 GPa时的1 700 K，由此可知，自旋转变引起的峰随着压强的增加向高温方向移动，同时峰的宽度也随之增加。图5(b)给出了随压强的变化关系，其峰高和峰宽随压力的变化趋势与类似。

图4 n随温度和压力的变化关系Fig.4 Ratio n as the function of pressure at different temperature

图6 MS态的体积V(a)、热膨胀系数α(b)、速度v(c)随压力的变化关系Fig.6 Volume V (a), thermal expansion coefficient α (b), velocity v (c) of MS state as the function of pressure

图6(b)和图6(c)分别给出了热膨胀系数α和速度v随压力的变化关系。由MS态时热膨胀系数的计算公式(4)可知，α主要受到LS态的占比n对温度的一阶导数的影响，其中α与成反比，由图5(a)可知，随着HS态向LS态发生自旋转变，在自旋共存区间内的变化趋势为先突然减小然后增大，因此热膨胀系数α的变化趋势为先突然增大然后减小，这与体积的变化趋势不同；由式（5）、式（6）、式（8）可知，MS态的速度主要由LS态的占比n对压强的一阶导数决定，且v与成反比，由图5(b)可知，在自旋共存区间内，的变化趋势为先突然增大后减小，因此速度v在自旋共存区间的变化趋势为先突然减小然后增大。在300 K时，MS态(Mg0.5Fe0.5)CO3在约40 GPa时产生异常变化，而在1 500 K时，这类异常变化出现在约50 GPa，该变化趋势与的变化一致，表明热力学参数的变化受n的影响。α、V、v由于自旋转变引起的异常随温度的增大向高压方向移动，并且变得更加平滑，但在2 500 K左右时，α、V、v的异常仍然比较明显。本研究计算结果与前人的实验研究均符合较好，在相同温度下，仅峰的位置存在较小的偏差，这是由于实验中Fe2+的浓度较高，其自旋转变的压力也随之相应地提高。(Mg0.5Fe0.5)CO3的MS态体积和速度的异常变化也会影响其在地球深部条件下的热传导特征。

3 结 论

采用第一性原理计算方法研究了含铁菱镁矿(Mg0.5Fe0.5)CO3中铁的掺入以及自旋转变对其热力学性质的影响，得出以下主要结论。

（1）含铁后菱镁矿晶体的体积发生明显变化，HS态在低温端体积增大，而在高温端体积比不含铁时减小，LS态的体积明显减小，表明铁的自旋转变会导致含铁菱镁矿晶体体积减小。

（2）含铁后菱镁矿的热膨胀系数显著减小，但是铁的自旋转变导致热膨胀系数有所增大，但仍然比不含铁时要小；格临爱森常数与体变模量的变化是热膨胀系数变化的主要决定因素。

若要进一步明确含铁菱镁矿的热力学性质对地幔在碳循环中的作用，还需全面考虑地幔中其他矿物对菱镁矿的影响，如铁方镁石和钙钛矿等在地球深部同时存在时所表现出的物理化学性质。