来源论文: https://arxiv.org/abs/2606.10101v1 生成时间: Jun 10, 2026 16:49

执行摘要

镧系金属铈（Ce）因其在温压调控下表现出极具传奇色彩的等结构 $\gamma \rightarrow \alpha$ 相变（体积坍缩近 15%）而成为凝聚态物理和强关联材料研究的圣杯体系。传统的密度泛函理论（DFT）在面对局域 $4f$ 电子所构成的多体强关联能谱时经常失效，难以定量预测其热力学平衡体积及物态方程（Equation of State, EOS）。

本研究工作实现了一次方法学上的重要突破：它将基于第一性原理的动力学平均场理论（DFT+DMFT）与精确的格点动力学（声子谱）相结合，系统性地评估了关联驱动的声子重整化作用（Correlation-driven Phonon Renormalisation）以及声子振动自由能（Vibrational Free Energy）对铈物态方程的定量贡献。通过结合最先进的嵌入式动力学平均场理论（eDMFT）、连续时间量子蒙特卡洛（CT-QMC）不纯物求解器以及**固定自能近似（Fixed Self-Energy Approximation）**下的有限位移法（Finite Displacement Method），研究成功捕捉到了 $\gamma$ 铈在强关联作用下声子模式的奇异硬化（Mode-selective Hardening）现象。

为了克服极高计算成本带来的数据稀疏性，研究人员进一步引入了基于**主成分分析结合多层感知机（PCA-MLP）的物理启发式机器学习插值方法，实现了声子色散关系和振动熵在全格点常数区间的连续高精度重构。结果表明，当同时计入局域磁矩的电子熵（Electronic Entropy）与振动熵（Vibrational Entropy）**的贡献时，计算所得的平衡晶格常数达到 5.06 Å，与实验值（5.16 Å）高度吻合，彻底解决了传统 DFT 严重高估结合能、低估体积的顽疾。这一工作标志着强关联复杂 $f$ 电子材料热力学计算步入高精度定量时代。

1. 核心科学问题，理论基础，技术难点与方法细节

1.1 核心科学问题：铈体积坍缩的物理之谜与晶格动力学缺失

elemental 铈（Ce）的 $\gamma \rightarrow \alpha$ 相变伴随着面心立方（fcc）晶格常数的剧烈收缩，但其晶体对称性保持不变（均为 $Fm\bar{3}m$）。目前学术界对这一相变的微观驱动机制主要存在两种经典理论模型：

  1. Mott 定域化-巡游化转变（Mott Transition Scenario）：认为 $4f$ 电子在 $\gamma$ 相定域，而在 $\alpha$ 相发生解定域化参与化学键合，驱动体积坍缩。
  2. Kondo 体积坍缩模型（Kondo Volume Collapse Scenario, KVC）：认为 $4f$ 电子在两相中均保持相对定域，但传导电子与 $4f$ 局域磁矩之间的 Kondo 屏蔽作用在 $\alpha$ 相（低体积、高杂化）剧烈增强，通过 Kondo 单态凝聚降低系统自由能。

定量回答上述科学问题的核心瓶颈在于：过去的研究过度关注电子自由能，而忽视了晶格振动自由能（声子）的定量贡献。实验研究对于声子在相变中的角色存在争议。非弹性中子散射（INS）显示两相间的声子态密度差异微小，而高压衍射和非弹性 X 射线散射（IXS）则表明振动熵变占总相变熵变的三分之一到二分之一（约 $0.33 \sim 0.75\, k_B/\text{atom}$）。更重要的是，晶格力常数是由处于强关联状态下的电子网络介导的。传统的静态均合场方法（如 DFT, DFT+$U$）无法捕获动力学多体自能反馈给晶格的力学重整化，因而无法给出定量正确的声子谱和物态方程。

1.2 理论基础：DFT+DMFT 与 Luttinger-Ward 泛函

为了克服上述瓶颈，必须在非微扰层次上同时处理局域电子关联和电荷自洽性。本工作采用 DFT+DMFT（密度泛函理论结合动力学平均场理论） 架构。系统总的电荷密度和局域格林函数通过电荷自洽循环确定。其热力学势（自由能）可在 Luttinger-Ward（LW）泛函形式下严格表述：

