非晶态物质原子局域连接度与弛豫动力学*

武振伟 汪卫华

1) (北京师范大学系统科学学院, 北京 100875)

2) (中国科学院物理研究所, 北京 100190)

(2019 年 12 月 10 日收到; 2020 年 1 月 2 日收到修改稿)

1 引 言

非晶态物质是自然界普遍存在的一种物质形态, 并广泛应用于人们的日常生活和工业生产中.非晶材料一方面具有与晶体材料类似的机械刚性,能够像一般固体一样承受一定的剪切应力, 另一方面在原子结构层面上又具有与液体类似的长程无序性, 这使得非晶态物质结构与性质的研究, 至今仍然是一个富有挑战性的开放性课题[1−3]. 非晶态(或玻璃态)通常由高温熔体的快速淬冷过程获得,在玻璃化过程中, 过冷液体的黏滞系数或弛豫时间会在很窄的温度区间内发生超过十个数量级的变化, 而其原子结构在一定程度上却没有显著变化.即玻璃转变现象并不存在传统相变理论描述中常见的对称破缺过程, 因此玻璃转变是否有结构起源一直是一个悬而未决的科学问题, 这也使得玻璃及玻璃转变的本质成为凝聚态物理和材料科学的核心难题之一, 长期吸引了众多领域科学家的研究和关注.Science创刊 125 周年之际, “玻璃态物质的本质是什么?”入选该出版社公布的21世纪亟待解决的125个最具挑战性的最前沿科学问题之一[4].

理解非晶态物质的本质及玻璃转变过程势必要对过冷液体的原子结构及结构-动力学关系有深刻了解[5]. 在描述或表征液体和玻璃的原子结构时,目前科学研究中所使用的序参量大多基于局域原子结构(最近邻原子)的几何特征给出[6,7]. 由于对种类繁多的局域原子结构的拓扑不变性认识不清,目前人们在刻画非晶态物质的结构-动力学关系时,往往只能给出定性层面的研究结果, 难以建立具有良好一一对应性的定量性质的构效关系, 从而无法对新材料的研发提供十分有力的理论支撑. 这也是至今人们对玻璃转变及过冷液体动力学非均匀性[8]是否存在结构起源这个问题存疑的重要原因, 极大地阻碍了人们对非晶及玻璃转变的本质的认识和理解. 随着研究的一步步加深, 人们逐渐意识到非晶态物质中原子中程序对体系性质的重要影响, 在构建具体体系的结构-动力学关系时, 原子中程序可能起到更加至关重要的作用[9−12]. 因此, 如何定义好的包含原子中程结构信息的序参量逐渐引起了人们的重视[13−22]. 如何有效地表征非晶态物质中的原子中程序, 以及揭示中程序与非晶态物质性能的关联 (如动力学急剧变化、热力学响应极值以及力学化学性能等), 对进一步认识非晶态物质本质和特性有重要意义. 基于局域短程序或原子团簇通过排列组合的方式搭建中程序是一种途径, 但此类方法一般构建方式繁杂, 且不易具象或数学表达, 从而对随后的结构-动力学分析造成一定困难.因此能够定义出具有更高变换不变性且易于表达的拓扑类型的序参量迫在眉睫. 将拓扑的概念引入物理学研究并不是新鲜事物, 早在20世纪70年代, 英国物理学家 Michael Kosterlitz和 David Thouless就将拓扑涡旋成对与否作为序参量来研究二维超流中的相变问题, 基于拓扑相变的概念他们随后成功发展出了Kosterlitz-Thouless相变理论, 80 年代早期, Thouless等又将拓扑观念引入电子结构研究, 与Duncan Haldane关于一维自旋链的研究工作一同导致了拓扑物态新概念的产生, 由于拓扑学原理的引入对整个物理学的发展产生了革命性影响, 三人被授予2016年诺贝尔物理学奖[23]. 由此可以推断, 在非晶态物质的结构与性质研究中尝试引入具有拓扑性质的序参量将具有重要意义. 这种性质的序参量首先要包含一定原子中程结构上的信息, 其次要具有更广义的不变性来归纳统一复杂多变的局域原子结构[24]. 在此基础上的序参量可对液体及非晶材料的原子结构进行更加简明有效的分类, 为后续建立更具定量意义的结构-动力学关系提供基础.

综上所述, 新的方法有望进一步简化对非晶结构的描述, 为更好地探索和理解非晶态物质的本质及特性打开一扇新的大门. 本文将在第2部分简要介绍非晶态物质研究的基本理论和计算模拟方法;第3部分集中综述原子的局域连接度这一新的原子结构序参量的引入, 及其在理解金属过冷液体诸多动力学行为上的应用; 第4部分做简单总结展望.

