第一性原理研究钼酸铅中氧填隙光谱

First Principles Study of Interstitial Oxygen Spectra in Lead Molybdate Crystals

摘要: 由氧填隙引起的光致变色效应出现在淡黄色的钼酸铅晶体中，这严重影响了晶体在光学仪器中的应用。因此，了解钼酸铅晶体中由内在缺陷引起的特性是非常重要的。在这项工作中，通过考虑有限尺寸(FNV)校正，获得了缺陷形成能量。在考虑电子–声子耦合的情况下，计算得到了氧填隙的光谱。计算结果表明，钼酸铅中氧填隙的光谱峰值分别为1.73 eV和2.15 eV，与实验值2.14 eV (580 nm)非常吻合。对应于 O_I^′ 和 O_I^″的发射峰的位置为1.23 eV和1.97 eV与实验结果比较接近。吸收带580纳米与我们计算得到的氧填隙的吸收带相吻合，所以我们预测580纳米的吸收带来自氧填隙缺陷。

Abstract: A photochromic effect caused by interstitial oxygen presents in the yellowish lead molybdate crystals, which seriously affects the application of the crystals in optical instruments. So it is essential to understand the properties caused by intrinsic defects in lead molybdate crystals. In this work, the defect formation energy has been obtained by taking the finite-size (FNV) correction into consideration. The optical spectra also have been gotten under the consideration of electron-phonon coupling. The calculated results indicate that the spectra of the interstitial oxygen in lead molybdate peaks are of 1.73 eV and 2.15 eV, respectively, which are in very good agreement with the experimental value of 2.14 eV (580 nm). The positions of the emission peak corresponding to the O_I^′ and O_I^″ defects are 1.23 eV and 1.97 eV, which are consistent with the experimental results (1.8 eV and 2.1 eV). The absorption band 580 nm coincides with the absorption band obtained from our calculation of the interstitial oxygen, so we predict that the 580 nm absorption band comes from the interstitial oxygen.

1. 引言

钼酸铅(PbMoO₄)是一种广泛应用的材料，在光催化、太阳能电池、污染物降解和农业等不同领域都有应用，同位素100Mo的存在也使得钼酸铅可以被用作双β衰变的检测器 [1] - [7] 。然而，这些应用都取决于钼酸铅的纯度。钼酸铅中的缺陷对材料的特性有很大影响。例如，由钼酸铅中的色心引起的光致变色效应对其在光学仪器中的应用是有限制的。值得一提的是，钼酸铅的晶体特性也受到晶体生长方法的强烈影响。掺杂剂、主要成分的偏差和紫外线辐射都会对钼酸铅的缺陷产生影响。因此，对钼酸铅中缺陷的研究对于控制晶体生长条件和避免其中的负面影响是必要的 [5] [8] - [14] 。

MoO₃和PbO是PbMoO₄生长的原始材料。MoO₃和PbO具有不同的气压，PbO具有较高的挥发性 [15] 。因此，如果钼酸铅生在富氧的环境中，会形成大量的氧填隙。经过过去十年对钼酸铅材料的广泛研究，材料中空位缺陷对材料光学性质的影响已相对成熟。氧填隙缺陷是钼酸铅中典型的内在点缺陷，直接或间接地控制着材料的许多特性。然而，钼酸铅中氧填隙缺陷的研究结果仍然是目前研究中的一个空白，特别是对带电缺陷的缺陷形成能计算和缺陷的吸收谱与发射谱的计算 [16] [17] [18] [19] 。

M. Buryi等人利用EPR和TSL实验研究了钼酸铅单晶中的电荷捕获和能量转移过程 [20] 。实验结果对钼酸铅晶体中的电子跃迁有着很好的指导意义。在这项工作中，我们使用密度泛函理论(DFT) [7] [21] 计算了氧填隙的电子结构。我们采用Freysoldt-Neugebauer-Van de Walle (FNV) [22] 校正来抵消缺陷状态的周期性图像，以获得准确的缺陷形成能量。带隙是通过杂化泛函方法调整的。通过考虑电声耦合效应，研究了氧填隙的光学性质。