$$\Omega[G, \Sigma] = -\text{Tr} \ln \left( G_0^{-1} - \Sigma \right) - \text{Tr} \left( G \Sigma \right) + \Phi_{\text{LW}}[G] + E_{\text{dc}}$$

其中 $G_0$ 为无相互作用格林函数，$G$ 为自洽格林函数，$\Sigma$ 为动力学局域自能，$\Phi_{\text{LW}}[G]$ 为 Luttinger-Ward 泛函中所有的双外线不可约图的贡献，$E_{\text{dc}}$ 为双计数修正项（Double Counting Correction）。

双计数修正

由于 DFT 部分已经包含了一部分局域库仑相互作用，必须通过双计数修正项进行扣除。本研究采用名义双计数（Nominal Double Counting）修正方案：

$$V_{\text{dc}} = U \left( n_f - \frac{1}{2} \right) - J \left( \frac{n_f}{2} - \frac{1}{2} \right)$$

其中，设定 $n_f$ 为固定值 1.0（即名义 $4f$ 占据数），这被证实对于铈体系具有极高的鲁棒性，有效避免了电荷在 $s\text{-}d\text{-}f$ 通道之间的不真实电荷转移。

连续时间量子蒙特卡洛（CT-QMC）

不纯物求解器是 DMFT 的心脏。本工作使用**连续时间杂化展开量子蒙特卡洛（CT-HYB）**算法来求解有效单杂质 Anderson 模型（SIAM）。在倒空间虚频轴（Matsubara 频轴）上，逆温度设为 $\beta = 20\,\text{eV}^{-1}$（对应 $T \approx 580\,\text{K}$）。库仑作用参数取为 $U = 6.0\,\text{eV}$，Hund 耦合 $J = 0.7\,\text{eV}$。相互作用矩阵在不纯物求解器中采用密度-密度（Ising）形式处理。

1.3 技术难点与方法细节：力常数计算与固定自能近似

在强关联电子框架下，晶格动力学（声子）的计算面临极为严苛的算力挑战。声子谱的计算基于有限位移法（Frozen Phonon / Finite Displacement Method）。在超胞（Supercell）内，将原子偏离平衡位置移动一个微小位移 $\mathbf{u}$，通过计算系统受力 $\mathbf{F}$ 来提取力常数矩阵（Interatomic Force Constants, IFC）：

$$C_{i\alpha, j\beta} = - \frac{\partial F_{j\beta}}{\partial u_{i\alpha}}$$

在常规 DFT 中，该过程计算成本适中。但在 DFT+DMFT 中，超胞尺寸的增大会导致不纯物求解器所需的时间呈指数级增长，特别是当晶体对称性被原子位移打破后，原本等价的原子不再等价，必须在每个位移构型下重新对多个不等价不纯物进行自洽求解。

固定自能近似（Fixed Self-Energy Approximation）

为了克服这一瓶颈，研究引入了固定自能近似。该近似认为，在小振幅原子位移下，局域多体动力学相互作用（自能 $\Sigma(i\omega_n)$）的变化极其微弱，对力常数的主要贡献来自于单粒子能谱的重整化。其具体操作步骤如下：

  1. 在完美的平衡 fcc 晶格结构下，运行完整的电荷电荷自洽 DFT+DMFT 计算，获得收敛的电子自能 $\Sigma(i\omega_n)$。
  2. 将原子位移引入超胞，建立位移构型。
  3. 在计算这些位移构型的受力时，锁定不纯物自能为平衡晶格的值，只允许电荷密度和单粒子态密度在有约束的条件下自洽。这避免了在位移超胞中进行极耗时的 CT-QMC 计算。
  4. 计算超胞中的离子受力，进而计算动力学矩阵并对角化获取声子色散关系。该近似将计算开销降低了一个数量级以上，使得采用 $4\times4\times4$ 的 $q$ 点网格采样成为可能。

2. 关键 Benchmark 体系、计算数据与性能分析

2.1 晶格常数依赖的内能分析 (DFT vs DFT+DMFT)