武振伟, 北京大学力学与工程科学系理学博士, 北京大学物理学院量子材料科学中心博雅博士后, 现就职于北京师范大学系统科学学院. 曾获博士研究生国家奖学金和北京大学优秀博士学位论文奖等奖励. 长期从事非晶态物理领域的研究工作, 以计算机模拟结合统计力学、图论等研究方法, 探索过冷液体和玻璃态物质等典型无序系统的结构与结构动力学关系. 代表性工作为具有拓扑属性的原子结构描述参量的提出及其构效关系研究.

汪卫华, 中国科学院院士, 博士生导师, 主要从事新型大块非晶及其他亚稳材料的制备、结构、物性研究, 研制出多种厘米级的大块非晶合金等亚稳材料, 并系统研究了其微结构及物理性能, 建立了一系列有自己特色的制备和研究亚稳材料的方法; 近年来, 致力于非晶合金材料的高通量设计、制备及表征工作, 在国际著名期刊Science,Nature,Nature Materials,Physical Review Letters,Advanced Materials,Physical Review B,Applied Physical Letters等上发表科研论文 500余篇, 在Adv.Mater.,Prog. Mater. Sci.,Mater. Sci. Eng. R,Acta Mater.等上发表邀请综述性文章6篇, SCI引用超过两万次.

2 非晶合金原子结构及动力学行为简述

理论上, 非晶态物质结构的研究主要借助于计算机模拟技术及一些常用的平衡态或非平衡态统计物理方法. 非晶局域原子结构的表征通常采用Voronoi多面体分析法, 该方法能给出无序体系中每个原子的最近邻原子结构特征, 是对非晶材料整体原子结构进一步分析和研究的基础. 非晶态动力学行为的表征主要通过计算原子的各类时空关联函数, 如能反映出系统结构弛豫快慢的中间散射函数等. 研究表明, 二十面体原子团簇在非晶态物质中广泛存在, 并被用作一种重要的结构单元来描述和表征非晶合金的原子结构. 在玻璃化转变过程中, 二十面体团簇的比份含量会随着液体的过冷而迅速增加, 并与体系的玻璃形成能力[25]、 热力学稳定性及力学行为[26]之间存在密切关系[6]. 因此, 在目前的研究中, 具有局域五次对称性[7,27,28]和空间平移非对称性的二十面体结构是一类普遍接受的非晶合金原子结构序参量. 进一步的研究还表明,在金属过冷液体中, 二十面体原子团簇具有非常强烈的自聚集效应[13], 这一特性与其他种类的局域原子团簇有着明显不同[29]. 二十面体原子团簇会通过自聚集效应形成一些链状或网状的中程尺度或大尺度空间原子结构, 这些非局域的原子结构便形成了非晶材料中的一类中程序[13]. 这类中程序在金属玻璃中相当普遍, 可将其作为中程序的一类典型代表. 最近又有研究结果表明, 在金属过冷液体中, 二十面体团簇中程序的建立和空间结构特征, 对系统的结构弛豫行为及热稳定性都有重要影响[30−33]. 目前, 非晶态材料结构与动力学关联的研究还主要是从短程序(最近邻原子结构特征)分析入手, 多侧重于对系统整体行为的影响的描述. 实际上, 非晶态物质在原子结构上“隐含”有中程序的存在[14−18,20], 但人们对无序系统中程序的表征及其对玻璃转变和玻璃材料性能的影响的认识还很缺乏, 特别是微观层面上的结构-动力学关联, 目前的研究大多数还停留在对中程序进行探索性描述的阶段. 因此, 如何有效探究和表征非晶态物质原子结构的中程序, 揭示中程序与非晶材料物性的关联, 如动力学弛豫[34], 热力学响应极值[35−37], 以及力学、化学性能等[38−40], 对进一步认识非晶态物质的本质和特性有重要意义.

2.1 分子动力学模拟

计算模拟技术的发展为研究液体及玻璃态物质的原子结构提供了非常实用有效的方法和途径.它不仅可以对样品中的原子结构进行简单明了的3D可视化, 并且还能系统深入地研究原子结构与其微观动力学行为之间的对应关系. 即使在某些情况下计算机模拟会受到时间和空间上的限制, 但只要所选取的研究对象合理, 并确保数值模拟结果不会受到时间或尺寸效应的影响, 依据所得的模拟数据, 仍然可以对诸多科学问题做出合理有效的讨论.

