摘要:在0—2.0 eV的碰撞能范围内采用含时量子波包方法和二阶分裂算符传播子对初始量子态为(v₀ = 0, j₀ = 0)的D + DBr反应进行了态分辨理论水平上的动力学计算. 计算了抽取反应通道D + DBr → Br + D₂和置换反应通道$ \rm D' + DBr \to D + D'Br $的反应概率、积分截面、微分截面、产物的振动和转动分布以及速率常数等动力学性质, 并与相应的理论和实验结果进行了比较. 结果表明: 计算的速率常数与实验结果十分符合. 微分截面的结果表明: 对于抽取反应, 前向的剥离反应机制在反应过程中占据主导地位; 对于置换反应, 头碰头的反弹机制占据主导地位.
关键词: D + DBr/
反应概率/
积分截面/
速率常数

English Abstract

全文HTML

--> --> -->

氢原子和卤化氢的反应在发展和理解基础的化学反应上有着重要意义, 因此其抽取和置换反应在最近几十年当中受到了广泛的关注. 对于H + HBr及其同位素反应, 大量的实验工作集中在测量该体系的热速率常数、同位素效应、积分截面、以及产物的振动和转动电子态分布等动力学信息^[1-10]. 1989年, Aker等^[7,8]采用连贯的反斯托克斯拉曼散射光谱测量了1.6 eV能量下H + HBr反应抽取和置换反应通道的积分截面和产物的振动和转动态分布. 1990年, Umemoto等^[6]采用一种脉冲辐射共振吸附技术在214—300 K的温度范围内测量了H + HBr及其同位素反应的速率常数. 采用激光闪光光解HBr分子和连续波共振荧光检测H原子方法, Seakins和Pilling^[3]与Talukdar等^[4]也分别测量了H + HBr反应不同温度范围的速率常数.
在理论方面, H + HBr反应的势能面构建也在最近的几十年中有了极大的进展. 基于多参考组态相互作用方法(MRCI)计算的104个从头算能量点, Lynch等^[11]构建了第一个extended London-Eyring-Polanyi-Sato (eLEPS)函数形式的从头算势能面. 然而基于该势能面得到的速率常数往往低于实验值并且100多个能量点也不足以描述复杂的势能面地貌. 2003年, Kurosaki和Tayanagi^[12]采用MRCI方法和AVTZ基组构建了BrH₂体系基态和两个激发态的势能面, 并且在基态的势能面上考虑了自旋轨道耦合修正(MB1). 然而基于MB1势能面的速率常数明显大于实验值. 不久之后, 他们又对MB1势能面进行了修正, 降低了约0.3 kcal/mol的势垒高度, 并将新的势能面命名为MB2^[13]. 之后又在MB2势能面的基础上修改了Br + H₂渐进区T型结构的势阱, 并命名为MB3. Kurosaki和Tayanagi^[12]的势能面得到了广泛的应用. Quan等^[14-16]采用MB1, MB2和MB3势能面和含时量子波包方法进行了一系列的动力学计算并报道了反应概率、积分截面和速率常数等物理量. 结果表明对于H + HBr反应, MB3势能面明显优于MB1和MB2势能面, 基于MB3势能面的速率常数和实验结果十分吻合. Fu和Zhang^[17]和Fu等^[18]也采用MB3势能面和含时量子波包方法进行了动力学计算, 并报道了量子态分辨的积分和微分截面. 然而, Fu和Zhang^[17]基于MB3势能面得到的速率常数明显大于实验值. 这与Quan等^[14-16]得到的结果有很大的差异.
2011年, Jiang等^[19]构建了H + HBr反应的一个新的势能面. 采用MRCI方法进行了从头算计算, 在计算中考虑了自旋轨道耦合修正并把基组外推到了完全基组极限. 此外, 基于新构建的势能面, Jiang等^[19]采用含时量子波包方法进行了动力学计算. 结果表明基于该势能面得到的速率常数与实验结果十分吻合. 之后, Xie等^[20]采用该势能面进行了H + HBr反应的态-态动力学计算, 报道了反应概率、积分截面、微分截面、产物的振动和转动态分布等动力学信息, 并与相应的实验结果进行了对比. 结果表明, 基于该势能面得到的产物振动和转动态分布与实验值存在较大差距. 同年, Jiang等^[21]又构建了关于BrH₂体系的透热势能面并计算了Br(²p_3/2, ²p_1/2) + H₂反应的积分截面. 2012年, Xie等^[22]采用Jiang等^[21]构建的透热势能面在量子态分辨的理论水平上计算了Br(²p_1/2) + H₂反应的猝灭过程. 2012年Zhang等^[23]也构建了BrH₂体系的透热势能面并与Jiang等^[21]的透热势能面进行了比较, 同时对D + DBr反应进行了动力学计算, 报道了积分截面和速率常数等动力学信息. 2014年, Zhang等^[24]采用他们构建的透热势能面进行了H + HBr反应体系的同位素计算. 将得到速率常数与实验结果进行了对比, 发现理论结果与实验结果有或高或低的偏差. 最近的势能面是由Li等^[25]在2019年构建的, 在从头算计算中采用MRCI方法、AVTZ和AVQZ基组计算了约63000个从头算能量点, 并采用神经网络方法对势能面进行了拟合. 势能面的构建过程中, 考虑了自旋轨道耦合修正并将基组外推到了完全基组极限. 此外, 基于新构建的势能面, Li等^[25]采用含时量子波包方法进行了H + HBr反应的动力学计算. 得到的抽取通道的产物振转态分布以及速率常数都与实验结果十分符合.
综上, 之前的理论研究大部分集中在H + HBr反应上, 对于同位素的理论研究相对较少. 据我们所知, 到目前为止还没有关于其同位素反应在量子态分辨理论层次上的研究. 基于Li等的势能面, 本文采用含时量子波包法和二阶分裂算符传播子对D + DBr反应进行初态为(v₀ = 0, j₀ = 0)的动力学计算.

