爆发性逾渗模型的Shortest-path指数和Backbone指数

，

(合肥工业大学电子科学与应用物理学院，安徽合肥230009)

几何模型在统计物理和复杂系统的研究中扮演着重要的角色，为相变和临界现象的理论研究提供了载体[1]，著名的例子包括Random Cluster模型、元胞自动机模型、逾渗模型等[2-8]。逾渗概念的提出距今已有约60年，最初用于模拟自然界中的森林火灾、流行病的传播、多孔介质中石油的流动等现象[6]。 随后的研究显示该模型存在非平凡的相变，因而引起统计物理学家的广泛关注[7-8]。 普通逾渗中，人们以概率p往一个格子的连边(edge)上放键(bond)，以1-p的概率不放。导致无穷大集团在无穷大格子上出现的最小p值记为pc，称为逾渗阈值。爆发性逾渗采用了不同于普通逾渗的放键规则[9]。每打算往格子上放一根键，事先随机选择两条edge作为候选放键位置，但最终只有一条edge的位置被放键，另一条edge的位置空置。 “二选一”时，常采用“乘积规则”，即分别计算两条候选edge的端点所属的两个集团大小的乘积，其中使得乘积值较小的edge将被放键。 相比于普通逾渗，该规则使得巨大集团的形成推迟到来，而中等大小的集团的数目增多。 一旦接近相变点，集团大小可能会因极少数目键的加入而呈现爆发式增长，因而得名爆发性逾渗。爆发性逾渗最初被怀疑经历一级相变，随后的研究表明其仍为连续相变，但拥有极其微小的序参量临界指数β[10]。爆发性逾渗模型的一个典型应用实例是Clusella等[11]曾利用它来探讨计算机网络的鲁棒性。

爆发性逾渗与普通逾渗分属不同的普适类，各自拥有一套独立的临界指数。借助于库伦气、共形场以及SLE等理论[12-13]，二维普通逾渗的临界指数大多能够解析地给出，但shortest-path指数dmin和backbone指数dB是例外，它们目前只有数值结果[14-15]。这两个指数的定义可简单地表达为：

lS～rdmin,lB～rdB,

(1)

式中r代表当系统处在pc处时同属一个无穷大集团的两个点a、b之间的欧几里得距离；lS代表a、b之间的最短路径，在编程语言中，它指从a点出发通过键路径到达b所需的最少的MC步数。在某些文献中它也被称作化学距离(chemical distance)[16]。为了解释lB，我们设想键为导体棒，edge本身为绝缘体，现施电势差于a、b之间，则所有的载流导体棒都称为backbone键，由backbone键组成的路径上所有格点的数目即为lB。由于观测量lS和lB较为复杂，在实际的蒙特卡洛模拟中，人们往往借助一些其他更易于测量的物理量来间接获得这两个指数，且相应的power-law关系中的自变量不再是r，而是线性系统尺寸L。例如，为获得dB，文献[14]的观测量不是lB，而是backbone关联函数。近年来，在普通逾渗问题的研究中，人们找到了两个易于对dmin、dB实现精确测量的物理量[17-19]，分别是完全(complete)构型上的图形距离S和bridge-free构型上的最大集团大小Cbf，它们在pc处满足的标度关系是

S～Ldmin,Cbf～LdB,

(2)

这两个观测量的定义将在第1节中作简要介绍。

对于爆发性逾渗，截至目前，shortest-path指数和backbone指数既无解析解，也无数值解。在本工作中，我们假定公式(2)对爆发性逾渗依然有效，并利用蒙特卡洛模拟研究了周期性边界条件下正方晶格上满足乘积规则的爆发性逾渗，通过测量S和Cbf并进行数据分析，获得了关于这两个临界指数的首个估计。

1 蒙特卡洛模拟

本文在正方格子上模拟了遵守乘积规则的爆发性逾渗模型，其逾渗阈值为pc=0.526 565(5)[20]，实际的模拟点为p=0.526 563 6，其位于文献值的一个误差棒之内。这样选择的原因在于，我们的另一个尚未发表的工作显示，该模型的临界点可更精确地表达为pc=0.526 563 6(10)。在模拟中，正方晶格的线性尺寸涵盖从L=4到L=8 192的多个尺寸，例如，L=4, 6, 8, 12, …, 6 144, 8 192。每个尺寸的样本数不尽相同，总的来说大系统尺寸的样本数少，而小系统尺寸的样本数较多，例如，L≤512的各系统，样本数均达到2千万左右，而L≥6 144的大系统，样本数仅10万左右。

本工作中的两个观测量分别是在完全构型和bridge-free构型上取样而得，下面介绍这两种构型产生的过程和相应观测量的取样方法。

