温稠密氢氘混合物热力学、 光学和电子输运特性理论研究

刘蕾 张伟 李治国

摘要：温稠密状态下的氘氚混合物是惯性约束聚变必然要经历的中间物质状态. 其在温稠密状态下的热力学特性和输运参数在聚变实验设计研究中起着重要作用. 但由于氚具有放射性，实验中可以用氢氘混合物来模拟替代研究氘氚混合物. 本文采用了密度泛函理论分子动力学模拟和基于自洽流体变分理论的化学模型，构建了氢氘混合物在0～200 GPa和300～ 10 000 K 宽范围状态方程和声速数据库. 此外，结合Kubo-Greenwood公式计算得到了氢氘混合物电子电导率以及光学反射率和折射率，建立起混合物绝缘-金属化相变与实验上易探测的光学物理量反射率和折射率之间的联系，为从实验角度探索氢及其同位素的金属化提供了理论指导.

关键词：温稠密物质； 氢氘混合物； 分子动力学； 状态方程； 声速； 电导率

中图分类号：  O414.12  文献标识码：A  DOI：10.19907/j.0490-6756.2023.054003

收稿日期：  2023-07-13

基金项目：  国家自然科学基金（12204387， 12174357）； 四川省科技计划 （2023NSFSC1370）； 西南科技大学博士基金（21zx7113）

作者简介：   刘蕾（1992-）， 男， 四川南充人， 博士研究生， 特聘副教授， 主要从事极端条件下物质特性的研究. E-mail： lei_liuchn@163.com

Theoretical study on the thermodynamic， optical， and electronic  transport properties of warm dense hydrogen-deuterium mixtures

LIU Lei  1 ， ZHANG Wei  1 ， LI Zhi-Guo  2

（1.School of Mathematics and Physics， Southwest University of Science and Technology， Mianyang 621010， China;

2.Institute of Fluid Physics， Chinese Academy of Engineering Physics， Mianyang 621900， China）

The deuterium-tritium mixture under warm dense state is an intermediate material state that is inevitably encountered in inertial confinement fusion. The thermodynamic properties and transport parameters of this mixture play a significant role in the design of fusion experiments. However， due to the radioactivity of tritium， hydrogen-deuterium mixtures can be used as a substitute in experiments to simulate and study deuterium-tritium mixtures. In this study， density functional theory molecular dynamics simulations and a chemical model based on self-consistent fluid variational theory were employed to construct the equation of state and sound velocity database for hydrogen-deuterium mixtures in the range of 0～200 GPa and 300～10 000 K. Additionally， using the Kubo-Greenwood formula， the electronic conductivity， optical reflectivity， and refractive index of hydrogen-deuterium mixtures were calculated. The relationship between the insulator-to-metal phase transition in the mixture and the experimentally detectable optical physical quantities， including optical reflectivity and refractive index， was established. This theoretical guidance helps in exploring the metalization of hydrogen and its isotopes from an experimental perspective.

Warm dense matter; Hydrogen-deuterium mixture; Molecular Dynamics; Equation of state; Sound speed; Electronic conductivity

1 引 言

温稠密物质是介于凝聚态物质和热稠密等离子体之间的一种物质状态，其粒子数密度在10  22 ～10  25 /cm  3 ，温度在0.1～100 eV范围，非理想耦合参数1 <Γ<100  ［1］ . 在过去的几十年里，温密物质（Warm Dense Matter， WDM）一直是科学界的热门课题，在惯性约束聚变（Inertial Confinement Fusion， ICF）、重离子聚变和天体演化等方面具有较强的工程应用背景和重要的科学意义，因为它是上述过程物质存在和发展必然要经历的中间物质状态. 处于温稠密状态下的物质将会出现一系列的物理化学反应，如分子解离、原子复合、电子激发、分子及原子电离等，形成包括分子、原子、离子、电子以及团簇在内的高度瞬态的复杂混合体，其往往呈现出离子间强耦合和电子部分简并的特性. 对跨越如此宽温度和压力范围的复杂量子统计系综性质的描述需要准确刻画非线性、强关联、强耦合、部分简并、多体和量子等多种物理效应，因此传统的固体理论和等离子体理论都难以有效应用  ［2］ .