2. 计算方法

2.1. 计算模型

钼酸铅晶体(PMO)的原胞由12个原子组成，空间群结构为I4₁/a。在本次计算中我们使用的计算模型为钼酸铅原胞2 × 2 × 2扩胞之后得到的96个原子的晶胞 [23] 。我们通过在合适的位置加入一个氧原子在计算氧填隙的性质(如图1所示)，氧填隙缺陷可能失去0个电子(O_I)，1个电子( ${\text{O}}_{\text{I}}^{-}$ )或者2个电子( ${\text{O}}_{\text{I}}^{2-}$ )。

Figure 1. Lead molybdate with interstitial oxygen

图1. 含有氧填隙缺陷的钼酸铅晶胞

2.2. 计算参数设置

所有缺陷结构的结构优化使用Vienna Ab-initio Simulation Package (VASP) [24] 完成。电子之间的相互作用选择projector argument wave (PAW) [25] 方法进行描述，交换–相关势由generalized gradient approximation of Perdew-Becke-Ernzerh (GGA-PBE) [26] 方法描述。选用的价电子为：Mo：4d⁵5s¹，Pb：6s²6p²，O：2s²2p⁴。在计算中选择了800 eV的截断能和3 × 3 × 3的K点设置。

2.3. 带边修正

使用DFT计算的材料带隙值往往比真实值偏小，考虑到后续计算中电子跃迁的问题需要，必须要对此问题进行修正。我们使用杂化泛函(HSE06) [27] 来修正带隙问题，此时交换–相关函数为：

$E_{xc}^{HSE}=\alpha E_x^{Fock}-\left(1-\alpha \right)E_x^{PBE}$ (1)

该方程描述了新的电子的交换–相关作用，可以通过调节参数 $\alpha$ 来获得更准确的带隙值。参数 $\alpha$ 对应Fock交换相关势和PBE交换相关势的权重。经过反复测试，我们在此次计算中设定 $\alpha =0.15$ ，计算得到的钼酸铅带隙值与实验值较为吻合。计算的对比结果我们列在表1中。

Table 1. Comparison of PMO bandgap values calculated using different functional functions with experimental values (in eV)

表1. 使用不同泛函计算的PMO带隙值与实验值的比较(单位：eV)

计算得到HSE与PBE的价带顶(VBM)与导带底(CBM)位置之后，我们进行带边修正，缺陷转化能级位置如图2所示。其中common reference level为公共参考能级。

3. 结果和讨论

3.1. 结构优化

在对含有氧填隙的晶体进行结构优化之后，氧填隙结构发生了明显的结构弛豫。在引入一个氧原子之后，填隙周围最近邻的原子离填隙的距离会发生变化，这是因为填隙周围的氧原子受到库仑力和机械力的作用。与完美晶体的结构相比，氧填隙一般呈现二价负电性( ${\text{O}}_{\text{I}}^{2-}$ )，因此一旦产生填隙，由填隙引起的两类势将出现，其一是氧填隙引起的电动势，这将吸引填隙周围Pb²⁺和Mo⁶⁺，而将带负电的O²⁻推开：其二是由于填隙本身产生的机械势，填隙结构带带电后，将于周围的Mo⁶⁺成键，此时填隙与Mo⁶⁺之间的相互作用大于产生负电产生的电势，使得周围的Mo⁶⁺和O²⁻离子更加接近填隙位置，且此时填隙本身的机械势大于电势，推开了Pb²⁺。晶格弛豫情况在表2中给出。

Figure 2. Schematic diagram of band edge correction for interstitial oxygen

图2. 氧填隙的带边修正示意图

Table 2. The distance variation between the vacancy and the nearest oxygen atom before and after structural relaxation (in Å)

表2. 结构优化前后填隙周围最近邻原子到氧填隙的位置(单位：Å)

3.2. 缺陷形成能修正