方法简介分子动力学(MD)模拟的基本原理是基于经典牛顿力学来预测体系中各个粒子在未来时刻的运动轨迹. 在任意时刻t, 当系统中各个粒子的位置和速度信息被给定后, 如果可以计算出作用在每个粒子上的力, 那么就可以运用牛顿方程来得到在t+∆t时刻各个粒子的位置、速度以及加速度信息, 这个过程随着时间的演化一直重复下去, 就得到了系统中各个粒子的运动轨迹. 在模拟计算过程中, 温度、压强一类的物理量可以通过具体的统计力学方法由所有粒子的位置、速度或加速度信息计算得出. 平衡态物理假设下这些物理量的时间平均值等价于其系综平均值. 分子动力学模拟中每个粒子的受力分析至关重要, 其计算结果的可靠性和准确性依赖于所选取的描述系统中粒子之间相互作用的势函数(potential). 在金属合金体系的模拟中, 嵌入原子势(EAM potential)的应用最多, 它可以相当好地描述和刻画非晶合金的结构及其系统内部的各种微观动力学过程. 针对在模拟工作中研究较多CuZr体系, 常用的EAM potential由美国Ames实验室开发和维护, 相应的势函数文件可从LAMMPS软件的potentials文件夹中找到, 或者从网站 http://www.ctcms.nist.gov/potentials/Cu-Zr.html下载得到. 目前这套势函数有2007版和2009版两个版本, 并已有专门的学术论文[41]对这一体系的不同势函数做了详细讨论.从讨论结果来看, 07版在对体系动力学性质的描述上与实验结果更加一致(如对自散射函数、不同波矢下的结构弛豫时间以及体系中Cu原子的自扩散系数的描述, 都与实验结果符合得非常好); 在对原子结构的描述上, 07版势函数倾向于促成体系中二十面体原子团簇的形成, 这在一定程度上使得系统的缓慢动力学行为更加符合实验结果的预期. 综合各方面的因素[41], 本文研究中所采用的势函数版本为07版[42].

样品制备非晶态物质的数值模拟常用的分子动力学模拟软件为开源代码LAMMPS[43]. 针对模拟工作中常见的CuZr系统, 研究体系的大小一般为10000个原子, 使用NPT系综和周期性边界条件, 模拟时间步长为 1 fs或 2 fs. 模拟的初始准备工作一般为: 1) 将目标研究比例的Cu原子和Zr原子随机分布于一定体积的模拟盒子中, 或者将满足原子比例要求的Cu晶体和Zr晶体置于模拟盒子当中; 2) 进行能量极小化过程, 一般使用共轭梯度退火算法, 这一过程是为了避免在模拟一开始由于可能的原子位置过近而产生的巨大原子间排斥力而使得系统在演化开始即崩溃而导致模拟中断; 3) 给系统中的各个原子赋予满足初始温度(如 2000 K)要求的初始速度, 一般按照高斯分布或玻尔兹曼分布给定, 然后启动分子动力学模拟,让模拟样品在高温液态弛豫足够长的时间, 一般为 2—4 ns, 这个时间相比于系统在 2000 K 下的结构弛豫时间(< 0.1 ps)是非常长的. 这样便得到了研究体系在高温液态下的原子构型, 这是后续一切模拟过程的出发点. 当要在某个特定的温度下对体系的性质展开研究时, 可将上述过程所得到的高温液态样品通过缓慢的降温过程淬火到目标温度,这个步骤可以为连续降温, 降温速度一般为0.01—1 K/ps不等; 也可为阶梯式降温, 即以一定的降温速率冷却到一定温度后, 让系统先行在这个温度下弛豫一段时间后再进行进一步的降温并重复这个过程. 一般后者所用的降温方法可以得到能量更低的样品构型. 当系统演化至目标温度下时,让样品构型在这个温度点上弛豫一段时间后再进行接下来的数据采集或原子构型分析, 弛豫等待时间一般要比系统在这个温度下的结构弛豫时间要长, 以保证一些平衡态的统计规律可以适用, 而不会使分析结果受到体系老化效应(aging effect)的影响.

2.2 原子结构表征

静态结构因子由于非晶态物质的原子结构不再具有长程有序性, 所以晶体材料中类似晶格常数的参量不再适用于描述非晶态原子结构. 但仍然可以从统计的角度出发, 探索式地认识非晶态物质中原子的短程序乃至中程序, 并以此为基础, 逐步建立一些针对非晶材料的结构-性能关系. 此小节将对在非晶态物质结构研究中经常用到的静态结构因子作简要介绍和讨论.

对于一个体积为V且含有N个粒子的系统,体系密度的微观表达式为其对应的静态结构因子则定义为? 在傅里叶空间的关联函数

其中(2)式对于无序结构与晶体结构同样适用. 一般计算结果可参照图1.

图1 静态结构因子S(q)的数值计算结果二维示意图(a)无序结构; (b)有序晶体结构Fig. 1. Structure factor S(q) obtained in computer simulation: (a) For disordered structure; (b) for ordered crystal structure.