温稠密氢其及其同位素氘、氚在理解极端条件下的物质行为方面发挥着关键作用，在惯性约束聚变  ［3］ 和天体物理学  ［4］ 中有着广泛的应用. 温稠密氢、氘、氚状态方程（Equation of State， EOS）及其输运特性是辐射流体动力学模拟的输入参数，对理解和设计惯性约束聚变实验至关重要  ［5］ . 同时，由于温稠密氢是气态巨型星的主要成分，这些参数也将有助于了解巨行星的形成、内部结构和演化规律  ［6］ . 一些基于动高压技术的实验，包括轻气炮  ［7-9］ 、磁驱动  ［10］ 、激光驱动  ［11］ 和爆炸驱动  ［12］ ，被用于获得氢或氘单质的状态方程、电导率和光学反射率. 但是实验不仅代价昂贵，而且获得的实验数据是离散的，无法实现以较小代价获得能用于建模的宽温度和压力区间数据. 理论计算可以实现这一目的，但其准确性和可靠性需要实验数据进一步检验. 基于两体势的经验模型无法准确描述温稠密物质中的多体效应，而模拟弱耦合和强简并等离子体的经典方法不能正确地估计由电子极化引起的非线性屏蔽  ［13］ . 基于密度泛函理论的分子动力学（Density Functional Theory Molecular Dynamics， DFT-MD）对电子进行完全量子力学处理，为计算局部热力学平衡中的电子结构和离子构型提供了一种强大的方法，是大家公认目前能够较为准确描述温稠密物质中复杂相互作用的一种方法  ［5］ .

与纯体系相比，混合物通常比其单质成分具有更丰富的相图，并且只有在这些多组分体系中才能观察到更加丰富的相互作用  ［14，15］ . 因此关于混合物的研究激发了当今实验和理论物理学家的浓厚兴趣. 特别地，我们专注于氢和氘的混合物，因为其是可以通过考虑同位素效应来模拟和理解ICF内爆特性的理想替代材料  ［5］ . 但目前关于该混合物的实验和理论研究还非常少. 2004年，中国工程物理研究院陈其峰等人  ［16］ 借助于二级轻气炮对这种混合气体进行了冲击实验，并将样品单次冲击至140 MPa. 2012年，他们使用相同的加载工具并结合多次冲击压缩技术，将该混合物的EOS扩展到36 GPa  ［17］ . 2019年，他们又在同一发实验中同时获得了该混合物一、二次冲击压缩的状态方程、折射率、反射率和分子极化率  ［9］ . 在理论方面，Kress等人  ［18］ 基于DFT-MD和自由轨道分子动力学（OFMD）获得了氘氚混合物的离子输运性质，包括扩散系数和粘滞系数. 刘蕾等人  ［19］ 利用DFT-MD研究了等比例温稠密氢氘混合物中的质子交换效应和离子输运特性，并以DFT-MD的计算结果为基准检验了单组份等离子体模型. 本文基于DFT-MD开展了大规模分子动力学模拟构建了温稠密氢氘混合物在压力范围0～200 GPa和温度范围300～10 000 K宽区域热力学特性数据库，并结合Kubo-Greenwood公式计算得到了该混合物的光学反射率和折射率以及电子电导率.

2 理论计算方法