带电缺陷形成能的计算精度一直是业界难题，当一个晶胞中存在一个带电点缺陷之后，由于晶体的周期性排列，则带电缺陷也会周期性排列，带电缺陷之间就会产生不能忽略的相互作用，因此必须对缺陷形成能进行修正。

带电点缺陷的缺陷形成能可以表示为 [28] ：

$E^f\left(D^q\right)=E\left(D^q\right)-E\left(\text{perfect}\right)+{\mu }_O+q\left(E_{VBM}+E_f-\Delta V\right)-E_{corr}$ (2)

其中 $E\left(D^q\right)$ 表示含有缺陷晶体的总能量， $E\left(\text{perfect}\right)$ 是完整晶体的总能量， ${\mu }_O$ 是氧的化学势， $E_{VBM}$ 指的是完美晶体的价带顶， $E_f$ 是费米能级位置， $\Delta V$ 为完整晶体和含有缺陷晶体的电势差。 $E_{corr}$ 即缺陷形成能的修正项，可以表示为：

$E_{corr}=E_{lat}^q-q{\Delta }_{\frac{q}{0}}$ (3)

其中 $E_{lat}^q$ 是晶格屏蔽能， $q{\Delta }_{\frac{q}{0}}$ 表示的短程势：

$E_q^{lat}={\displaystyle \int {\text{d}}^3\stackrel{\to }{r}}\left[q_d\left(\stackrel{\to }{r}\right)+n\right]\left[{\tilde{V}}_q^{lr}\left(\stackrel{\to }{r}\right)-V_q^{lr}\left(\stackrel{\to }{r}\right)+n{V}_q^{lr}\left(\stackrel{\to }{r}\right)\right]$ (4)

其中 $q_d\left(\stackrel{\to }{r}\right)$ 为缺陷位置的电荷分布，n为背景电荷，加入背景电荷之后即可消除由于带电缺陷引起的静电相互作用， ${\tilde{V}}_q^{lr}\left(\stackrel{\to }{r}\right)$ 为周期性结构中的长程电势分布， $V_q^{lr}\left(\stackrel{\to }{r}\right)$ 为孤立点缺陷的长程电势分布，下面我们对方程(4)进行详细解释：

在缺陷附近，模型电荷与实际电势分布差异较大，在远离缺陷位置，两者电势分布非常接近，所以对于缺陷附近，我们选择加入短程电势。因此，首先应该确定缺陷位置的电荷分布：

$q_d\left(\stackrel{\to }{r}\right)=q{\left|{\Psi }_d\left(\stackrel{\to }{r}\right)\right|}^2$ (5)

进一步考虑带电缺陷周围电荷的电势，我们引入：

$V_{\frac{q}{0}}\left(\stackrel{\to }{r}\right)=V_{defect,q}\left(\stackrel{\to }{r}\right)-V_{defect,0}\left(\stackrel{\to }{r}\right)$ (6)

其中 $V_{defect,q}\left(\stackrel{\to }{r}\right)$ 带电量为q的缺陷电势分布， $V_{defect,0}\left(\stackrel{\to }{r}\right)$ 为带电量为0的缺陷电势分布，二者相减本质即为电荷q产生的电势，又考虑到周期性结构与补偿电荷的存在，通过对实空间中的电势分布进行傅里叶变换，可以得到倒格子空间中的电势分布：

$V_{\frac{q}{0}}^{rec}\left(\stackrel{\to }{G}\right)={\displaystyle \int {\text{d}}^3{\stackrel{\to }{r}}^{\prime }}{V}_{\frac{q}{0}}\left({\stackrel{\to }{r}}^{\prime }\right){\text{e}}^{-i\stackrel{\to }{G}·{\stackrel{\to }{r}}^{\prime }}$ (7)

${\tilde{V}}_{\frac{q}{0}}\left(\stackrel{\to }{r}\right)=\frac{1}{\text{Ω}}{\displaystyle \underset{G\ne 0}{\sum }{V}_{\frac{q}{0}}^{rec}\left(\stackrel{\to }{G}\right){\text{e}}^{i\stackrel{\to }{G}·\stackrel{\to }{r}}}$ (8)