图2 液态 Cu50Zr50 在 1000 K 下的静态结构因子 S(q) (图中对不同计算方式得到的结果进行了对比)Fig. 2. Structure factors S(q) of liquid Cu50Zr50 at 1000 K,obtained with different protocols.

这里对这两种不同的计算方式所得到结果在图2中做了对比, 可以看出, 两种计算方法所得到的的峰/谷位置基本相当, 但还是会在峰/谷的具体强度上略有差异, 这种差异在小q的区间上尤为明显, 这在一定程度上表明, 由定义(2)式通过原子构型直接计算所得能够更好地显示出非晶态原子结构的中程有序性.

复杂网络与图论分析本文综述的部分研究内容与图论(Graph theory)在描述和表征非晶态物质原子结构上的应用相关, 因此这里就图论做一点简单介绍. 图论在处理复杂网络相关问题时得到了广泛应用[44], 本文所介绍的部分研究内容就借鉴了其中的一些方法和定义来处理非晶合金中局域原子团簇连接及空间堆积情况中的问题.

网络(数学上称为图, Graph)是由一些基本元素所构成的集合, 其中一般包含顶点(或称节点)以及它们之间的连接关系. 数学上, 图G由一个非空的有限顶点集合V和伴随的一个集合E组成, 其中集合E由集合V的二元子集构成. 集合V中的元素成为顶点, 集合E中的元素成为边. 这里简单介绍两个数学上关于图G的术语: 图G的势(Cardinality of G)定义为集合V中的元素个数; 某个节点的度 (Degree of)定义为与节点有边相连的其他节点的总个数. 由定义可知, 可以用图的势来刻画具体问题当中的某个团簇的尺寸大小, 而节点的度可以用来描述体系中一个粒子周围的特征环境信息. 例如, 图3所示的一张简单的图, 它的势为 5, 其中节点的度为 2, 相应的公式表达为 Card(G) = 5, Deg() = 2.

当一个图很简单时(如图3), 可以通过观察来获取许多重要信息, 比如节点个数、连接情况及各个节点的度等[44]. 但当处理的网络涉及成千上万个节点以及错综复杂的连接关系时, 必须诉诸于某些数学上的统计方法, 并借助于计算机来处理如此庞大的数据. 因此, 需要将图转化为计算机更易处理的矩阵表示. 图的矩阵表示称为图的邻接矩阵(adjacency matrix), 邻接矩阵A的构造非常简单,对于无向无权图[44],A的元素, 当节点i和j之间有连接时为 1, 无连接时为 0, 且有由于不考虑节点自身与自身之间的连接, 所以总有= 0. 综上所述可知, 针对本文所关心的研究体系, 相应网络的邻接矩阵A总是主对角线元素为0的对称矩阵(图3). 从图3的邻接矩阵表示中, 可以简单地得到关于节点的一些重要信息, 比如对邻接矩阵的第i列或j行的元素进行求和, 便可以得到相应的Deg(i)或Deg(j).

图3 图或网络的示例以及它的邻接矩阵[44]Fig. 3. A small example graph (or network) with its adjacency matrix[44].

2.3 动力学行为表征

均方位移均方位移(mean square displacement, MSD)是描述粒子相对初始位置的偏移量随着时间演化的物理量, 其定义如下:

散射函数中间散射函数 (intermediate scattering function, ISF)[45−47]与密度-密度关联函数的傅里叶变换相关, 其定义为

它给出了体系中粒子在波矢q所给定的空间尺度内的所有动力学信息. 当公式中的求和部分有时, 公式退化为该函数的非相干部分, 又称为非相干散射函数 (incoherent intermediate scattering function)或自散射函数(self-intermediate scatter-ing function, SISF), 通常标记为一般, 在q较小时,比衰减得要慢 (图5(a)), 但当q较大时,和之间的差异并不是十分明显(图5(b)).

图4 CuZr合金中 Cu 原子和 Zr原子在 1000 K 下的均方位移Fig. 4. Mean-square displacement for Cu/Zr atoms in liquid CuZr alloy at 1000 K.

图5 相干散射函数 F(q, t)和自散射函数 Fs(q, t) 在不同q 下的比较 (a) q = 0.6 Å－1; (b) q = 2.8 Å－1Fig. 5. Comparison between F(q, t) and Fs(q, t) at different q values: (a) q = 0.6 Å－1; (b) q = 2.8 Å－1.

2.4 原子振动谱分析

原子振动谱的一种获得方式是进行正交模分析 (normal-mode analysis), 可以通过微小位移扰动 (small displacement perturbation, 一般针对势函数表达式非常复杂的体系)或解析的办法(体系势函数表达式相对简单)首先得到体系的Hessian矩阵[48], 该矩阵的表达式为