本文使用CP2K/QICKSTEP程序  ［20］ 进行密度泛函分子动力学模拟. CP2K采用高斯波和平面波的杂化波方法，该方法允许对大系统进行快速和精确的密度泛函理论计算. Goedecker-Teter-Hutter（GTH）赝势  ［21］ 与双ζ价极化基组（DZVP）一起被采用. 平面波截止能量被设置为500 Ry，而高斯基组的截断能量设定为40 Ry. 自洽场（Self-Consistent Field， SCF）猜测使用始终稳定的预测-校正器（ASPC）方案，SCF收敛准则设定为10  -5   a.u ..电子交换相关势采用基于广义梯度近似（Generalized Gradient Approximation， GGA）的Perdew-Burke-Ernzerhof（PBE）泛函  ［22］ . 众所周知，GGA-PBE泛函无法准确描述体系中的色散力，这限制了其在低密度下的范德华（vdW）系统中的应用  ［23］ . 因此，我们在这里采用非局域交换相关泛函vdW-DF1  ［24］ . 所有的模拟都采用了周期边界条件，并使用Gama点对布里渊区域进行采样. 模拟是用正则系综进行的，离子温度由Nosé-Hoover恒温器控制  ［25］ ，电子温度通过费米-狄拉克分布考虑  ［26］ ，将电子温度（ T  e  ）设置为与离子温度（ T  i  ）相等来维持系统处于局域热力学平衡. 使用包含128个氢和128个氘原子的立方晶胞，并且晶胞的大小由质量密度决定. MD模拟中使用的时间步长为1.0 fs，动态模拟至少持续1 ps的时间以达到热平衡，总模拟时间为10 ps. 此外，每个温度压力点的动态电导率 σ （ ω ）是借助于Kubo-Greenwood公式  ［27］ 对最后2 ps中的20个关联度较低的结构进行计算再取平均后得到的. 在电导率计算时设置了更密集的3×3×3 Monkhorst-Pack  k 点来对布里渊区进行采样. 直流电导率 σ 是通过对动态电导率取 ω SymbolnB@

0的极限得到. 基于计算的电导率再结合Kramers-Kronig关系可以计算得到介电常数，从而计算得到光学折射率和反射率.

3 结果讨论与分析

3.1 热力学特性

3.1.1 状态方程  在图1和图2中我们展示了由PBE泛函和考虑了范德瓦尔斯效应的vdW-DF1泛函计算的氢氘混合物在密度为0.1～1.7 g/cm  3 和温度为1000～12 000 K宽区域状态方程. 从图中可以看出，在密 度低于1.5 g/cm  3 时，压力并不随着温度单调升高，等容线上存在一个 （?p/?T）<0  的不稳定区域. 该不稳定区域在氢氦混合物中也曾出现过  ［28］ . 这是由于混合物中氢或氘分子解离吸收能量导致压力降低. 在低密度区域时，分子解离主要是因为温度的升高所导致，随着密度的升高，压制解离效应也逐渐明显，因而 （?p/?T）<0 区域逐渐向低温区域移动（见图中的灰色实线）. 当密度大于1.5 g/cm  3 时，该不稳定区域向温度低于1000 K区域移动，故压力随着温度单调升高.

此外，我们还应用自洽流体变分理论  ［29］ （Self-Consistent Variational Theory， SFVT）计算了该混合物的状态方程（见图3）以检验两体势的化学模型在该温度和密度区间的适用性. 沿着各等容线，压力随着温度的升高而升高，不存在上述不稳定区域，说明SFVT模型中所考虑的分子解离相变是连续的.

图4比较了三种不同方法计算的沿着四条等容线状态方程. 从图中可以发现：温度小于6000 K 且密度小于1.5 g/cm  3 时三种方法计算的压力相差较大，且vdW-DF1给出的压力是最高的. 这可能是因为在低温低密度区域范德瓦尔斯效应不可以忽略，PBE泛函和SFVT模型均没有考虑此效应导致预测的压力偏小. 当温度大于6000 K时，三种方法给出的结果趋于一致. SFVT模型虽然只考虑了两体相互作用，但在本文所考虑的温度和压力范围内与DFT-MD预测的EOS大体符合，说明多体相互作用在此温度压力区间几乎可以忽略. 我们通过DFT-MD计算得到的氢氘混合物宽区域物态方程数据不仅可以为惯性约束聚变和天体模型构建提供参数而且可以用于检验经验物态方程模型.