式(7)为考虑到周期性结构的电势分布。

计算孤立点缺陷形成能，周期性的带电缺陷之间有不真实的库仑相互作用，使用下式对该相互作用进行消除：

$E^{inter}=\frac{1}{2}{\displaystyle \int {\text{d}}^3\stackrel{\to }{r}}\left[q_d\left(\stackrel{\to }{r}\right)+n\right]\left[{\tilde{V}}_{\frac{q}{0}}\left(\stackrel{\to }{r}\right)-V_{\frac{q}{0}}\left(\stackrel{\to }{r}\right)\right]$ (9)

为了使计算收敛，计算对带电体系设置了了等量异号均匀分布的背景电荷。背景电荷实际体系并不存在，需要剔除它的影响：

$E^{intra}={\displaystyle \int {\text{d}}^3\stackrel{\to }{r}}n{V}_{\frac{q}{0}}\left(\stackrel{\to }{r}\right)=-q\left(\frac{1}{\text{Ω}}{\displaystyle \int {\text{d}}^3\stackrel{\to }{r}}{V}_{\frac{q}{0}}\left(\stackrel{\to }{r}\right)\right)$ (10)

其中

$V_{\frac{q}{0}}\left(\stackrel{\to }{r}\right)=V_q^{lr}\left(\stackrel{\to }{r}\right)+V_{\frac{q}{0}}^{sr}\left(\stackrel{\to }{r}\right)$ (11)

$V_q^{lr}\left(\stackrel{\to }{r}\right)=\frac{1}{\epsilon }{\displaystyle \int {\text{d}}^3{\stackrel{\to }{r}}^{\prime }}\frac{q_d\left({\stackrel{\to }{r}}^{\prime }\right)}{\left|\stackrel{\to }{r}-{\stackrel{\to }{r}}^{\prime }\right|}$ (12)

考虑了周期性结构之后，短程势应该为超原胞中所有点缺陷附近电势的叠加

${\tilde{V}}_{\frac{q}{0}}^{sr}\left(\stackrel{\to }{r}\right)={\displaystyle \underset{R}{\sum }{V}_{\frac{q}{0}}^{sr}\left(\stackrel{\to }{r}+\stackrel{\to }{R}\right)}+C$ (13)

其中 $\stackrel{\to }{R}$ 为倒格子基矢，则其他原胞的短程势即可忽略不计，因此，式(4.12) $\approx V_{\frac{q}{0}}^{sr}\left(\stackrel{\to }{r}\right)+C,\stackrel{\to }{r}\in \text{Ω}$

$\begin{array}{c}{V}^{sr}\left(\stackrel{\to }{r}\right)={\tilde{V}}^{sr}\left(\stackrel{\to }{r}\right)-C={\tilde{V}}_{\frac{q}{0}}\left(\stackrel{\to }{r}\right)-V_q^{lr}\left(\stackrel{\to }{r}\right)-C\\ ={\tilde{V}}^{els}\left(\text{charged}\right)-{\tilde{V}}^{els}\left(\text{neutral}\right)-{\stackrel{\to }{V}}_q^{lr}-\text{Δ}V\end{array}$ (14)

其中， ${\tilde{V}}_{\frac{q}{0}}\left(\stackrel{\to }{r}\right)$ ，C以及介电常数 $\epsilon$ 都可以通过DFT计算得到， $E_{corr}=E^{inter}+E^{intra}$ ，至此，可以计算缺陷形成能的修正项。

缺陷电荷模型：

$q_d\left(\stackrel{\to }{r}\right)=qx{N}_{\gamma }^{-1}{\text{e}}^{-\frac{\left|\stackrel{\to }{r}-{\stackrel{\to }{r}}_0\right|}{\gamma }}+q\left(1-x\right)N_{\beta }^{-1}{\text{e}}^{-\frac{{\left|\stackrel{\to }{r}-{\stackrel{\to }{r}}_0\right|}^2}{{\beta }^2}}$ (15)