图6 针对CuZr玻璃的不同粒子数及不同计算方法所获得的振动态密度D(w)的比较, 及对玻色峰的展示[49]Fig. 6. Vibrational density of states for CuZr glass with different protocols, and the test for the present of a boson peak[49].

图7 CuZr金属玻璃10 K下通过对其原子的速度自关联函数进行傅里叶变换所得到的振动态密度, 数据表明10 K下系统内部已经几乎不存在可能影响到D(w)的老化过程Fig. 7. Vibrational density of states obtained by calculation of the time Fourier transformation of the velocity auto-correlation function. It can be seen that there is no apparent aging effect at 10 K.

3 过冷液体原子局域连接度与多尺度动力学关联

从液体到玻璃的转变通常伴随着急剧的动力学变化但却无明显原子结构上的改变, 这使得过冷液体及玻璃态中结构-动力学关系的研究成为凝聚态物理和材料科学里极具挑战性的开放性问题之一[4]. 研究表明, 液体在被过冷的过程中会展现出诸多异常的缓慢动力学行为, 这些异常的动力学行为可由系统对微小扰动的响应过程显现, 弛豫响应函数的衰减方式一般为拉伸e指数形式 (stretched exponential decay)[45,54,55], 这种弛豫方式通常被认为是与过冷液体的动力学非均匀性 (dynamic heterogeneity)相关[45,54,55]. 然而, 当系统进入玻璃态(非平衡态)以后, 体系中原子的输运机制会发生明显改变, 实验数据表明, 在诸如结构玻璃[56,57]、胶体悬浮液和凝胶[58,59]、或玻璃形成溶剂中的纳米颗粒[60,61]等非平衡系统中, 原子弛豫响应函数的衰减速度一般要快于e指数形式, 显示出非平衡态物质中原子所具有的独特的输运行为, 这种输运行为跟液体或过冷液体中典型的扩散或亚扩散的输运方式明显不同, 并且这种不同并不能基于动力学非均匀性的图像来加以解释. 目前, 非平衡态物质所特有的这种快弛豫模式(faster-than-exponential relaxation) 一般解释为系统内部的某些内应力释放过程[58,59]或者玻璃转变附近体系中溶剂原子的协同运动所导致[60,61], 对这一快弛豫过程在微观原子结构层面上的研究还没有深入开展. 尤其值得关注的是, 在系统整体行为表现为正常弛豫方式的近平衡过冷液体中, 是否也存在某些具有特殊局域对称性的原子同样可以表现出这种快弛豫模式, 对这一问题的探究将有助于人们对过冷液体本身性质及玻璃转变的理解. 实验上, 非弹性X射线(inelastic X-rays)衍射技术可以在一定程度上探究液体和玻璃的结构与动力学, 但探测的空间尺度和实验精度并没有到达能够细致刻画微观动力学过程的程度.因此, 能时刻跟踪记录原子位置和状态的分子动力学模拟技术的应用, 是对探究过冷液体在不同空间尺度上结构-动力学关系的有效理论补充.

此外, 人们在探索非晶态物质的时间关联函数时发现的另一个非常有趣的现象是: 针对某些非晶形成体系 (主要指 strong glass former), 在其时间关联函数的微观动力学部分之后, 笼蔽效应(caging effect)所导致的弛豫平台之前, 会出现一个小的谷信号 (a dip or valley). 有研究称[62−66]这个小的阻尼振荡行为在某些情况下可能与频域空间的玻色峰(Boson peaks)相关. 在频域空间中,玻色峰一般表现为光、X射线、或中子散射探测[67,68]中的一些非弹性散射过程. 最近的研究还发现, 非晶材料中慢的结构弛豫过程与快的玻色峰频率附近的原子振动之间存在一定的相关关系[69−72], 这就使得关联函数中短时间尺度上的小的阻尼振荡(dip)跟长时间尺度上反常的快弛豫模式或超扩散现象是否存在一定的原子结构上的联系成为一个根本且十分有趣的问题. 如果弛豫函数中的短时阻尼振荡可以跟原子的局域结构建立联系, 那他们之间的依赖关系也是同样值得探讨的话题. 因此, 对短时阻尼振荡及长时快弛豫行为的结构-动力学关系的研究, 将有助于揭示非晶态物质中各尺度动力学行为存在关联关系的微观原子结构基础.