3.1.2 声 速  用4阶多项式将氢氘混合物的状态方程数据进行拟合便可以获得整个温度压力区域内任何一点的状态方程数据. 接着根据基本热力学关系就可以计算得到等熵声速，计算细节可参见文献［7］，计算结果如图5所示. 在研究的温度压力区间，氢氘混合物的声速在7～18 km/s之间，与相同温度压力下的纯氢声速比较接近  ［30］ . 氢氘混合物的声速呈现出强烈的规律性：声速随着密度呈线性增加，说明混合物越稠密声速越大；随着温度的升高，声速也是增大的. 氢氘混合物的声速可以为气态巨行星地震波模型的构建提供支撑  ［7］ .

3.2 光学和电子输运特性

氢在的金属化相变是科学界一直以来关注的热点. 实现氢的金属化由两条途径：一是通过静态压缩实验，二是通过动态压缩实验  ［31］ . 理论预测静高压实验在压力为400 GPa左右能够实现氢的金属化  ［32］ . 但是想要达到这一压力无疑是非常困难的. 而据报道通过动态压缩氢同时提高压力和温度，氢的金属化压力将大大降低，大约在140 GPa左右  ［8］ ，但是目前动态冲击压缩测量氢的电导率实验数据太少. 我们通过第一性原理分子动力学预测氢氘混合物在高温高压宽区域的电导率从而为动态冲击压缩实验的开展提供指导. 因为激光冲击压缩实验只能测量波阵面反射率并借助与经验公式来推断氢的电导信息. 我们通过理论计算可以建立折射率和反射率等光学信息与电导率之间的联系，进而就可以直接从反射率大致推断冲击压缩样品的电导率. 文献报道vdW-DF1方法能够相对准确描述氢的金属化相变  ［33］ ，故本文也采用此方法进行计算.

通过DFT-MD模拟再结合库伯Symbolm@@格林伍德公式可以得到不同频率（不同波长）下的混合物的电导率. 图6中展示了在氢氘混合物在固定密度和不同温度下电导率和光学反射率随频率的变化关系. 在高温高压条件下氢氘混合物的折射率是一个复数即： n ～ =n+ik . 其中 n 为折射率的实部， k 为折射率虚部. 图7为固定密度和不同温度下折射率实部和虚部随频率的变化关系.

将直流电导率定义为频率趋于0时对应的电导率. 图8是氢氘混合物在不同密度下直流电导率随温度变化关系. 从图中可以看到随着密度或温度增大，氢氘混合物的电导率都增大，说明温度和压力是使氢氘混合物金属化的两条途径. 从图8还可以看到在同一密度下，随着温度的升高，直流电导率会出现一个突增，而且电导率突增对应的温度区间与上文所提到的等容线上不稳定区间大体一致，说明分子的解离伴随着电导率的突增. 随着密度的升高电导率的突增区间逐渐向温度较低处移动，这也与等容线上的不稳定区间是一致的. 在同一密度下，随着温度的增加，电导率会趋于一平台值，说明混合物中分子已经完全解离，电导率达到饱和.

图9是氢氘混合物在不同密度下反射率随温度变化关系. 反射率取值对应于动态冲击实验中常用的诊断波长532 nm时（对应频率为2.3 eV）的数值. 反射率随密度和温度的变化趋势与电导率相似. 这说明氢氘混合物 金属化的同时伴随着反射率的增加. 氢的最小金属化电导为2000 S/cm  ［8］ ， 比较图8和9就可以推断出氢氘混合物发生金属化相变时光学反射率大致在0.4～0.5之间.

图10是氢氘混合物在不同密度下折射率实部与虚部随温度变化关系. 在密度较低时（小于0.3 g/cm  3 ），折射率实部随温度近似呈线性关系. 随着密度的升高，折射率实部先是急剧升高再逐渐趋于一不变值达到饱和. 对比图8和图10可知，氢氘混合物处于绝缘状态时折射率虚部很小，几乎可以忽略. 随着氢氘混合物的金属化，折射率的虚部逐渐增大，变得不可忽略，这说明折射率虚部的增大与氢氘混合物的金属化也是同时进行的. 氢氘混合物发生金属化相变时，折射率虚部值的大小逐渐变得与实部值大小相当. 通过比较还可以知道氢氘混合物发生金属化相变时，折射率实部值大致在3.5～4.5之间，折射率虚部在1.5～2.5之间.