含时量子波包方法在计算原子与双原子的反应散射以及包含多原子的散射计算中有着广泛的应用^[26-30]. 这里, 仅简单的介绍一下原子与双原子分子反应散射的动力学方法. 在空间坐标系下, 采用反应物坐标, 对于给定的总角动量J, D + DBr反应体系的哈密顿量可以写为

${\hat H} = - \frac{{{\hbar ^2}}}{{2{\mu_R}}}\frac{{{\partial ^2}}}{{\partial {R^2}}} - \frac{{{\hbar ^2}}}{{2{\mu_r}}}\frac{{{\partial ^2}}}{{\partial {r^2}}} + \frac{{{{( {\hat J - \hat j} )}^2}}}{{2{\mu _R}{R^2}}} + \frac{{{{\hat j}^2}}}{{2{\mu _r}{r^2}}} + \hat V,$

其中$\hat J$是总角动量，$\hat j$是DBr分子的转动角动量，μ_R和μ_r分别为R和r方向的约化质量.
在进行动力学计算之前, 具有确定初始振动和转动量子态的初始波包需要提前构建. 初始的含时波包是由体坐标系下的平动、振动和转动基矢构成$\left\{ {u_n^\upsilon \left( R \right){\varphi _\upsilon }\left( r \right)Y_{jK}^{JM\epsilon }\left( {R, r} \right)} \right\}$, 初始含时波包可以写为

$\begin{split} \psi _{{\upsilon _0}{j_0}{K_0}}^{JM\epsilon } = \; &\mathop \sum \limits_{n,\upsilon ,j,K} F_{n\upsilon jK,{\upsilon _0}{j_0}{K_0}}^{JM\epsilon }\left( t \right)u_n^\upsilon \left( R \right){\varphi _\upsilon }\left( r \right) \\ & \times Y_{jK}^{JM\epsilon }\left( {R,r} \right),\end{split}$

其中, 下标n为平动基矢, $ \epsilon $是体系的宇称, $ \left({\upsilon }_{0}{j}_{0}{K}_{0}\right) $是体系的初始量子态, M是总角动量在空间坐标系下在z轴的投影.
在传播过程中, 二阶分裂算符^[31]进行初始波包的传播. 假设t时刻的波函数为$ \psi (R, r, t) $, 那么经过时间步长$ {\varDelta }_{t} $, 波函数变为$ \psi (R, r, t+{\varDelta }_{t}) $