1.1 产生完全构型

为了产生爆发性键逾渗的完全构型，本文采用如下算法：

Step1: 从一个空的格子开始，每个格点i自成一个集团，集团编号与格点自身的编号相同；

Step2:从所有未被键占据的edge中随机选择两条候选edge(分别记为e1和e2)，计算e1的两个端点所属的两个集团大小的乘积；对e2作类似的操作；

Step3:使得乘积较小的edge位置将被放键，另一根edge处将被空置；若乘积相等，则从这两根edge中随机地选择一根edge进行占据；

Step4:将新加入键的两个端点所属的两个集团进行合并，并将合并后的集团编号进行更新；

Step5:重复Step2～4，直到格子上键总数达到预设值(等于pc×edge总数并取整)。

值得一提的是，在Step4中，我们事实上仅更新参与合并的两个集团中较小的那个集团的编号，使其与较大的那个集团编号一样。

键构型产生后，采用Swendsen-Wang算法统计各个集团的大小[21]。在具体实施时，采用breadth-first搜索方案。在搜索时，从某个集团的具有最小编号的格点(称为第0代seed site)出发，搜索它的所有最近邻邻居，其中位于集团上的那些邻居构成新一代seed sites(记为第1代)；接下来搜索第1代seed sites的所有最近邻邻居，那些位于集团上且此前未被搜索过的邻居点构成第2代seed sites；依此方案可以继续搜索第3代，第4代，…，第m代。 其中m是该集团中最后一个被搜索的点所属的seed sites代数。很显然，构型上的每一个集团都存在一个确定的m，值最大的记为m0。S=〈m0〉是m0的系综平均值，方便起见，我们称其为图形距离。

1.2 产生bridge-free构型

将完全构型上的所有占据键根据其连通性特征分为两类：桥键和非桥键。所谓桥键是指若被删除，其两个端点将分属不同的集团；若被删除后仍属于同一集团，则该键为非桥键。若采用一定的算法删除掉构型上所有的桥键，就得到bridge-free构型。为了获得backbone指数的估计，我们测量bridge-free构型上最大集团的大小，并取其系综平均值Cbf。

b为树枝键；j为枢纽键；未作标记为非桥键。 图1 3种构型及3类键的例子Fig.1 An example for illustrating three types configurations

在实际的模拟中，为了更方便地删除桥键，采用了类似文献[18-19]的做法，将所有桥键进一步分类为树枝键和枢纽键。树枝键是指若其被删除，其两个端点中至少有一个属于单点集团；一根树枝键的删除有可能导致新的树枝键产生，因而在蒙特卡洛模拟中删除树枝键的过程是迭代式进行的。排除掉所有树枝键后的构型，称为leaf-free构型。接下来在leaf-free构型上采用“回溯”(backtracking)法识别非桥键并进行标记。回溯法识别非桥键的基本策略是，从集团上的一个点出发采用breadth-first方案搜索最近邻邻居，并将每一代seed sites与下一代邻居点link起来。比如，用于link的数组father(j)=i表示j号格点是作为上一代seed sitei的下一代邻居被搜索到并加入link列表的。 一旦某一代seed sitej′搜索到位于同一集团且已被访问过的点i′，则意味着一个由键组成的封闭的环路被找到，环路径上的每一根键必为非桥键。为了标记这个环上的所有键为非桥，需要从i′,j′出发沿link列表反向查找各自的father，这个过程迭代地进行，直到找到共同的father为止。其他的环状路径也是采用这种方式进行查找和标记。当整个识别过程结束时，leaf-free构型上所有被标记为非桥的键组成的子构型即为bridge-free构型。 为方便直观理解，本文提供了3种构型及3类键的一个例子(图1)。

2 数据分析

理论上，式(2)只在热力学极限(系统为无穷大)条件下严格成立。考虑到蒙特卡洛方法只能模拟有限大小的系统，因此必须考虑有限尺度效应。式(3)对S和Cbf进行拟合：

F(L)=LdO(a+bLy1)，

(3)

式中F代表观测量S或Cbf，dO代表临界指数dmin或dB，y1代表有限尺度修正的leading项修正指数。拟合的流程按照least-square标准进行。一个可以被接受的拟合应保证χ2/dDF在O(1)数量级，其中χ2代表数据与拟合公式之间的残差，dDF是自由度，定义为数据点的数目与拟合参数个数之差。实际操作中，一般将χ2/dDF≈1的拟合视为合理的。拟合结果呈现在表 1中，Lmin表示被用于拟合的数据来自最小尺寸为Lmin的那些系统。