4 总 结

本文借助于大规模DFT-MD模拟和SFVT理论模型，系统地研究了氢氘混合物在温稠密条件下的热力学、光学和电子输运性质. 建立了氢氘混合物在0～200 GPa和300～10 000 K宽范围状态方程、声速以及电子输运参数数据库. 研究发现DFT-MD计算的等容线上存在 （?p/?T）<0 不稳定区域，这是由于混合物中氢和氘分子的解离所导致的. 而SFVT化学模型则没有捕捉到该不稳定区域，说明该模型中所考虑的分子解离相变是连续的. 在温度小于6000 K的区域，未考虑范德瓦尔斯效应的PBE和SFVT与考虑了范德瓦尔斯效应的vdW-DF1相差较大，当温度大于6000 K时三种方法计算的结果趋于一致，说明在温度较低时范德瓦尔斯效应对EOS影响较大. 氢氘混合物的声速随密度和温度的增加而升高，声速数值在7～ 18 km/s 之间.  电导率、反射率和折射率均随着温度的升高而升高，在密度较高时（大于0.3 g/cm  3 ）， 这三个参数的数值随着温度急剧升高，随后达到一饱和值. 数值突增的温度压力区间刚好与等容线上的不稳定区域相吻合，说明该突增是分子解离和绝缘-金属化相变所导致的. 通过对比混合物的电导率和反射率以及折射率，可以推断氢氘发生绝缘-金属化对应的光学反射率和折射率实部分别在0.4～0.5之间和3.5～4.5之间.

参考文献：

［1］   Koenig M， Benuzzi-Mounaix A， Ravasio A，  et al . Progress in the study of warm dense matter ［J］. Plasma Phys Control Fusion， 2005， 47： B441.

［2］  Falk K. Experimental methods for warm dense matter research ［J］. High Power Laser Sci Eng， 2018， 6： e59.

［3］  Lambert F， Recoules V， Decoster A，  et al . On the transport coefficients of hydrogen in the inertial confinement fusion regime ［J］. Phys Plasmas， 2011， 18： 056306.

［4］  Guillot T. Interiors of giant planets inside and outside the solar system ［J］. Science， 1999， 286： 72.

［5］  Hu S X， Goncharov V N， Boehly T R，  et al . Impact of first-principles properties of deuterium-tritium on inertial confinement fusion target designs ［J］. Phys Plasmas， 2015， 22： 056304.

［6］  Swift D C， Eggert J H， Hicks D G，  et al . Mass-radius relationships for exoplanets ［J］. Astrophys J， 2012， 744： 59.

［7］  Li G J， Li Z G， Chen Q F，  et al.  Multishock to quasi-isentropic compression of dense gaseous deuterium-helium mixtures up to 120 GPa： probing the sound velocities relevant to planetary interiors ［J］. Phys Rev Lett， 2021， 126： 075701.

［8］  Nellis W J， Weir S T， Mitchell A C. Minimum metallic conductivity of fluid hydrogen at 140 GPa ［J］. Phys Rev B， 1999， 59： 3434.

［9］  Liu L， Chen Q F， Gu Y J，  et al . Measurement of multiple physical parameters of dense gaseous hydrogen-deuterium mixture under double-shock compression： evaluating theoretical models from multiple views ［J］. Appl Phys Lett， 2019， 115： 231905.

［10］  Knudson  M D， Desjarlais M P. High-precision shock wave measurements of deuterium： evaluation of exchange-correlation functionals at the molecular-to-atomic transition ［J］. Phys Rev Lett， 2017， 118： 035501.

［11］ Loubeyre P， Brygoo S， Eggert J，  et al . Extended data set for the equation of state of warm dense hydrogen isotopes ［J］. Phys Rev B， 2012， 86： 144115.

［12］ Mochalov  M A， Ilkaev R I， Fortov V E，  et al . Quasi-isentropic compression of a nonideal plasma of deuterium and its mixture with helium at pressures up to 250 GPa ［J］. J Exp Theor Phys， 2021， 132： 985.