研究团队首先对铈的基态内能 $E(a)$ 进行了系统性扫描，分析了不同理论层级的预测能力（见论文图 1）：

理论方法预测平衡晶格常数 $a_0$ (Å)实验值 $a_{\text{exp}}$ (Å)偏差幅度物理机制简析
LDA+SOC4.524.85 ($\alpha$-Ce) / 5.16 ($\gamma$-Ce)-6.8% (对 $\alpha$)局域密度近似严重过结合（Overbinding），无法描述关联局域化
GGA+SOC4.694.85 / 5.16-3.3% (对 $\alpha$)广义梯度修正部分缓解过结合，但仍极大地偏离高体积相
LDA+DMFT+SOC (Internal E)4.885.16 ($\gamma$-Ce)-5.4% (对 $\gamma$)成功捕捉到关联驱动的体积膨胀，使平衡位置向实验 $\gamma$ 相移动
GGA+DMFT+SOC (Internal E)5.285.16 ($\gamma$-Ce)+2.3% (对 $\gamma$)GGA 的起伏特征使其在 DMFT 叠加后表现出稍微过度的体积膨胀

数据表明，标准的 DFT 方法（无论是 LDA 还是 GGA）都完全无法再现 $\gamma$ 相的平衡体积，本质上是因为静态平均场方法高估了 $4f$ 电子的化学键合活性。而 LDA+DMFT 成功引入了非微扰强关联，降低了系统的结合能，使内能极小值点向大体积（高晶格常数）方向发生了显著跃迁。

2.2 声子色散关系的关联重整化行为 (DFT vs DMFT)

图 2 的全声子色散曲线和局部精细扫描（图 2b 与 图 2d）揭示了关联电子在力常数上的非凡重整化作用。研究发现了关联诱导的声子模式特异性硬化（Mode-selective Hardening）：

  • 在常规 DFT (LDA) 下的膨胀晶格中（图 2b，当晶格常数 $a$ 从 5.05 Å 膨胀至 5.25 Å）：高、中、低频声子全线崩溃式软化（Softening），在 K-$\Gamma$ 路径上表现出剧烈的下弯，这是由于典型的化学键随距离拉伸而减弱所致。
  • 在 DFT+DMFT 下的膨胀晶格中（图 2d，相同体积范围）：
    • 高频声子依旧顺从常规的膨胀软化趋势。
    • 令人震惊的现象发生在低频声子区域：位于 Brillouin 区边缘的 $X$ 点和 $L$ 点的横向声学（TA）模式出现反常硬化。例如，在 $X$ 点，随着晶格常数从 5.00 Å 增加到 5.25 Å，TA 频率竟然从 8.28 meV 逆势上升至 9.28 meV，硬化幅度超过 12%（大于 1 meV）。

物理机制深度解析

这一现象的微观物理起源非常精妙：在 DFT 框架下，部分巡游的 $4f$ 电子贡献了额外的电荷屏蔽（Screening）效应，降低了原子横向剪切运动时的离子间排斥力和力常数。而在 DMFT 框架下，强烈的局域库仑相互作用 $U$ 触发了动力学局部自能，将 $4f$ 电子的谱权重从费米能级大幅度转移开（形成相干的上下 Hubbard 带，见附录图 S4 的 DOS 谱），极大地削弱了 $4f$ 电子对剪切运动（横向位移）的屏蔽能力。因此，裸露出来的离子间硬芯排斥力使横向声子硬化。这一机制在常规静态 DFT 中是根本无法描述的。

2.3 声子自由能对 EOS 的定量修正 (Fig 6 & Table I)