本节将针对Cu50Zr50这一常见的非晶合金模型体系, 对新近的一些微观原子层面上的结构-动力学关系的研究和讨论进行简要介绍. 对Cu50Zr50中一类最近邻具有二十面体特殊配位的原子(icosahedrally coordinated atoms, IC atoms) 进行相干/非相干散射函数的测量研究发现, IC原子的自散射函数曲线中, 存在位于微观弹道过程与笼蔽效应平台之间的阻尼振荡行为(damping dip),并且, 这一阻尼振荡发生时间和振荡幅度都与研究粒子的原子结构性质-局域连接度密切相关[49]. 与此同时, 随着原子局域连接度的上升, 其相应的自散射函数的长时弛豫方式逐步从拉伸e指数衰减转变为压缩e指数衰减, 表明过冷液体中原子结构的变化也可导致其从扩散到超扩散的动力学行为的转变[49]. 本节所综述的研究内容说明, 在非晶态物质中, 原子的长时弛豫动力学行为和短时振动激发行为都可与原子的微观结构特征建立紧密联系.

3.1 原子局域连接度的分析与统计

Voronoi分析法是表征无序系统中局域原子结构的一般方法, 它可以给出指定原子所对应所有最近邻原子的配位方式和对称性的信息[73]. 本文的主要讨论对象IC粒子所对应的Voronoi指标为分析表明, Cu50Zr50在 1000 K 时体系中IC粒子的数量约占原子总数的3.67%, 其中3.6%为Cu原子. 所以绝大部分的IC粒子的化学元素种类为Cu. 针对IC原子, 其局域连接度(local connectivity)的定义为: 与指定IC原子有直接“连接”关系的其他IC原子的个数(图8), 这里的“连接”规定为两个IC原子互为最近邻[30,31]. 也就是,如果想要知道一个指定的IC原子的连接度, 那么就要考察其12个最近邻原子中有多少个原子同样属于一个二十面体原子团簇的中心原子. 可以想象, 原子的连接度其实是一种非绝对局域(包含部分中程原子结构信息), 并具有拓扑属性的结构序参量. 非局域性体现在: 序参量要对中心原子的所有最近邻原子的最近邻进行考察, 也即它部分包含中心原子的次近邻结构信息; 拓扑属性体现在: 它不关心团簇具体的几何结构如何, 也不关心团簇连接的具体空间位置如何, 而只关心其周围可能的相互连接数目的多少. 一般意义下, 一个好的结构序参量应尽量定义简单, 并能尽可能多地包含一些重要的物理信息, 比如能够自动地内涵“长度”或“角度”的信息等. 长度信息一般与体系的密度(或局域密度)密切相关, 比如在大家所熟知的一级相变中, 体系密度的变化在相变描述上扮演绝对重要的角色[35,36]; 角度信息一般与局域原子结构或体系整体的对称性, 亦或是某个物理量的相位信息有紧密联系, 例如在描述超流(superfluidity)或者电子结构的非平庸拓扑性质的时候, 波函数的相位信息成为影响材料物性的关键因素[23]. 研究表明, 原子连接度的变化会相应地伴随有局域原子团簇体积的变化[75]和团簇对称性[31]的变化, 因此它是一种可以同时反映出原子结构在长度和角度上的变化的序参量. 图9给出了IC原子的局域连接度的分布[49].可以看出, 1000 K 时局域连接度为k= 0 和 1 的粒子的数目最多, 同时对比不同温度下的数据可以发现, 局域连接度的分布具有温度依赖性, 这也从一定程度上反映出体系的原子结构随温度变化的演化[49].

图8 局域连接度定义的示意图(指定原子的最近邻原子中, 与指定原子具有相同局域对称性(用相同的颜色表示)的原子的总数即为指定原子的连接度)Fig. 8. Illustration of the definition of particles with different connectivities k : Particles in blue are the center of an icosahedral-like cluster.

3.2 局域连接度与短时弛豫动力学

首先给出了具有不同局域连接度k的IC原子的自散射函数

其中求和针对所有具有相同k值的IC原子,是粒子j在t时刻的位置矢量. 计算时q的模值取为2.8 Å–1, 约等于体系静态结构因子主峰的位置.如图10所示的静态结构因子由(2)式通过样品原子构型直接计算得出, 可以看出,(特别是峰位)与体系温度的依赖关系并不十分强烈[49].

图11给出了具有不同k值的IC原子的自散射函数数据表明, 在短时尺度内每一个k值所对应的曲线都会在 0.2 ps左右存在一个小幅振荡. 如前所述, 一般情况下这是 strong glass former中才有的特征短时时间关联函数行为, 但这里的研究[49]发现, 即使是在Cu50Zr50这样的 fragile glass former中, 某些具有特殊局域对称性的粒子(如这里的IC原子) 同样可以展现出此类特殊的短时动力学特征. 另外从图11中的数据可以看出, 体系中所有原子(All atoms)、Cu 原子 (Cu atoms)、以及其他非 IC-Cu原子(NON-IC-Cu)的曲线并没有展现出这个特征的短时动力学信号, 这说明这个谷信号是体系中IC原子所特有的一种动力学弛豫模式[49]. 这个特征动力学信号还与IC原子的局域连接度k值密切相关. 随着k值的增加, 振荡信号的谷底位置向高频移动, 且伴随有振荡幅度的增加. 当粒子运动的时间尺度很小时, 粒子主要在其初始位置附近做微小振动, 这种微小振动的行为特征显然与体系的振动简正模式(vibrational normal-mode)密切相关. 因此计算了具有不同k值IC原子的振动态密度 (density of states, DOS),