［13］ Wang C， He X T， Zhang P. The equation of state， electronic thermal conductivity， and opacity of hot dense deuterium-helium plasmas ［J］. Phys Plasmas， 2012， 19： 042702.

［14］ Soubiran F， Mazevet S， Winisdoerffer C，  et al . Optical signature of hydrogen-helium demixing at extreme density-temperature conditions ［J］. Phys Rev B， 2013， 87： 165114.

［15］ Brygoo S， Loubeyre P， Millot M，  et al . Evidence of hydrogen-helium immiscibility at Jupiter-interior conditions ［J］. Nature， 2021， 593： 517.

［16］ Chen Q F， Cai L C， Chen D Q，  et al . The equation of state for the mixtures of dense hydrogen and deuterium ［J］. Physica B， 2004， 348： 299.

［17］ Gu Y J， Chen Q F， Zheng J，  et al . The equation of state， shock-induced molecule dissociation， and transparency loss for multi-compressed dense gaseous H  2 +D  2  mixtures ［J］. J Appl Phys， 2012， 111： 013513.

［18］ Kress J D， Cohen J S， Horner D A，  et al . Viscosity and mutual diffusion of deuterium-tritium mixtures in the warm-dense-matter regime ［J］. Phys Rev E， 2010， 82： 036404.

［19］ Liu L， Li Z G， Dai J Y，  et al . Quantum molecular dynamics study on the proton exchange， ionic structures， and transport properties of warm dense hydrogen-deuterium mixtures ［J］. Phys Rev E， 2018， 97： 063204.

［20］ VandeVondele  J， Krack M， Mohamed F，  et al . Quickstep： fast and accurate density functional calculations using a mixed Gaussian and plane waves approach ［J］. Comput Phys Commun， 2005， 167： 103.

［21］ Goedecker S， Teter M， Hutter J. Separable dual-space Gaussian pseudopotentials ［J］. Phys Rev B， 1996， 54： 1703.

［22］ Perdew J P， Burke K， Ernzerhof M. Generalized gradient approximation made simple ［J］. Phys Rev Lett， 1996， 77： 3865.

［23］ van Mourik T， Gdanitz R J. A critical note on density functional theory studies on rare-gas dimers ［J］. J Chem Phys， 2002， 116： 9620.

［24］ Dion M， Rydberg H， Schrder E，  et al . Van der Waals density functional for general geometries ［J］. Phys Rev Lett， 2004， 92： 246401.

［25］ Nosé S. A unified formulation of the constant temperature molecular dynamics methods ［J］. J Chem Phys， 1984， 81： 511.

［26］ Mermin N D. Thermal properties of the inhomogeneous electron gas ［J］. Phys Rev， 1965， 137： A1441.

［27］ Greenwood  D A. The Boltzmann equation in the theory of electrical conduction in metals ［J］. Proc Phys Soc， 1958， 71： 585.

［28］ Vorberger J， Tamblyn I， Militzer B，  et al . Hydrogen-helium mixtures in the interiors of giant planets ［J］. Phys Rev B， 2007， 75： 024206.

［29］ Chen Q， Zhang Y， Cai L，  et al . Self-consistent variational calculation of the dense fluid helium in the region of partial ionization ［J］. Phys Plasmas， 2007， 14： 012703.

［30］ Li G J， Gu Y J， Lan Y S，  et al . Compression of gaseous hydrogen into warm dense states up to 95 GPa using multishock compression technique ［J］. Phys Rev B， 2023， 107： 014309.

［31］ Nellis W J. A perspective on hydrogen near the liquid-liquid phase transition and metallization of fluid H ［J］. J Phys Chem Lett， 2021， 7972.

［32］ Dias R P， Silvera I F. Observation of the Wigner-Huntington transition to metallic hydrogen ［J］. Science， 2017， 355： 715.

［33］ Schottler M， Redmer R.  Ab initio  calculation of the miscibility diagram for hydrogen-helium mixtures ［J］. Phys Rev Lett， 2018， 120： 115703.