通过综合电子与振动熵的贡献，研究定量划分了稳定 $\gamma$-铈的各项能量组分。论文图 6 给出了 580 K 下不同自由能组合对平衡晶格常数的修正：

  1. 仅考虑电子内能 $E$： $$\text{Min}[E] \Rightarrow a_0 = 4.88\,\text{Å}$$
  2. 考虑电子自由能 $F_{el} = E - T S_{el}$（加入局域 $4f$ 磁矩熵贡献）： $$\text{Min}[F_{el}] \Rightarrow a_0 = 4.99\,\text{Å}$$ (格点常数增加了约 2.2%)
  3. 仅在内能上加声子自由能 $E + F_{ph}$： $$\text{Min}[E + F_{ph}] \Rightarrow a_0 = 4.96\,\text{Å}$$
  4. 全自由能耦合 $F_{tot} = E - T S_{el} + F_{ph}$： $$\text{Min}[F_{tot}] \Rightarrow a_0 = 5.06\,\text{Å}$$ (与实验真实 $\gamma$-Ce 值 5.16 Å 偏差降至 1.9%)

这表明，振动熵对平衡晶格体积提供了约 1.5% 的额外膨胀贡献。如果不同时在强关联层次上自洽考虑这两者，强关联材料的物态方程定量预测就无法实现。

3. 代码实现细节、复现指南与开源工具链

为了方便科研人员复现此工作，以下梳理了其核心计算工作流、参数配置以及辅助机器学习插值的完整实现路径。

3.1 核心计算工作流拓扑

整个工作流由三个计算内核耦合而成：

[ WIEN2k (FP-LAPW 基底) ]
           │
           ▼ (电荷、轨道自洽)
[ eDMFT (不纯物求解器 CT-QMC) ] ───► 锁定收敛自能 Σ(iω_n)
           │
           ▼ (固定自能近似，进行小位移超胞计算)
[ Phonopy Wrapper ] ───► 获取受力 F_i ───► 提取动力学矩阵 D(q)

3.2 关键参数配置

(1) WIEN2k 结构参数与能带截止

  • 基底网格: $R_{\text{MT}} \times K_{\text{max}} = 7.0$，确保 LAPW 基底在紧凑 fcc 空间内的完备性。
  • 相对论修正: 必须开启完整自洽自旋轨道耦合（SOC）。而在声子受力扫描中，采用无 SOC 结构（如正文所述，SOC 对声子动力学矩阵的相对修正量 $<5\%$）。

(2) eDMFT (case.indmfl) 关联窗口选择

# 关联轨道配置: Ce 4f
1                           # 关联原子数
1   3   3                   # 原子索引, 轨道角动量 l=3, 维度
4f-correlated               # 轨道标签
-0.45  0.45  -2.5  2.5      # 能带杂化窗口选择 (eV)

(3) CT-QMC 求解器参数 (case.inldmft)

U           6.0             # 库仑排斥能 (eV)
J           0.7             # 洪德耦合 (eV)
beta        20.0            # 逆温度 (1/eV), 对应 ~580 K
nom_dc      1.0             # 开启名义双计数修正, n_f=1.0
warmup      1000000         # QMC 预热步数
max_steps   50000000        # 每个 QMC 迭代的最大 Monte Carlo 采样步数

3.3 机器学习 PCA-MLP 声子插值代码实现

由于高精度 DFT+DMFT 声子谱的计算成本巨大，研究在 $a \in [4.30\,\text{Å}, 5.25\,\text{Å}]$ 之间仅计算了 15 个离散晶格常数的点。这里提供一个完整的 Python/PyTorch 代码示例，展示如何构建基于主成分分析（PCA）降维与神经网络（MLP）逼近的声子高精度插值算法（对应论文 S9 节及图 S10）：

import numpy as np
from sklearn.decomposition import PCA
import torch
import torch.nn as nn
import torch.optim as optim

# 1. 模拟生成原始声子数据
# 假设我们有 15 个格点常数下的声子色散谱，每个谱沿高对称路径采样了 999 个点 (3个分支，共 2997 个频率值)
num_lattices = 15
num_features = 2997

# 模拟格点常数 [4.30 - 5.25]
a_samples = np.linspace(4.30, 5.25, num_lattices).reshape(-1, 1)
# 模拟声子谱矩阵 (y_raw)
np.random.seed(42)
y_raw = np.sin(a_samples) * 15.0 + np.random.normal(0, 0.1, (num_lattices, num_features))