$\begin{split} &\psi \left( {R,r,t + {\Delta _t}} \right) \\ =\;& {{\rm{e}}^{\tfrac{{ - {\rm{i}}V{\varDelta _t}}}{{2\hbar }}}}{{\rm{e}}^{\tfrac{{ - {\rm{i}}{T_j}{\varDelta _t}}}{{2\hbar }}}}{{\rm{e}}^{\tfrac{{ - {\rm{i}}\left( {{T_R} + {T_r}} \right){\varDelta _t}}}{\hbar }}}{{\rm{e}}^{\tfrac{{ - {\rm{i}}{T_j}{\varDelta _t}}}{{2\hbar }}}}{{\rm{e}}^{\tfrac{{ - {\rm{i}}V{\varDelta _t}}}{{2\hbar }}}}\\& \times \psi \left( {R,r,t} \right),\end{split}$

其中${T}_{j}=\dfrac{{(\hat{J}-\hat{j})}^{2}}{2{\mu }_{R}{R}^{2}}+\dfrac{{\hat{j}}^{2}}{2{\mu }_{r}{r}^{2}}$, $ {T}_{R}=-\dfrac{{\hslash }^{2}}{2{\mu }_{R}}\dfrac{{\partial }^{2}}{\partial {R}^{2}} $, $ {T}_{r}=-\dfrac{{\hslash }^{2}}{2{\mu }_{r}}\dfrac{{\partial }^{2}}{\partial {r}^{2}} $, V是势能.
在计算中, 采用快速傅里叶变换^[32]计算径向动能算符, 在离散变量表象^[33]里计算势能算符. 此外, 为了满足边界条件, 在边界处设置一个吸收势, 其形式为

$D\left( r \right) = \begin{cases}\exp \left[ { - {\varDelta _t} \cdot {C_a} \cdot {{\left( {\dfrac{{r - {r_a}}}{{{r_b} - {r_a}}}} \right)}^n}} \right], & {r_a} \leqslant r \leqslant {r_b},\\\exp \left[ { - {\varDelta _t} \cdot {C_b} \cdot {{\left( {\dfrac{{r - {r_b}}}{{{r_{{\rm{end}}}} - {r_b}}}} \right)}^n}} \right]\exp \left[ { - {\varDelta _t} \cdot {C_a}} \right], & {r_b} \leqslant r \leqslant {r_{{\rm{end}}}}.\end{cases}$

最终, 通过反应物坐标变换方法(RCB)^[34,35]来提取态-态散射信息. 利用态-态散射矩阵可以得到反应概率、积分截面和微分截面等动力学信息.

$P_{{\upsilon _0}{j_0}{K_0}}^J\left( E \right) = \frac{\hbar }{{{\mu _R}}}{\rm{Im}}\left[ {\left\langle {\psi \left( E \right){\rm{|}}\delta \left( {r - {r_0}} \right)\frac{\partial }{{\partial r}}{\rm{|}}\psi \left( E \right)} \right\rangle } \right],$

${\sigma _{{\upsilon _\upsilon }{j_\upsilon } \leftarrow {\upsilon _0}{j_0}}} = \frac{{\text{π}}}{{\left( {2{j_0} + 1} \right)k_{{\upsilon _0}{j_0}}^2}}\mathop \sum \limits_{{K_\upsilon }} \mathop \sum \limits_{{K_0}} \mathop \sum \limits_J \left( {2J + 1} \right){\left| {{S}_{{\upsilon _\upsilon }{j_\upsilon }{K_\upsilon } \leftarrow {\upsilon _0}{j_0}{K_0}}^J} \right|^2},$

$\frac{{{\rm{d}}{\sigma _{{\upsilon _\upsilon }{j_\upsilon } \leftarrow {\upsilon _0}{j_0}}}\left( {\vartheta ,E} \right)}}{{{{{\rm{d}}\varOmega }}}} = \frac{1}{{2{j_0} + 1}}\mathop \sum \limits_{{K_\upsilon }} \mathop \sum \limits_{{K_0}} {\left| {\frac{1}{{2ik_{{\upsilon _0}{j_0}}^2}}\mathop \sum \limits_J \left( {2J + 1} \right){d}_{{K_\upsilon }{K_0}}^J\left( \vartheta \right){S}_{{\upsilon _\upsilon }{j_\upsilon }{K_\upsilon } \leftarrow {\upsilon _0}{j_0}{K_0}}^J} \right|^2},$