2.1 Shortest-path指数dmin

在对S的拟合中，当y1作为自由参数时，发现y1≈-0.35。为了减少一个拟合参数，将y1固定在y1=-0.35进行拟合，拟合结果见表1。为了判断y1的轻微变化对临界指数dmin的估计值的影响，将y1固定在其他值进行了拟合(未显示于表格中)。 此外，还尝试过在式(3)中添加subleading修正项cL-2进行拟合，但当添上该项时，修正幅值b和c以及leading修正指数y1均不能被MC数据所分辨(不确定度与中心值相当，甚至大于中心值)。综合各种各样的拟合，我们给出关于shortest-path指数的最终估计结果为dmin=1.189(3)，为显示关于误差估计的可靠性，取dmin等于中心值、中心值加上3个误差棒、中心值减去3个误差棒分别画了S/Ldmin关于L-0.35的数据点线图(图2)。理论上，dmin等于中心值时数据点应呈现出一条近似的直线，而dmin等于另外两个值时，在L很大的区域应呈现出向上或向下的弯曲。图2清楚地展示出了这种预期中的趋势，表明关于dmin的估计是可靠的。

表1 观测量S的拟合结果

图2 S/Ldmin关于L-0.35的点线Fig.2 Points and curve of S/Ldmin versus L-0.35

Lminχ2/dDFdBaby12563847.4/74.5/61.545 0(5)1.544 0(8)1.002(4)1.009(6)-3.0(3)-3.9(6)-1-13845123.9/63.4/51.546 1(5)1.545 7(8)0.992(4)0.995(6)-696.0(100)-885.0(283)-2-27683.5/51.546 6(5)0.987(4)--

2.2 Backbone指数dB

对Cbf首先实施让y1作为自由参数的拟合，但y1和b均不能被MC数据所分辨。为此，将y1固定在常见的修正指数y1=-1和-2处进行拟合，以评估修正指数y1的变化对临界指数dB的影响。我们还尝试了抛弃所有修正项，直接使用O(L)=aLdO进行拟合，拟合结果见表2。综合这些拟合以及其他未显示于表格中的拟合，得到关于dB的最终估计为dB=1.546(5)。为了展示结果的可靠性，取3个不同的dB值(中心值以及加、减3个误差棒)去画Cbf/LdB关于L-1的数据曲线。显然，当系统无穷大(L-1→0)且dB为中心值时，Cbf/LdB应等于常数a(≈1)，数据点在图中应展示出一条近似的水平直线，而另外2个dB值则导致曲线在L-1→0的区域呈现出向上或向下的弯曲。这种预期中的现象在图3中被清楚地观察到，表明关于dB的估计也是可靠的。

图3 Cbf/LdB随L-1变化的点线Fig.3 Points and curve of Cbf/LdB versus L-1

3 讨论和结论

通过蒙特卡洛模拟和细致的数据分析，获得了周期性边界条件下正方晶格上满足乘积规则的爆发性逾渗模型的Shortest-path指数和Backbone指数分别为dmin=1.189(3)和dB=1.546(5)。这两个指数(dmin,dB)在爆发性逾渗类系统中从未被解析地确定或数值地估计过， 因此本文的工作提供了关于它们的第一个估计。本工作预期达到的目的包含两个方面：一是为本领域同行提供可靠的科学数据，抛砖引玉，期待能为人们将来解析地研究爆发性逾渗提供检验基础；二是对测量dmin、dB的新方法的一次重要应用，在今后的工作中，该方法还将应用到更多的模型，如Ising模型和Potts模型。

对比后不难发现，爆发性逾渗的Shortest-path指数和Backbone指数均不同于普通逾渗，后者相应的结果分别为dmin,普=1.130 77(2)[15]和dB，普=1.643 4(2)[14]。这表明两个模型拥有不同的集团结构特征，尤其是在集团的紧致性(compactness)方面。尽管爆发性逾渗拥有更大的pc值(或占据键浓度)，但其最大集团内部格点彼此之间的连通度(或集团的紧致性)略低于普通逾渗的最大集团。此外，本文的结果也从临界指数dmin、dB的角度再次印证爆发性逾渗模型与普通逾渗模型的确分属于不同的普适类。由本工作出发还引申出一系列有待回答的问题。如前所述，本文研究的爆发性逾渗模型采用了乘积规则和周期性边界条件。 一个自然的问题是：由该模型获得的结果是否具有普遍性？例如，当乘积规则改为求和规则，或者周期性边界条件改为开放边界条件时，dmin、dB的值是否会改变？探究这些问题，将有助于加深人们对具有非平凡规则的统计模型相变普适类性质的理解。