逐图读 O’Hern 等人 2003 年的论文，找出堵塞点

先看协议：构型从哪里来

第 II A 节和式 (1) 定义了模型：两个颗粒只有重叠时才相互排斥，势能为 $\varepsilon(1-r_{ij}/\sigma_{ij})^{\alpha}/\alpha$。二维用大小比 1.4、各占一半的两种圆盘，防止结晶；粒子数从 4 到 4096。

构型从无穷高温出发：圆心在周期性盒子里完全随机地放，再用共轭梯度一路往能量低处走。停止条件有两个：每个粒子的能量低于 $10^{-16}$，或者相邻两次迭代的能量差小于 $10^{-15}$。第一个条件就是“未堵塞”的判据，对应压强低于 $10^{-10}$。

读到这里要记住：论文里的“堵塞”是这个协议下的结论。换一种生成方法，比如先慢慢降温，得到的阈值会不同。论文后面所有的分布和外推，都是对这种随机起点的统计。

来源：C. S. O’Hern

一个构型的阈值怎样测（第 II B 节）

对每个零温构型：压强为零就把所有颗粒按同一比例放大，压强不为零就缩小，每步之后重新求极小，直到压强刚好变为非零。步长在 $10^{-6}$ 到 $10^{-4}$ 之间，越接近阈值越小，并且逐次减小步长，确认最后落在同一个构型上。

这一节的关键结论是：体积模量和剪切模量在同一个堆积分数上变为零。也就是说，“开始抗压”和“开始抗剪”是同一件事，阈值对每个构型都有确切定义。作者还检查了动力学矩阵：阈值之上至少 1e-6 处，除了松动颗粒和整体平移，没有零频模式。

图 3 有一处要小心。我们读的 arXiv v1 中，图注说参考线斜率是 2.0 和 2.5；而表 I 的压强指数是 1.01 和 1.50，第四章从 $p\propto\Delta\phi^{\alpha-1}$ 推出的也是这两个数。2.0 和 2.5 恰好是能量的指数 $\alpha$。我们按表 I 检验压强；正式发表版的图注我们没有看过。

来源：C. S. O’Hern

从堵住比例到阈值分布（图 5 到图 8）

图 5 画的是在给定堆积分数下，随机起点中已经堵住的比例 $f_j(\phi)$。它是一条从 0 升到 1 的台阶曲线。图 6 是它对堆积分数的导数，也就是阈值的分布 $P_j(\phi_c)$。作者把三种接触指数的分布叠在一起，在误差内完全重合：阈值只由几何决定，与接触力的软硬无关。

图 7 是分布的半高全宽随粒子数的变化。读图时注意两段：粒子很少时宽度先随粒子数变大，越过最大值后才按式 (3) 下降。式 (3) 的 $w_0=0.16\pm0.04$、$\Omega=0.55\pm0.03$ 只适用于后一段。

图 8 处理的是分布峰的位置 $\phi_0$。它随系统长度 $L=N^{1/d}$ 按式 (4) 趋近极限值 $\phi^{*}$，拟合给出 $\nu=0.71\pm0.08$、$\delta_0=0.12\pm0.03$。可见 $\phi^{*}$ 是外推出来的，不是哪一次模拟直接测到的。

来源：C. S. O’Hern

接触数与表 I：哪些数依赖假设

第 II E 节用式 (6) $Z-Z_c=Z_0(\phi-\phi_c)^{\zeta}$ 拟合图 9，正文给出 $\zeta=0.50\pm0.03$。表 I 里的 $Z_c=3.98$ 是另一回事：作者先去掉松动颗粒（约 5%），再假定 $\zeta=0.5$，才得到这个“精确”的值。$Z_c$ 与 4 相符，并不能反过来独立证明指数是 1/2。

Goodrich、Liu 与 Nagel 2012 年把计数做得更细：周期性盒子里 N 个圆盘，要有正的体积模量至少需要 $2N-1$ 个接触，并且当 $p^{1/2}N$ 接近 1 时，接触数的平方根律会失效。读 2003 年的表 I 时，要把这条有限尺寸限制带上。

我们的复算：在阈值处的承载网络上，超出最少接触数的接触在 32 和 64 个盘时都是 0，在 128 和 256 个盘时最多 1 个；松动盘平均占 4.0% 到 5.0%，与论文的约 5% 一致。接触数指数在 128 个盘时得到 0.577 ± 0.009，没有复现 0.49。

来源：C. S. O’Hern · C. P. Goodrich

完整例题：用式 (3) 和式 (4) 检查我们的阈值统计

问题：我们在周期性盒子里对 32、64、128、256 个圆盘各取 24 个随机起点，测得阈值分布的半高全宽为 0.0259、0.0148、0.0129、0.00654，平均阈值为 0.827、0.835、0.838、0.839。它们与式 (3)、式 (4) 相容吗？无限大系统的阈值大约是多少？

第一步，宽度。式 (3) 的中心值给出 $0.16\times N^{-0.55}$：0.0238、0.0162、0.0111、0.00758。取 $w_0$ 和 $\Omega$ 不确定度的两端，256 个盘允许的范围是 0.0048 到 0.0112，我们的 0.00654 在里面；其余三个尺寸也都在各自的范围里。注意我们是用标准差乘 $2\sqrt{2\ln2}$ 估计全宽，假设了分布是高斯形。

第二步，外推。二维时 $L=\sqrt{N}$，256 个盘 $L=16$，式 (4) 的中心值给出 $0.12\times16^{-1/0.71}=0.0024$。平均阈值随尺寸上升，说明峰从下方趋近极限，所以 $\phi^{*}\approx0.839+0.0024=0.841$。按 $\delta_0$ 和 $\nu$ 的不确定度，这个修正在 0.0011 到 0.0045 之间。

结论：宽度与式 (3) 相容；外推值约 0.841，落在图 5(a) 中大系统 0.84 附近。这个估计有两处近似：式 (4) 讲分布的峰，我们用 24 个样本的平均值代替；256 个盘本身也只是中等大小。

来源：C. S. O’Hern

读这类论文的清单

一、先找协议：构型怎样生成、用什么判断堵住、收敛阈值多大。二、分清直接测量和拟合：J 点的位置、接触数的极限都来自外推或假设。三、找适用范围：式 (3) 只对大系统成立；按图 3 和图 9 的图注，二维的标度数据来自 1024 个圆盘。四、对照图注和正文：像图 3 这样的出入，以表和推导为准，并记下自己没有核实的部分。

五、复算时写明什么算相符。我们要求宽度落在式 (3) 不确定度给出的范围里，接触数指数落在表 I 的 $\pm0.04$ 之内。前者通过；后者 0.577 对 0.49，没通过。两个结果都记录在案。

来源：C. S. O’Hern · C. P. Goodrich