图10 Cu50Zr50 在不同温度下的静态结构因子 S(q)[49]Fig. 10. The q-dependence of the partial structure factors for three temperatures considered[49].

为了能够更好地理解k值与短时动力学之间的构效关系, 基于简正模式分析的有关结果,Wu等[49]在研究中提出了一个唯象的拟合公式.基于用两个谐振子的线性叠加来描述粒子运动的基本思想,在短时区间上的表达式便可解析表达为

图11 具有不同 k 值的粒子的自散射函数 (SISF)[49], 图中黑色实线为唯象模型((9)式)的拟合结果, 图的右上角给出SISF曲线的整体形状; 其他类型粒子的SISF曲线也一并在图中给出, 以方便对比Fig. 11. Short-time behavior of the self-intermediate scattering function of particles with different local connectivity k(symbols)[49]. The wave-vector is q = 2.8 Å–1 and T =1000 K. The solid lines are fits to the data with Eq. (9).Also included is Fs(q, t) for the Cu atoms in an icosahedral cluster (dashed red line), the Cu atoms not in an icosahedral cluster (blue dashed line), and all Cu atoms (green).The black dashed line is the correlation function averaged over all atoms. The upper inset shows the same data in a larger time interval.

图12 具有不同 k 值原子的振动态密度 D(w)Fig. 12. Vibrational density of states of particles with different local connectivity k.

图13 (a) wH 和 wL 与局域连接度 k 均成正相关关系 [49]; (b) 随着 k 值的增加, 权重因子 CL 上升而 CH 降低; (c) 高频模式wH(q)的波矢q无关性表明其局域模式的特征; (d) wL(q)所具有的色散关系表明它是一种扩展性质的模式, 图中黑色虚线为相应数据点的线性拟合Fig. 13. (a) Both the high and low frequency modes, wH and wL, increases with increasing k[49]; (b) the fraction of motion CL/H increases for wL and decreases for wH; (c) the high frequency mode wH(q) is approximately q-independent, characteristic of localization of the vibrational modes; (d) the low frequency mode wL(q) increases monotonically with increasing q, characteristic of collective dynamics.

3.3 局域连接度与长时输运行为

本节将主要针对IC原子的长时动力学弛豫行为做简要分析和讨论. 图14(a)给出了1000 K下具有不同连接度k的IC粒子在波矢q= 2.8 Å–1时的自散射函数长时弛豫曲线[49]. 这一阶段的时间关联函数可用Kohlrausch-Williams-Watt (KWW)公式来拟合,

图14 (a) 具有不同 k 值的粒子在 q = 2.8 Å–1 下的长时动力学弛豫曲线, 图中黑色实线为 KWW 公式拟合所得[49]; (b) 不同波矢q下的形状因子b与k值的关系, b值的变化预示着动力学行为从拉伸e指数弛豫到压缩e指数形式的转变, 转变的有无及具体位置与波矢q密切相关Fig. 14. (a) Long-time decay of the correlation functions at q = 2.8 Å–1 for particles with different k values[49]. The black solid lines are the Kohlrausch-Williams-Watt (KWW) fits. (b) The k dependence of the exponent b. The variation of b reveals a dynamic crossover from stretched (b < 1) exponential relaxation to compressed (b > 1) one. It can be seen that the cross-over from stretched to compressed exponential depends on q.

图15 弛豫时间 t 与波矢 q 之间的关系, 为了更好地区分 1/q scaling 和 1/q2 scaling, 这 里 把 数 据 重 新 表述 成 了qt与q之间的关系[49]Fig. 15. Wave vector q dependence of the relaxation time t of the final decay of Fs(q, t) for particles having different local connectivities[49]. Here we show qt as a function of q to make it simpler to see the 1/q law and to distinguish it from the 1/q2 law.