其中$ \psi \left(E\right) $是通过$ \psi \left(t\right) $的快速傅里叶变换得到, $ {{S}}_{{\upsilon }_{\upsilon }{j}_{\upsilon }{K}_{\upsilon }\leftarrow {\upsilon }_{0}{j}_{0}{K}_{0}}^{J} $是散射矩阵, $ {{d}}_{{K}_{\upsilon }{K}_{0}}^{J}\left(\vartheta \right) $是转动矩阵, $ \vartheta $是散射角.
指定初始量子态的速率常数是通过热平均积分截面的平动能得到, 其公式为

${k_{{v_0}{j_0}}}\left( T \right) \!=\! \sqrt {\frac{{8{k_{\rm{B}}}T}}{{{\text{π}}\mu }}} {\left( {{k_b}T} \right)^{ - 2}} \cdot \! \int_0^\infty \!\! {E{{\rm{e}}^{ - E/{k_{\rm{B}}}T}}} {\sigma _{{v_0}{j_0}}}\!\left( E \right){\rm{d}}E,$

其中$ {k}_{\mathrm{B}} $是玻尔兹曼常数.

本文的动力学计算是基于2019年Li等^[25]报道的势能面. 图1给出了D + DBr反应抽取反应通道和置换反应通道的最小能量路径. 由图1可知, 在抽取反应通道入口处有一个微小的势垒, 其高度约为1.54 kcal/mol. 该反应是一个放热反应, 其放热能约为18.77 kcal/mol, 并且在反应路径上没有明显的势阱. 对于置换反应通道, 在反应路径上有一个高度为12.14 kcal/mol的势垒.
图 1 势能面的最小能量路径　(a) 抽取反应通道; (b) 置换反应通道
Figure1. Minimum energy path of the potential energy surface: (a) Abstraction channel; (b) exchange channel.

在进行动力学计算之前需要对计算参数进行调试. 本文测试了不同参数J = 0的反应概率, 并最终得到了最优的收敛参数. 计算中使用的参数如下:
格点范围和大小 $R \in \left[ {1.5, 20.0} \right], {N_R} = 179$, $ {N}_{R}^{\mathrm{i}\mathrm{n}\mathrm{t}}=159 $, $ r\in \left[\mathrm{1.5, 20.0}\right] $, $ {N}_{r}=179 $, $ {N}_{r}^{\mathrm{a}\mathrm{s}\mathrm{y}}=49 $ $ {j}_{\mathrm{m}\mathrm{i}\mathrm{n}}=0—{j}_{\mathrm{m}\mathrm{a}\mathrm{x}}=180 $, $ {N}_{j}=180 $;
初始波包 $\mathrm{exp}\Big[-\dfrac{{(R-{R}_{c})}^{2}}{2\varDelta _{R}^{2}}\Big]\mathrm{c}\mathrm{o}\mathrm{s}\left({k}_{0}R\right)$, ${k}_{0}= $ $ \sqrt{2{E}_{0}{\mu }_{R}} $, $ {R}_{c}=11.0 $, $\varDelta _{R}=0.15$, $ {E}_{0}=1.0\;\mathrm{e}\mathrm{V} $;
投影截面位置 $ R=12.0 $;
总角动量J在z轴投影数目 $ {K}_{\mathrm{m}\mathrm{a}\mathrm{x}}=9 $;
总的传播时间 $ {T}_{\mathrm{t}\mathrm{o}\mathrm{t}}=13000 $;
传播时间步长 $ {T}_{\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{p}}=10 $.
2

-->

3.1.反应概率