其中 ${\stackrel{\to }{r}}_0$ 为缺陷的中心位置，等式右边第一项是指数函数，第二项为高斯函数，高斯衰减速度大于指数衰减速度，因此，当电荷比较局域，

我们使用高斯衰减，x为选择函数的权重。 $N_{\gamma }=8\pi {\gamma }^3$ ， $N_{\beta }={\pi }^{\frac{3}{2}}{\beta }^3$ ，是归一化常数。图3为氧填隙缺陷形成能修正短程势与长程势的拟合，图中虚线即是短程势，横坐标中间位置表示远离缺陷的位置，远离缺陷位置趋近于0，表示修正达到了良好的效果。

Figure 3. Correction of defect formation energy

图3. 缺陷形成能修正

我们在表3中给出修正前后的缺陷形成能。

Table 3. Interstitial oxygen defect formation energy with and without correction (in eV)

表3. 氧填隙缺陷修正前后的缺陷形成能(单位：eV)

3.3. 氧填隙的光学性质

在获得缺陷形成能和缺陷转化能级的位置之后，通过考虑电子–声子耦合效应，我们可以近似的计算出氧填隙缺陷对应的吸收光谱和发射光谱。

根据弗兰克–康登原理(Frank-Condon Principe) [29] ，我们使用公式(4)来计算发射谱 [30]

$L\left(\epsilon \right)={\displaystyle \underset{n}{\sum }{\text{e}}^{-s}\frac{S^n}{n!}{g}_{\sigma }\left(E_{ZPL}-n{\epsilon }_0-\epsilon \right)}$ (16)

其中S表示黄昆因子， $S_g=\frac{\Delta E_g}{{\Omega }_g}$ ( ${\Omega }_g$ 为有效频率)， $g_{\sigma }$ 为高斯函数， $E_{ZPL}$ 为零声子线。 ${\epsilon }_0$ 平均声子能量。另外，用于计算有效质量的参数 ${\left(\text{Δ}Q\right)}^2={\displaystyle {\sum }_{\alpha ,i}{m}_{\alpha }\text{Δ}{R}_{\alpha ,i}^2}$ ， $\text{Δ}Q=M^{\frac{1}{2}}\text{Δ}R$ ， $\alpha$ 代表原子， $i=\left\{x,y,z\right\}$ ， $\text{Δ}R$ 为位移矢量。计算光谱的有效参数列于表4。

Table 4. Calculation of the effective parameters of the interstitial oxygen defect spectra

表4. 计算氧填隙缺陷光谱的有效参数

对于吸收谱的计算：

$f\left(h\nu \right)={\displaystyle \underset{p-\infty }{\overset{\infty }{\sum }}\delta }\left(E_0-p\hslash \omega -h\nu \right)W_p$ (17)

$W_p=\exp \left[-\left(2n+1\right)S\right]{\left[\frac{\bar{n}+1}{\bar{n}}\right]}^{\frac{p}{2}}\times I_p\left(2S{\left[\left(\bar{n}+1\right)\bar{n}\right]}^{\frac{1}{2}}\right)$ (18)

$\bar{n}=\frac{1}{\exp \left(\frac{\hslash \omega }{kT}\right)-1}$ (19)

其中 $E_0$ 是转化能量， $\hslash \omega$ 是声子能量； $h\nu$ 是光子能量； $I_p$ 是贝塞尔函数。

计算光谱的本质即考虑电子跃迁的问题。弗兰克–康登近似的优点在于能够大大减小计算量，得到可靠的计算结果。电子跃迁的结果可以使用位型坐标图来表示，它们代表了缺陷能级和价带顶之间的电子跃迁结果，氧填隙的光吸收和发射的一维位型坐标图(Configuration coordinate)如图4所示。

