非晶固体中的振动不均匀性及其几何特征

《APS Open ScienceNEW》：Vibrational inhomogeneity in amorphous solids and its geometric features

　　

由于缺乏平移对称性，非晶固体的动力学响应预计是不均匀的，这与晶体不同。在本研究中，我们采用格林函数方法结合模拟数据中的海森矩阵来表征二维非晶固体的振动不均匀性。我们定义了振动不均匀性的序参量——频率分辨参与率，对于晶体而言，该值在所有频率范围内均为常数1，而在非晶固体中，在整个频率范围内约为0.7，这表明在大部分频率范围内振动不均匀性的程度相似。非晶固体中振动不均匀性的空间分布通过三种方式呈现：按频率、在粒子层面以及通过渗流分析，均表明振动不均匀性呈碎片化状态并在粒子尺度上发生。结果表明，格林函数方法是探测振动不均匀性的便捷工具，并为非晶固体中振动不均匀性的作用提供了另一种视角。

I. 引言

晶体的平移对称性确保了在平衡状态下所有粒子对振动的贡献相等（当每个晶胞只包含一个粒子时）[1]。对于非晶固体，这种平移对称性不存在，因此粒子可能不会同等贡献于振动，这取决于它们的局部环境。非晶固体中的振动不均匀性因其与低频局部动力学异常（如局域化模[2–4]、软斑[5]、动力学缺陷[6,7]和动力学激发[8,9]）的关联而引起了极大的关注。这些局部动力学异常也被认为与全局动力学异常（如玻色峰[6,8,9]）相关。

在文献中，振动不均匀性的识别主要基于动力学矩阵的本征模态分析。例如，当某个粒子的位移矢量在本征模态中大于或小于其他粒子时，就可以观察到振动不均匀性。一个序参量——参与率，可以直接应用于本载模态来估计不均匀性[2,4,10,11]，其定义如下：
(1)
这里，N是粒子总数，μ表示本载模态的索引，α表示自由度（DOF）的索引。参与率的上述定义面临两个问题，使其不能完美地表征振动不均匀性。第一个问题是本载模态物理意义的模糊性。例如，如前所述，由于对称性，晶体中的动力学响应应该是均匀的，但本载模态是强度在空间上波动的平面波。因此，晶体本载模态的参与率不是恒为1，而是一个分数。这一点将在后续章节中进行更详细的讨论。更重要的是，由于强度波动，无法从单个本载模态中获得动力学响应的正确空间分布。第二个问题是该定义对数值的敏感性。如果有M个粒子的强度为10，另有M个粒子的强度为1，则参与率为2M/(M+M)；然而，如果M个粒子的强度变为100而不是10，该值将急剧下降至2M/(100M+M)。这种敏感性使得该序参量无法很好地估计有多少粒子对振动有贡献。因此，要全面理解非晶固体中的振动不均匀性，仍然需要空间分布信息。

格林函数方法在非晶系统的理论分析中更为常用[12–14]，因为它能方便地表征外部激发下的动力学响应。从数学上讲，格林函数是动力学方程在施加均匀脉冲激发至每个粒子时的解[15,16]
(2)
这里，H是动力学矩阵，G(t)是时域格林函数，I是单位矩阵。因此，格林函数的对角元素G_αα(t)给出了在均匀脉冲激发下第α自由度的局部动力学响应。如前所述，对于非晶固体，尽管激发相同，但由于局部环境的变化，局部动力学响应可能因粒子而异。因此，格林函数方法提供了一种具有明确物理意义的直接探测振动不均匀性的工具。

在本研究中，我们采用格林函数方法探测模拟得到的二维（2D）非晶固体的振动不均匀性。我们定义了总振动态密度（DOS）、局部DOS（代表频率分辨振动不均匀性）和频率分辨参与率来表征非晶固体的动力学性质。我们发现，在大部分频率范围内，振动不均匀性保持在相似的水平，约为0.7。振动不均匀性的空间特征是碎片化，并发生在粒子尺度上。研究结果提示我们应从更广泛动力学和结构的角度关注非晶固体中的不均匀性，并跨越不同的频率尺度和时间尺度。