在0—2.0 eV的碰撞能范围内, 总角动量J = 0, 20, 40, 60, 80, 90的反应概率见图2. 通过Li等^[25]的势能面可知, 对于抽取反应和置换反应对应的势垒高度分别为1.54 kcal/mol (约0.067 eV)和12.14 kcal/mol (约0.525 eV). 由图2可知, 对于抽取反应, 其阈值接近为0; 对于置换反应, 其阈值约为0.47 eV; 明显低于势能面上的势垒高度, 这是由于无论是抽取反应还是置换反应都发生了隧穿效应, 从而降低了反应阈值. 对于抽取反应, 所有J的反应概率的外形都很类似. 都是先快速升高当峰值接近0.1左右之后随着碰撞能的升高反应概率降低. 这是因为当能量较低的时候只能通过隧穿效应发生反应, 随着能量升高, 隧穿效率增加反应概率增高. 而当碰撞能大于势垒高度, 巨大的放热能使得D原子的速度变快, DBr分子来不及与D原子形成化学键. 因此, 当碰撞能升高时, 反应概率反而降低, 这也是放热反应的典型特征. 对于置换反应, 反应概率随着碰撞能的升高而增高, 计算的最大J = 90, 由图2可知, 在研究的碰撞能范围内所有的结果都是收敛的. 然而, 对于抽取反应只有在碰撞能低于1.2 eV时, 得到积分截面和微分截面的结果才是收敛的.
图 2 若干总角动量J的反应概率　(a) 抽取反应; (b) 置换反应
Figure2. The reaction probabilities of several selected total angular momentum J: (a) Abstraction reaction; (b) exchange reaction.

2 -->

3.2.积分截面

抽取反应的总积分截面和振动态分辨的积分截面分别列入图3(a)和图3(b). 为了与其他理论结果比较, Zhang等^[24]的结果也列入到了图3(a). 如图3(a)所示, 本文的结果与Zhang等^[24]的结果十分符合尤其是当碰撞能比较低的时候, 然而随着碰撞能的增加, 本文的结果与Zhang等^[24]的结果差距逐渐变大. 这可能是由于计算中使用了不同的势能面. 本文采用的是2019年Li等^[25]最新报道的绝热势能面, 该势能面没有考虑不同电子态之间的非绝热效应. 而Zhang等^[24]结果是基于他们自己构建的非绝热势能面. 我们认为高能部分的差距是由于非绝热效应引起的. 图3(b)列入了积分截面的振动分辨, 发现了一个有趣的现象, 当碰撞能低于0.45 eV时, 出现了振动粒子数反转的现象. 产物D₂的分布主要集中在v' = 1上. 随着碰撞能的升高, 反转现象变得不明显进而消失. 这种反转分布符合Polanyi的规则^[36,37]. 随着碰撞能的升高, 其他振动激发态的通道逐渐打开. 尽管如此, D₂分子的振动态分布主要还是集中在基态和第一激发态上.
图 3 (a) 抽取反应总积分截面和文献[24]的结果比较; (b) 抽取反应振动分辨的积分截面
Figure3. (a) Comparison between present total integral cross section of abstraction reaction and the values obtained from Ref. [24]; (b) the vibrational state-resolved integral cross section of abstraction reaction.

图4给出了抽取反应在碰撞能为0.4, 0.6, 0.8和1.0 eV时产物的振转态分布. 由图4可知, 随着碰撞能的升高, 振动电子态的数目和转动量子数的数目都增加. 这反映了随着碰撞能的升高会有更多的碰撞能转化为内能. 此外, 还可以看出越高的振动态具有越少的转动量子数, 这是因为当总能固定的情况下, 内部的转动能量转化为振动能量. 我们还发现随着碰撞能的升高, 出现了共振峰结构. 例如, 当0.4 eV时, v' = 1是单峰结构, 0.6 eV时变为双峰, 0.8 eV时变为三峰, 1.0 eV时变为多峰结构. 在图1, 抽取反应的最小能量路径上没有发现明显的势阱. 本文认为, 尽管最小能量路径上没有势阱, 然而当产物D₂分子处于振动激发态的时候, 就会形成势阱, 并且随着振动态数目的变大形成的势阱也会变深, 相应的势阱支持的束缚态和准束缚态的寿命会变长并且不同电子态之间的波函数耦合也会变得强烈. 这也可以解释为什么1.0 eV时基态的共振峰的数目多于其他能量.
图 4 在若干选择的碰撞能下, 产物D₂的振转态分布
Figure4. The ro-vibrational distribution of D₂ at several selected collision energies.