从图4中不难看出， $E_{abs}=E_{ZPL}+\text{Δ}{E}_e$ ， $E_{lum}=E_{ZPL}-\text{Δ}{E}_g$ 。费米能级的位置设置在价带顶时，在图3(a)中，缺陷能级获得价带之间发生空穴跃迁( ${\text{<msup> O <mo> ′ </mo> </msup>}}_{\text{I}}+{\text{h}}^{\text{.}}\to {\text{O}}_{\text{I}}$ )，图3(b)中，类似的过程同样发生在缺陷能级与价带之间发生空穴跃迁( ${\text{<msup> O <mo> ″ </mo> </msup>}}_{\text{I}}+{\text{h}}^{\text{.}}\to {\text{<msup> O <mo> ′ </mo> </msup>}}_{\text{I}}$ )。

计算得到位型坐标之后，可以定量计算氧填隙归一化的吸收谱和发射谱。钼酸铅晶体的吸收和发射光谱展示在图5中，当费米能级设置在导带底时，可以看出氧填隙 ${\text{<msup> O <mo> ′ </mo> </msup>}}_{\text{I}}$ 和 ${\text{<msup> O <mo> ″ </mo> </msup>}}_{\text{I}}$ 的吸收峰位置在1.73 eV (717.5 nm)和2.15 eV (576 nm)，发射峰的位置则在1.23 eV (1008.5 nm)和1.97 eV (629.4 nm)。其中576 nm的吸收峰非常接近580 nm [20] 的实验值。这样的结果表明，580 nm的实验吸收峰由很大的可能性来自于二价填隙氧 ${\text{<msup> O <mo> ″ </mo> </msup>}}_{\text{I}}$ 。另外的结果表明，一价氧填隙的吸收和发射波长应该处在红外区。

(a) (b)

Figure 4. Configuration coordinate of ${\text{<msup> O <mo> ′ </mo> </msup>}}_{\text{I}}$ and ${\text{<msup> O <mo> ″ </mo> </msup>}}_{\text{I}}$

图4. 氧填隙 ${\text{<msup> O <mo> ′ </mo> </msup>}}_{\text{I}}$ 和 ${\text{<msup> O <mo> ″ </mo> </msup>}}_{\text{I}}$ 的位型坐标图

我们的计算结果中，缺陷形成能量已经得到很好修正，同时使用HSE方法计算得到与实验相一致的带隙值，在计算光谱的过程中，考虑了电声子耦合作用，得到与实验值相一致的结果。因此，我们预测，钼酸铅晶体中580 nm的吸收峰来自于氧填隙缺陷。我们计算的吸收和发射光谱描述了氧填隙缺陷空穴的抓捕和释放的过程，可以作为模型来期望更深入的了解白钨矿材料中的电荷捕获和电子跃迁的过程。

(a) (b)

Figure 5. The absorption and luminescence spectra of ${\text{<msup> O <mo> ′ </mo> </msup>}}_{\text{I}}$ and ${\text{<msup> O <mo> ″ </mo> </msup>}}_{\text{I}}$

图5. 氧填隙 ${\text{<msup> O <mo> ′ </mo> </msup>}}_{\text{I}}$ 和 ${\text{<msup> O <mo> ″ </mo> </msup>}}_{\text{I}}$ 的发射谱和吸收谱

4. 结论

采用DFT方法，结合FNV修正方法计算得到了准确的缺陷形成能，使用杂化泛函HSE对带边进行了修正。计算了氧填隙的吸收峰在1.73 eV (717.5 nm)和2.15 eV (576 nm)，发射峰的位置在1.23 eV (1008.5 nm)和1.97 eV (629.4 nm)，我们预测580 nm的吸收带来自于二价氧填隙。考虑到氧的填隙缺陷会大量存在于钼酸铅晶体中，所以氧填隙是钼酸铅晶体中绿色吸收光的主要来源，我们期待进一步深化对这类材料光学性能的理论研究，为实际应用提供更多的指导。

致谢

感谢刘廷禹教授在文章撰写过程中的指导。感谢小组成员在文章调查过程中的帮助。

基金项目

本研究没有任何公共、商业或非盈利资助机构的任何具体资助。