II. 模拟与理论

为了通过模拟获得非晶固体结构，我们退火了4096个双分散粒子，混合比为1:1，尺寸比为1.4:1，在二维空间中通过谐振势...相互作用（当r < r_c时，否则为零）。所有粒子具有相同的质量m，并在正则系综中施加周期性边界条件，固定面积分数φ。（长度、质量、能量和时间分别用小粒子直径σ和1表示的约化Lennard-Jones单位表示，其中ε=1。）系统首先在T=1下平衡50000时间单位，然后以ΔT=0.01逐步冷却至T=0.1，每步平衡时间τ=5000，设定冷却速率τ_c（范围从10^-2到10^-1）。更多模拟细节见文献[17]。

在总势能最小的构型E_min处，通过计算关于每个自由度α的二阶导数可以得到海森矩阵：H_αα = ?2E/?x_α?x_α。这里，自由度包括所有粒子的x和y坐标。一旦知道系统的海森矩阵H，就可以使用格林函数（G）方法获得动力学响应性质。它们在频域中的关系可以通过对公式(2)进行傅里叶变换直接获得，即
(3)
其中γ是确保因果性的小数，-1表示矩阵的逆，I是单位矩阵。总态密度可以容易地获得如下[18]
(4)
其中Tr计算矩阵的迹。α自由度上的局部DOS为
(5)
为了方便显示玻色峰，DOS曲线通常根据Debye定律进行归约后绘制：S = g(ω)/(ω2)和S = g(ω)/(ω^(d-1))。

公式(5)中定义的局部DOS提供了每个粒子在给定频率下对总DOS贡献程度的直接度量，如下恒等式所保证
(6)
局部DOS的优势在于揭示了特定频率下动力学响应的空间分布。然而，在我们看来，这种方法在表征局域化模[2–4]、软斑[5]甚至拓扑缺陷[19]等局部动力学对象方面相对有限。

为了识别给定频率ω下局部DOS分布的形态特征，我们定义频率分辨参与率如下[8]
(7)
频率分辨参与率P(ω)=N表示所有粒子都贡献于频率ω处的模态，因此该模态是扩展的且动力学均匀的，而P(ω)=1表示只有一个粒子贡献，因此该模态是局域化的且强烈动力学不均匀的。公式(7)中频率分辨参与率所捕捉的物理与公式(1)中本载模态参与率的物理不同。频率分辨参与率表征给定频率下振动强度的分布，即振动不均匀性，而本载模态参与率表征单个本载模态的振动不均匀性。更一般地，P(ω)不能简单地从本载模态参与率（P_μ）中获得：
(8)
图1的底部面板清楚地显示了这一不等式。在给定频率下，可能存在许多本载频率近简并的本载模态。原则上，它们对该频率振动的相对贡献是不同的。然而，公式(1)中的归一化丢弃了不同本载模态的权重信息。另一方面，频率分辨参与率不依赖于任何本载模态分解，因此能够正确地表载给定频率下振动响应的不均匀性。

图1. 左上图（晶体情况）和右上图（非晶固体情况）面板显示了从格林函数方法获得的归约态密度S（灰色）和从海森矩阵本载频率分布获得的S（灰色）。在右上图面板中，玻色峰频率被确定为约ω_BP。底部两个面板显示了晶体（左）和非晶固体（右）情况下的参与率，通过公式(1)的方法（红点）和公式(7)（黑线）获得。底部面板中的灰色虚线是本载模态的频率分箱参与率。

III. 结果与讨论

A. 格林函数与本载模态分析结果的比较