图5给出置换反应的总积分截面和振动分辨的积分截面. 它们的形状与图2中的反应概率类似, 都是随着碰撞的增加单调递增. 并且随着碰撞的增加, 其他振动激发态的通道逐渐打开. 产物的振动态分布也是随着振动激发态的增加逐渐递减. 主要的振动态分布集中在基态和第一激发态上.
图 5 置换反应总的积分截面和振动分辨的积分截面
Figure5. The total and vibrational state-resolved integral cross section of exchange reaction.

图6给出了置换反应在碰撞为0.8, 1.2, 1.6和2.0 eV时的产物分布. 与图4类似, 随着碰撞能的升高, 振动态的数目和转动量子数的数目都增加, 表明更多的碰撞能转化为内能. 在0.8和2.0 eV时的基态出现了有趣的双峰结构. 1.2和1.6 eV基态的峰偏向于较大转动量子数的位置.
图 6 在若干选择的碰撞能下, 置换反应的产物的振转分布
Figure6. The ro-vibrational distribution of products for exchange reaction at several selected collision energies.

2 -->

3.3.微分截面

图7给出了不同碰撞能下抽取反应和置换反应的微分截面. 对于抽取反应, 在低碰撞能时主要有后向和侧向散射信号, 这表明在低碰撞能时, 头碰头的反弹机制反应机制占据主导地位. 随着能量的升高, 后向散射变的不明显, 前向散射变得越加明显, 这表明反应机制由头碰头的反弹机制转变为剥离反应机制. 侧向的信号在Li等^[25]的分析中是因为势能面上两个电子态交叉的势阱引起的, 是不可靠的. 对于置换反应, 在低碰撞能时, 散射信号主要集中在后向散射. 随着碰撞的增加, 后向散射信号变得明显, 同时侧向和前向的散射信号也变得明显. 不难看出对于置换反应, 反弹反应机制一直占据主导地位. 随着碰撞能的增加, D的速率变得越来越大, 使其带走质量相对较重的Br原子成为可能, 因此出现了侧向和前向的散射信号.
图 7 (a) 在若干碰撞能下, 抽取反应的微分截面; (b) 在若干碰撞能下, 置换反应的微分截面
Figure7. (a) The differential cross sections of abstraction reaction at several selected collision energies; (b) the differential cross sections of exchange reaction at several selected collision energies.

2 -->

3.4.速率常数

抽取反应和置换反应在不同温度范围内的速率常数分别示于图8(a)和图8(b)中. 为了与相应的理论和实验结果比较, 文献[1,6,24]的结果也列入到了图8(a)中. 对于抽取反应, 速率常数随着温度的升高线性增加. 此外, 本文结果与Zhang等^[24]的非绝热结果以及实验值都十分吻合. 尽管在低温时出现了一些微小的差距. 这可能是因为本文的结果是基于绝热的势能面, 另一方面与实验结果比较需要包含转动激发态的贡献, 而本文只考虑了基态的结果. 图8(b)给出了0—3000 K温度范围内置换反应的速率常数. 如图8(b)所示, 速率常数随着温度升高而增加.
图 8 (a) 在200?1000 K温度范围内, 抽取反应的速率常数以及文献[1,6,24]的结果; (b) 在0?3000 K温度范围内, 置换反应的速率常数
Figure8. (a) The rate constant of abstraction reaction and the values obtained from Refs.[1,6,24] in the temperature range from 200 to 1000 K; (b) the rate constant of the exchange reaction in the temperature range from 0 to 3000 K.

基于Li等的势能面采用含时量子波包方法和二阶分裂算符传播子对D + DBr反应的抽取通道和置换通道进行了动力学计算. 在量子态分辨的理论层次上报道了反应概率、积分截面、微分截面、产物振转态分布以及速率常数等动力学信息. 积分截面结果与Zhang等^[24]结果十分吻合, 尽管在高碰撞能部分有一些微小的差距. 速率常数以及实验结果与Zhang等^[24]的结果进行了比较. 结果表明, 本文结果与实验结果十分吻合. 微分截面信息表明, 对于抽取反应在低碰撞能时, 反弹机制占据主导地位. 剥离反应机制在高碰撞能时占据主导地位. 对于置换反应, 反弹机制一直占据主导地位, 然而随着碰撞能升高时, 逐渐出现前向和侧向的散射信号.