# 2. 对声子谱执行 PCA 主成分分析
# 减去均值，提取前 K=3 个主成分 (如论文所述，主成分1占 98% 以上方差)
mean_y = np.mean(y_raw, axis=0)
y_centered = y_raw - mean_y
pca = PCA(n_components=3)
c_coefficients = pca.fit_transform(y_centered) # 形状为 (15, 3)
components = pca.components_                # 形状为 (3, 2997)

print(f"前3个主成分解释的方差比例: {sum(pca.explained_variance_ratio_):.4f}")

# 3. 构建多层感知机 (MLP) 以预测投影系数 c_k(a)
class MlpRegressor(nn.Module):
    def __init__(self):
        super(MlpRegressor, self).__init__()
        self.network = nn.Sequential(
            nn.Linear(1, 16),
            nn.Tanh(),
            nn.Linear(16, 16),
            nn.Tanh(),
            nn.Linear(16, 3) # 输出 3 个主成分的系数
        )
        
    def forward(self, x):
        return self.network(x)

# 训练模型
x_train = torch.FloatTensor(a_samples)
y_train = torch.FloatTensor(c_coefficients)

model = MlpRegressor()
criterion = nn.MSELoss()
optimizer = optim.Adam(model.parameters(), lr=0.01)

for epoch in range(1000):
    optimizer.zero_grad()
    outputs = model(x_train)
    loss = criterion(outputs, y_train)
    loss.backward()
    optimizer.step()
    if (epoch + 1) % 200 == 0:
        print(f"Epoch [{epoch+1}/1000], Loss: {loss.item():.6f}")

# 4. 连续插值与重构任意格点常数（例如 a = 5.00 Å）下的声子谱
def predict_phonon_spectrum(a_query):
    model.eval()
    with torch.no_grad():
        a_tensor = torch.FloatTensor([[a_query]])
        c_predicted = model(a_tensor).numpy() # 预测系数
    # 逆 PCA 变换
    reconstructed_spectrum = np.dot(c_predicted, components) + mean_y
    return reconstructed_spectrum.flatten()

a_test = 5.00
predicted_spectrum = predict_phonon_spectrum(a_test)
print(f"成功生成晶格常数 a={a_test} Å 下的连续声子色散关系，点数: {len(predicted_spectrum)}")

3.4 相关开源软件与库资源链接

  1. WIEN2k: http://www.wien2k.at/ (全电子 FP-LAPW 计算平台)
  2. eDMFT: http://physics.rutgers.edu/orjan/software/ (Kristjan Haule 开发的嵌入式自洽 DMFT 代码库)
  3. Phonopy: https://phonopy.github.io/phonopy/ (广泛使用的晶格动力学分析工具，提供力常数提取和热力学积分接口)
  4. Scikit-Learn (PCA): https://scikit-learn.org/ (提供降维组件)

4. 关键引用文献与局限性评论

4.1 关键参考文献及其科学贡献

  • [Ref 27] K. Haule and T. Birol, Phys. Rev. Lett. 115, 256402 (2015)
    • 贡献：首次在电荷自洽的 DFT+DMFT 层次上讨论了铈的电子自由能，确立了不纯物磁矩熵作为稳定 $\gamma$ 相主导因素的地位。本研究在其基础上引入了自洽的格点声子自由能修正。
  • [Ref 20] M. Krisch et al., Proc. Natl. Acad. Sci. U.S.A. 108, 9342 (2011)
    • 贡献：利用高分辨率非弹性 X 射线散射（IXS）精确测量了 elemental 铈在 $\gamma \rightarrow \alpha$ 相变过程中的全声子色散关系。本工作计算所得的 $\gamma$ 相声子谱直接与该实验数据进行了 Benchmark 对比（图 3），展现了近乎完美的拟合优度。
  • [Ref 16] K. Haule and G. L. Pascut, Phys. Rev. B 94, 195146 (2016)
    • 贡献：建立了在电荷自洽 DFT+DMFT 框架下直接计算作用在离子上受力的形式化数学架构，为后来的有限位移声子计算奠定了力学基础。

4.2 本工作局限性批判与科学评论