低频声波是晶体中的集体激发。由于平移对称性，在每个声波模态中，所有粒子都平等参与。因此，晶体的态密度在低频区域应遵循Debye定律[20]，频率分辨参与率在整个频率范围内应为1。在图1的左上角面板中，我们通过比较归约总态密度与本载频率分布来验证六角密堆晶体的格林函数方法。两条归约态密度曲线在整个频率范围内趋势一致。由于公式(3)中的有限展宽参数γ，格林函数结果比本载频率统计更平滑。理论上，较小的γ使结果更接近本载频率分布，而较大的γ产生更平滑的曲线。为了在精度和视觉清晰度之间取得平衡，我们在计算中使用γ=0.01。附录更详细地展示了γ对归约态密度曲线的影响。

在图1的左下角面板中，我们比较了通过公式(1)（P_μ）和公式(7)（P(ω)）获得的两种参与率。如前所述，频率分辨参与率在所有频率上恒为1，而本载模态参与率则分布在从0.5到1的宽范围内。鉴于确保所有粒子在给定频率下振动中平等参与的平移对称性，频率分辨参与率正确地表征了参与粒子的比例。相比之下，本载模态参与率描述了单个本载模态的不均匀性，但由于其分数值以及公式(8)中的不等式，无法可靠地指示特定频率下振动的参与率。

通过比较格林函数方法和本载模态分析为晶体获得的归约态密度和参与率，我们验证了格林函数方法，并强调了通过格林函数定义的参与率在表征振动不均匀性（即参与粒子的比例）方面的优势。图1的右上角和右下角面板显示了非晶固体的归约态密度和参与率。对于非晶固体，两种方法获得的归约态密度在高于ω_BP的频率处一致，但在玻色峰附近和以下（ω_BP）出现明显偏差。在零频率附近，由于模态稀疏，本载频率统计通常存在噪声。一种常见的补救方法是通过对多个独立构型进行系综平均[21]。这里，我们注意到格林函数方法的另一个小优势：未进行系综平均获得的归约态密度曲线已经与近40个构型的平均结果一样平滑，与文献[21]中2D或3D玻璃的归约态密度相当。

对于非晶固体，总归约态密度在ω_BP处明显表现出玻色峰，远低于范霍夫峰（根据图1中的晶体结果，位于约ω_VH处）。频率分辨参与率显示出几个值得注意的特征。首先，它在玻色峰频率处显示出一个下降。在玻色峰以下，P(ω) ~ 0.6。在玻色峰以上但范霍夫峰以下的频率处，P(ω)从0.6逐渐增加到0.7，表明在此范围内大量粒子可以被激发并参与振动。0.7的值低于晶体中声波的理想值1，表明存在空间振动不均匀性（见图2）。在范霍夫峰以上，频率分辨参与率降至零，表明集体动力学响应被抑制。

本征模参与率遵循相似的趋势，但取值低得多，正如晶体情况所示，因为它仅反映单个本征模的不均匀性，而非局域态密度。

图 2. 频率 、、和 下局域态密度的空间分布，从左上到右下。颜色从深到浅表示局域态密度更强。

B. 局域态密度与振动不均匀性的分布 图 2 表明，动力学响应的空间分布在宽频率范围内呈现碎片化，且发生在粒子尺度上。总有一些粒子，数量级为 ，其动力学响应强于其余粒子。为了全面展示振动不均匀性，我们在图 3 中给出了不同频率下局域态密度值的分布。在从 到 的宽频率范围内，该分布惊人地相似。大部分局域态密度值约为最大值的 左右，而局域态密度非常强的粒子极其罕见（小于 ）。因此，绝大多数粒子主导了总态密度，包括玻色子峰（BP）。尽管少数粒子主导了振动不均匀性，但其数量太少，无法显著影响总态密度。

图 3. 不同频率下局域态密度沿数值的分布。局域态密度值已根据每个局域态密度空间分布的最大值进行了归一化。

振动不均匀性以局域态密度的非均匀空间分布为特征。与频率分辨参与率（PR）相比，图 3 提供了更详细的信息，后者是总结分布的序参量。图 3 表明，图 2 所示的局域态密度分布相当均匀，且图 2 中局域态密度较大的粒子（浅色点）占少数。