如前所述, 压缩e指数式的结构弛豫行为常见于非平衡系统当中, 与体系中粒子的超扩散输运方式呈密切相关. 在玻璃态物质中, 这种超扩散过程的微观机制一般认为与系统中的局域内应力释放有关[58]. 而黏稠玻璃形成液体中的纳米颗粒的超扩散行为则一般认为与溶剂在玻璃转变过程中的原子协同运动相关[60,61]. 因为在颗粒悬浮系统当中, 纳米颗粒之间的间距一般比较大, 并有没承载内应力的结构基础. 针对这里所讨论的CuZr过冷液体, 研究认为内应力释放和协同运动这两种机制对高k值IC原子的超扩散行为都有一定贡献. 一是研究表明高k值的粒子的运动模式中扩展模成分增加(图13(b)), 为前面提到的可能的周围介质对高k值原子产生的等效对流传送作用提供了一定的数据支持; 二是已有研究表明[75]局域连接度增加会导致粒子的自由体积收缩并伴有局域模量上升, 这为局域应力集中提供了物理上的条件. 在应力集中区域周围活跃粒子流的等效剪切作用下,应力释放所伴随的比较激烈的原子弛豫过程便可发生, 进而原子快弛豫模式的动力学行为就可能会被观测到.

图16 1000 K 下标记为不同类型的原子的均方位移曲线[49]Fig. 16. Mean squared displacement for different type of atoms at 1000 K[49].

分子动力学模拟结果表明[78], 对非晶样品进行略微低于玻璃转变点的保温退火处理可获得最大的退火效率. 通过此种方法可获得能量状态较低的样品且样品中的二十面体原子团簇的占比会明显增加. 二十面体团簇具有自聚集效应[13], 大量二十面体团簇的出现同样意味着高局域连接度原子的出现, 这些高k值IC原子的超扩散行为对体系平均弛豫行为影响和贡献可能就会变得越来越显著. 最新的实验研究结果表明[79], 对金属玻璃进行0.9左右的退火预处理, 那么在后续的应力弛豫实验中就会相应地检测到超扩散行为的信号, 这种信号的出现与退火预处理过程密切相关.

综上所述, 在近平衡过冷液体中, 不同的原子结构可导致非常丰富的异常动力学行为的出现. 即使是在 fragile liquid 当中, 受特殊对称性“保护”的原子同样可以出现一般在strong liquid中才有的动力学行为. 粒子的短时动力学行为与输运方式都可与新的结构序参量-局域连接度建立联系. 局域连接度或可成为非晶材料中短时动力学行为与长时动力学行为常常具有关联的微观原子结构基础.从另一角度讲, 同一体系同一温度下由不同原子结构所引起的动力学行为的转变, 也可为理解过冷液体和玻璃材料中结构-动力学关系提供新的途径.

4 总 结

本文首先阐述了非晶态物质研究中的几个关键的科学问题, 并简要介绍了目前基于分子动力学模拟和非平衡统计的主要研究方法, 诸如无序系统局域原子结构的表征、动力学时空关联函数的计算、以及简正模式分析等. 为更好地理解非晶态物质中的结构-动力学关系, 在以上研究方法的基础上, 本文介绍了图论或复杂网络方法在描述非晶态物质原子结构及构建新型构效关系上的应用.Wu等[30]定义了粒子的连接度这一新的结构序参量, 这个序参量不同于以往基于局域几何结构所给出的原子结构描述参量, 从数学上图论里的一般理解出发, 局域连接度是具有拓扑性质的序参量, 它包含一定的中程原子结构信息. 基于这一新的结构描述参量的系列研究结果表明: 非晶态物质中原子的连接度不仅对粒子的动力学快慢[30]、局域原子团簇对称性及形成能(formation energy)[31]等诸多物理性质产生重要影响, 而且还在过冷液体的热稳定性[32]、晶化路径[80]、及异常形核速率[81]等系统整体行为的理解上起到关键作用. 就本文主要综述的研究内容[49]而言, 局域连接度不仅会调制原子的短时动力学弛豫行为, 而且还会对原子的长时输运方式产生重要影响. 局域连接度是一种可能的联系非晶态物质中多尺度动力学关联[69−72,82]的微观结构基础. 随着原子局域连接度的升高, 短时振动模式的振动频率上升, 振动模式的扩展成分增加, 局域成分减小; 与此同时, 在一定的空间尺度内, 原子的长时动力学弛豫行为从以扩散方式为主导逐步转变为以超扩散为主导, 相应的弛豫函数的衰减形式也会发生从拉伸自然指数到压缩自然指数的变化. 这是超扩散动力学行为第一次在近平衡过冷液体中被观测到并被给予详细描述, 表明过冷液体中某些具有特殊局域环境的原子同样可以展现出以往只在非平衡系统原子身上才能观测到的特殊动力学行为. 本文所综述的研究内容在一定程度上丰富了人们对非晶态物质中原子结构的描述方式及过冷液体中某些特殊原子的反常动力学行为的认识.