尽管本工作代表了强关联材料理论预测的最高水平，但在以下几个方向上依然存在明显的局限性：

  1. 声子计算中缺失自旋-轨道耦合（SOC）的反馈： 由于当前的 eDMFT-Phonopy 计算接口无法在超胞受力扫描中完全保持 SOC 的全电子自洽，作者在处理力常数时忽略了 SOC。虽然作者在附录（第 S5 节）中通过传统的 DFT-SOC 测试，证明了 SOC 对声子自由能的相对偏差仅在 $5\%$ 以内，但对于强关联 $f$ 电子体系而言，自旋-轨道轨道杂化（Spin-Orbit-assisted Hybridisation）是物理精细机制的重要部分，未来的全自洽研究必须将 SOC 深度整合进超胞受力计算中。
  2. Ising 密度-密度近似（Density-Density Approximation）的局限性： CT-QMC 求解器中采用的 Coulomb 相互作用算符排除了自旋翻转项（Spin-flip）和非对角杂化项。这种 Ising 近似虽然极大加快了量子 Monte Carlo 的采样速度，但在一定程度上会人为低估传导电子对 $4f$ 局域磁矩的 Kondo 屏蔽效应，特别是在中间杂化强度的区域。
  3. 固定自能近似（Fixed Self-Energy Approximation）的适用边界： 在极其剧烈的格点畸变或大振幅声子振动下，局部多体环境会发生显著重组。固定自能近似在这些情况下将会失效，无法再现强的声子-声子相互作用和极高温下的非谐（Anharmonic）物理过程。
  4. 对 U 参数的经验依赖性： 研究依然依赖于经验设定的 $U=6.0\,\text{eV}$。正如补充材料中所示，平衡晶格常数对 $U$ 的微小波动高度敏感。未来需要引入更加无偏的第一性原理方法（例如自洽的 GW+DMFT），实现相互作用参数 $U(\omega)$ 的从头计算，彻底摆脱经验参数的限制。

5. 补充内容：从铈外推到其它复杂 $f$ 电子体系的热力学展望

本研究揭示的“电子-晶格动力学强耦合”以及“关联诱导声子硬化”机制，具有极强的普适性。它不仅仅适用于 elemental 铈，更可以作为一把钥匙，打开整个锕系与镧系强关联金属及其氧化物物理机制预测的大门。

5.1 锕系金属（以钚 Pu 为例）的格点动力学挑战

钚（Pu）是元素周期表中最复杂的金属。它在常压下具有六个不同的固态相，其体积变化甚至比铈更加剧烈。其中，稳定 fcc 结构的 $\delta$-Pu 在热力学性质上与 $\gamma$-Ce 极其相似，都展现出了非平凡的局域 $5f$ 电子特征。 传统的 DFT 甚至无法给出 $\delta$-Pu 的声子稳定谱（常常预测出大面积的虚频虚模，预示着晶格不稳定性）。通过本工作确立的 DFT+DMFT + 固定自能近似 方案，科研人员将能够定量探究 $\delta$-Pu 中由于关联电子导致的横向声子非凡硬化，定量解释为什么极少量的镓（Ga）掺杂就能戏剧性地稳定原本极其脆弱的 $\delta$ 相物性。

5.2 熵协同效应的定量拆解

复杂强关联材料在有限温度下的自由能竞争，本质上是各种熵源之间的“协同与妥协”（Synergy and Competition）。通过本文的定量拆解，我们可以绘制强关联材料中总熵变的物理构图：

$$\Delta S_{\text{total}} = \Delta S_{\text{elec}} (T, V) + \Delta S_{\text{phonon}} (T, V) + \Delta S_{\text{magnetic}} (T, V)$$

在等结构相变过程中，当体积膨胀时，局部磁性自由度释放（局域矩形成），$\Delta S_{\text{magnetic}}$ 增加；与此同时，格点力常数弱化（尽管TA模式由于关联效应被抑制了软化，但总体声子谱重心依然下移），导致 $\Delta S_{\text{phonon}}$ 同样增加。本研究表明，强关联材料的相稳定化不是单一熵源驱动的，而是电子定域化伴随的磁矩解冻（电子熵）与晶格软化/硬化重组（振动熵）双轮驱动的结果。这一洞察对设计高性能超导材料、巨磁阻材料以及高性能核燃料具有重大的指导意义。