C. 玻色子峰涉及大多数粒子 为了可视化整个频率范围内所有粒子的动力学响应，我们在图 4 中提供了局域态密度随频率变化相对于所有 4096 个粒子的热图。不均匀性可从两个维度观察到：对于给定频率或单个粒子而言，局域态密度都是不均匀的。在图 4 中，粒子索引已根据其低频（ ，玻色子峰附近）强度进行了排序。正如左下角所示，确实有一些粒子，数量级为数百，对低频总态密度的贡献“更多”。

图 4. 所有粒子的局域态密度随频率变化的关系。颜色从深到浅表示局域态密度更强。粒子索引已根据低频（ ）强度排列，如式 (10) 所示。

为了定量表征哪一部分粒子（ ）贡献了给定的比例（ ），我们进行了渗流分析以选取最具主导性的粒子，如图 5 所示。前 和 的粒子贡献了玻色子峰频率附近总态密度的 和 。一旦固定 ，贡献给定比例 总态密度的粒子比例 可通过关系式 (9) 确定。此处，积分区间 对应于玻色子峰区域。算符 表示根据 从大到小对原子索引 进行重新排序； 是原子 的第 个自由度的局域态密度。重新排序后，我们有 (10)

从形态上看，贡献者粒子在空间中分布相当均匀，无论它们在极端渗流阈值下形成小簇 [8,9]，还是在更宽松的阈值下形成连通的复杂模式 [6]。在这种情况下，声称低频总态密度由少数粒子主导是没有依据的。

图 5. 对总态密度低频（玻色子峰附近）贡献最大的粒子的渗流分析。渗流阈值，从左上到右下，分别为总态密度在频率区域 的前 和 贡献。渗流后存活者（涂黑部分）的比例分别为 和 。渗流在式 (9) 中进行了定量描述。

上述结果确认了我们模拟的二维非晶固体中存在振动不均匀性。尽管频率分辨的振动不均匀性（即局域态密度）决定了总态密度的特征（如玻色子峰和范霍夫峰），但其空间分布并未表现出明显的特征来区分总态密度的不同区域。考虑图 1 中的晶体情况：从动力学或总态密度角度来看，范霍夫峰和低频区域是截然不同的，因为后者对应声波而前者不对应。这种区别无法通过动力学响应强度的分布来捕捉，因为在晶体中，动力学响应强度在所有频率下都是完全均匀的。类似地，非晶固体中的玻色子峰仅伴随着频率分辨参与率的轻微下降（图 1）。在我们看来，鉴于前文讨论的参与率敏感性，这种相关性过于微弱，不足以暗示因果关系。此外，图 2 和图 3 表明，参与率主要由一小部分粒子决定，而这些粒子的数量太少，无法影响总态密度。这些结果表明，人们应重新审视局域振动不均匀性是否能提供一条令人满意的路径来解释全局动力学异常。

在引言中，我们提到了几种动力学异常，如局域模式 [2–4]、软斑 [5]、动态缺陷 [6,7] 和动态激发 [8,9]。大多数此类异常是通过本征模分析识别的，需要分离类声子和非声子贡献。这些异常通常涉及少数行为与其他粒子不同的粒子。本工作并不质疑此类动力学异常（即振动不均匀性）的存在；然而，它们对总态密度特征（如玻色子峰）的贡献应予以慎重考量，因为这些特征实际上源于大多数粒子。

IV. 结论 总之，我们利用格林函数方法结合模拟数据，更好地探测了非晶固体中动力学响应的空间分布。提出了序参量——频率分辨参与率，用于表征特定频率下振动不均匀性的程度。通过三种互补的方法，即空间图、统计分布和给定频率范围内的渗流分析，确认并可视化了振动不均匀性。结果表明，虽然振动不均匀性由少数粒子主导，但玻色子峰等全局动力学异常则由绝大多数粒子支配。这些发现为非晶固体中振动不均匀性的作用提供了另一种视角。