尋之朋 郝大鵬

(中國礦業(yè)大學(xué)材料與物理學(xué)院,徐州 221116)

1 引言

滲流模型蘊(yùn)含著豐富的物理,是統(tǒng)計物理學(xué)領(lǐng)域的基礎(chǔ)模型之一[1-6].一方面,研究滲流模型對理解自然界中的相變和臨界現(xiàn)象、分形、標(biāo)度理論等均有著重要的科學(xué)價值.另一方面,滲流具有廣泛的應(yīng)用,如流體流經(jīng)多孔物質(zhì)、森林火災(zāi)、流行病的傳播(如最近爆發(fā)的新冠肺炎)等.因此,對滲流的深入研究具有重要的理論價值和現(xiàn)實(shí)意義.

建立各種格點(diǎn)滲流模型是研究滲流常用的方法之一.在指定格點(diǎn)模型上,每一個格點(diǎn)(或鍵)是獨(dú)立的并以一定的概率p被占據(jù)(不被占據(jù)的概率為 1-p),被占據(jù)的格點(diǎn)(或鍵)可以形成毗連的團(tuán)簇.當(dāng)占據(jù)概率p從 0 逐漸增大到某一臨界值時,格點(diǎn)模型上將開始出現(xiàn)大到能夠貫穿整個系統(tǒng)的團(tuán)簇,并且這時的系統(tǒng)通常表現(xiàn)出連續(xù)相變的特性,稱這樣的系統(tǒng)為滲流,并且稱此時的臨界概率為滲流閾值,用pc表示.滲流閾值pc是滲流的核心參數(shù)之一,只有準(zhǔn)確地確定滲流閾值,才能更好地研究其臨界行為以及確定臨界指數(shù).

無論從解析角度還是采用蒙特卡羅模擬等數(shù)值方法,許多格點(diǎn)模型被廣泛研究以獲得滲流閾值pc以及相應(yīng)的臨界指數(shù).在眾多格點(diǎn)模型中,含復(fù)雜近鄰的格點(diǎn)滲流模型由于多方面的應(yīng)用而成為該領(lǐng)域的熱點(diǎn)研究課題之一.例如,該類模型點(diǎn)滲流可以鏡像到格點(diǎn)模型上擴(kuò)展形狀(如圓形、球形)的吸附滲流過程[7,8].鍵滲流不僅與小世界網(wǎng)絡(luò)有著密切的聯(lián)系[9],還可以從滲流的角度去研究分析流行病的傳播過程[10,11].其次對這類滲流模型的深入研究,既能豐富相關(guān)的滲流理論,又對研究滲流閾值與格點(diǎn)模型的結(jié)構(gòu),尤其是格點(diǎn)配位數(shù)之間的關(guān)系具有重要的指導(dǎo)作用.此外,這類格點(diǎn)滲流模型的結(jié)構(gòu)介于離散滲流與連續(xù)滲流之間,對其進(jìn)行深入探討也有助于揭示離散滲流和連續(xù)滲流之間轉(zhuǎn)變的物理機(jī)制.

對含復(fù)雜近鄰的格點(diǎn)滲流模型的研究可以追溯到由Dalton,Domb 和Sykes[12-14]提出的“等效近鄰模型”.之后,研究人員在該方向開展了廣泛的研究,如菱形結(jié)構(gòu)緊密區(qū)域的長程點(diǎn)滲流[15],考慮最近鄰和次近鄰的體心立方格子的點(diǎn)滲流和鍵滲流[16],含第四近鄰格點(diǎn)的面心立方格子的點(diǎn)滲流[17]等.d’Iribarne 等[18-20]從圖論的角度系統(tǒng)計算了長程關(guān)聯(lián)(最大計算到了第10 近鄰格點(diǎn))的全部11 種阿基米德格子的點(diǎn)滲流.Malarz 和Galam[21]引入“復(fù)雜近鄰”的概念,即不同近鄰格點(diǎn)的各種組合而并不局限于緊密近鄰格點(diǎn),隨后該課題組對二維正方格子、三維簡單立方格子、四維超立方格子的點(diǎn)滲流開展了一系列的研究工作[22-26].最近,基于多種有效算法的發(fā)展,含復(fù)雜近鄰的格點(diǎn)模型的鍵滲流也得到了廣泛的研究[27-30].

滲流閾值pc與格點(diǎn)模型結(jié)構(gòu)(如配位數(shù)z)之間的依賴關(guān)系一直是探討的焦點(diǎn)之一.例如文獻(xiàn)[23]指出一些三維格點(diǎn)模型的點(diǎn)滲流閾值與格點(diǎn)配位數(shù)之間滿足冪律關(guān)系pc～z-a,指數(shù)a=0.790(26) .Galam 和Mauger[31],van der Marck[32],以及其他研究人員也分析研究了許多其他體系,得出類似的冪律關(guān)系[29,30].Domb[12]指出在相同的格點(diǎn)單元形狀下,含擴(kuò)展近鄰的格點(diǎn)模型點(diǎn)滲流在配位數(shù)z較大時的漸近行為與連續(xù)滲流閾值ηc有關(guān);該論點(diǎn)得到了其他學(xué)者的贊同和推進(jìn)[7,8,20],最近更深入的研究[33]揭示出該漸近關(guān)系為pc～2dηc/z.對于鍵滲流,該漸近關(guān)系趨于Bethe 格子的結(jié)果pc～1/(z-1);以及當(dāng)z為有限值時滿足有限尺寸修正zpc-1～a1z-x,其中在二維及三維情況下x=(d-1)/d[34,35].這些理論分析結(jié)果均得到了大量數(shù)值模擬工作的驗證[27,28,33,35].

但是,以上漸近行為及有限尺寸修正行為僅對含緊密近鄰的格點(diǎn)滲流模型是有效的,而對于含復(fù)雜(非緊密)近鄰格點(diǎn)滲流模型的分析卻遇到了困難.其主要原因是這類格點(diǎn)滲流模型表現(xiàn)出所謂的“簡并”現(xiàn)象,即一個配位數(shù)z對應(yīng)于多個格點(diǎn)模型結(jié)構(gòu),如圖1 所示的含復(fù)雜近鄰格點(diǎn)的二維正方格子(圖中“?”代表中心格點(diǎn),①代表最近鄰或第1近鄰格點(diǎn),②代表次近鄰或第2 近鄰格點(diǎn),······),“第1 近鄰+第2 近鄰”與“第1 近鄰+第3 近鄰”所對應(yīng)的配位數(shù)是相同的,均為8;但是這兩種結(jié)構(gòu)的滲流閾值卻是不相同的.因此,對于這類格點(diǎn)滲流模型,僅僅考慮配位數(shù)一個參量是不夠的,如何綜合考慮格點(diǎn)模型的拓?fù)浣Y(jié)構(gòu)來消除“簡并”是解決該問題的關(guān)鍵所在.

圖1 正方格子中心點(diǎn)(“ ? ”)的近鄰格點(diǎn),到第 5 近鄰Fig.1.Neighbors of a central site (“ ? ”) on a square lattice,up to 5 th nearest neighbors.

為了更深入地探討滲流閾值與格點(diǎn)模型的幾何或拓?fù)浣Y(jié)構(gòu)之間的關(guān)聯(lián),本文基于圖論[36]、有機(jī)化學(xué)[37,38],以及最近Malarz[39]關(guān)于點(diǎn)滲流中模型參數(shù)識別的基本思想,引入新的標(biāo)量參數(shù)ξ=其中zi和ri分別為第i近鄰格點(diǎn)的配位數(shù)及其到中心格點(diǎn)的距離,綜合地考慮格點(diǎn)滲流模型的拓?fù)浣Y(jié)構(gòu)來消除這種“簡并”,并將該標(biāo)量參數(shù)應(yīng)用于研究分析含復(fù)雜近鄰的格點(diǎn)模型的鍵滲流.同時,基于一種高效的單團(tuán)簇生長算法對含復(fù)雜近鄰的二維正方格子鍵滲流進(jìn)行了大量的蒙特卡羅模擬研究,得到了這些格點(diǎn)滲流模型高精度的鍵滲流閾值.所得到的滲流閾值與新的標(biāo)量參數(shù)很好地滿足冪律關(guān)系pc∝ξ-γ,數(shù)值擬合得指數(shù)γ≈1 .

2 理論基礎(chǔ)

團(tuán)簇的統(tǒng)計分布是研究滲流的一個中心物理量,通過該分布的特性即可確定滲流閾值及相應(yīng)的臨界指數(shù).首先定義團(tuán)簇的尺寸為一個團(tuán)簇中所含格點(diǎn)(或鍵)的個數(shù),用s表示.則團(tuán)簇的分布通常用ns(p) 表示,它是指在占據(jù)概率p下,尺寸為s個格點(diǎn)(或鍵)的團(tuán)簇的數(shù)目.在實(shí)際使用時,將對團(tuán)簇的分布進(jìn)行歸一化處理.在滲流閾值pc處,ns(pc)滿足關(guān)系式

其中,τ為Fisher 指數(shù),Ω為冪律關(guān)系的一階修正指數(shù).τ和Ω均被認(rèn)為具有普適性,只與系統(tǒng)的維度有關(guān).A0和B0為與具體的模型有關(guān)的兩個參數(shù),是非普適的.由(1)式可得出模型上每個格點(diǎn)(或鍵)屬于尺寸大于或等于s的團(tuán)簇的概率為

其中,A1=A0/(τ -2),B1=(τ -2)B0/(τ+Ω-2) .(2)式兩邊同時乘以sτ-2得

可以看出,當(dāng)s較大時,sτ-2P≥s隨s-Ω的變化將呈線性關(guān)系.利用該線性關(guān)系可以確定滲流閾值,因為當(dāng)p/=pc時,關(guān)系是非線性的.

當(dāng)p/=pc時,需要引入標(biāo)度函數(shù).在標(biāo)度極限s →∞和p→pc下,(p-pc)sσ趨于常數(shù),P≥s滿足

式中σ為另一普適指數(shù).對標(biāo)度函數(shù)f(x) 作泰勒展開

其中,C2=B2f′(0) .這里假設(shè)f(0)=1,將(5)式代入(4)式可得

其中D2=A2C2.(6)式表明,當(dāng)s較大時,sτ-2P≥s在p→pc時將趨于常數(shù),而在p遠(yuǎn)離pc時將偏離該常數(shù).這提供了確定滲流閾值的另一個途徑.

3 數(shù)值模擬結(jié)果與分析

單團(tuán)簇生長算法的基本規(guī)則[29,30]如下:在格點(diǎn)模型上任取一個格點(diǎn)作為種子,以該種子為起點(diǎn)開始生長單獨(dú)的團(tuán)簇.在周期邊界條件下,格點(diǎn)模型上任意一個格點(diǎn)都可以作為種子,本文在數(shù)值模擬中選取為坐標(biāo)原點(diǎn).從種子格點(diǎn)開始,團(tuán)簇以一定的概率p占據(jù)(不被占據(jù)的概率為 1-p)與之連接的鍵從而向近鄰格點(diǎn)生長.在系統(tǒng)尺寸為L×L=16384 × 16384 的正方格子上,針對每一個格點(diǎn)滲流模型,獨(dú)立地生長不少于 5×108個系綜(即獨(dú)立的團(tuán)簇).從統(tǒng)計上來講,這些團(tuán)簇具有不同的尺寸,尺寸較小的團(tuán)簇將很快停止生長,而有的團(tuán)簇則可以生長到很大的尺寸.因此在實(shí)際的數(shù)值模擬中,需要設(shè)定團(tuán)簇生長的上臨界截斷閾值尺寸來避免出現(xiàn)邊界環(huán)繞,同時控制程序運(yùn)行時間,當(dāng)團(tuán)簇的生長達(dá)到預(yù)設(shè)的上臨界截斷閾值尺寸時停止生長.然后,將尺寸落在區(qū)間 [ 2n,2n+1-1] (n=0,1,2,···)內(nèi)的團(tuán)簇的數(shù)量記錄到結(jié)果向量的第n個值中.對于達(dá)到上臨界截斷閾值尺寸仍在生長的團(tuán)簇,它將被統(tǒng)計到結(jié)果向量的最后一個值中.對于本文所模擬的所有格點(diǎn)模型,上臨界截斷閾值均設(shè)置為 216,這也就意味著結(jié)果向量的長度為17.采用單團(tuán)簇生長算法在數(shù)值模擬中不需要存儲記錄大量的中間迭代數(shù)據(jù),僅需要對結(jié)果向量進(jìn)行處理即可,這是該算法優(yōu)于其他方法的一個方面.此外,新的單團(tuán)簇生長算法[29,30]在原有的基礎(chǔ)上采用了有效的技術(shù)手段對算法進(jìn)行改進(jìn),成功地避免了每次單個團(tuán)簇生長之后將整個格點(diǎn)模型信息清除.優(yōu)化的基本思想如下:將每個格點(diǎn)的初始值設(shè)置為零,對于第n個團(tuán)簇(對應(yīng)于統(tǒng)計系綜次數(shù)),晶格模型上任何值小于n的格點(diǎn)都被視為未被占據(jù).當(dāng)一個格點(diǎn)在新團(tuán)簇生長的過程中被占據(jù)時,其值賦為n.這個過程節(jié)省了大量時間,因此我們可以模擬非常大的晶格尺寸,并且不必在每個團(tuán)簇(實(shí)際上許多團(tuán)簇都很小)生長之后清除整個晶格的信息,使計算效率得到了顯著的提高.在文獻(xiàn)[29]中,我們給出了該優(yōu)化算法的程序偽代碼.

對于普適的標(biāo)度指數(shù),目前已得到二維情況下的解析 解τ=187/91,σ=36/91和Ω=72/91,本文在數(shù)值分析中直接使用這些普適參數(shù)值.于是根據(jù)蒙特卡羅模擬數(shù)據(jù),相應(yīng)的物理量,如sτ-2P≥s,便可易于計算得到.

本文使用 S Q-a,b,···來表示含第a近鄰,第b近鄰,······的二維正方格子.基于理論公式(6),圖2為在占據(jù)概率分別為 0.250365 ,0.250367 ,0.250368,0.250369 及 0.250371 時 S Q-1,2 格子 鍵滲流sτ-2P≥s隨sσ的變化曲線.可以看出:當(dāng)團(tuán)簇尺寸s較小時,由于有限尺寸效應(yīng),sτ-2P≥s隨sσ急劇下降;當(dāng)s較大時,sτ-2P≥s隨sσ的變化逐漸表現(xiàn)出線性行為,并且隨著占據(jù)概率p趨于滲流閾值pc,曲線逐漸趨于水平,即趨于常數(shù).根據(jù)圖2 中線性部分的性質(zhì),對每一個占據(jù)概率下曲線的線性部分計算斜率,滲流閾值的中心值可以通過作圖

來進(jìn)行確定,如圖2 中的插圖所示.通過插圖中的線性部分與橫軸的交點(diǎn),可以確定滲流閾值的中心值為pc=0.2503683 .

圖2 SQ-1,2 格子鍵滲流在不同占據(jù)概率p 下sτ-2P≥s隨 s σ 的變化曲線,其中 τ =187/91 ,σ =36/91 .插圖表示主圖中所示結(jié)果線性部分的斜率隨占據(jù)概率p 的變化關(guān)系Fig.2.Plot of s τ-2P≥s vs.s σ with τ =187/91 and σ=36/91 for the bond percolation of the S Q-1,2 lattice under different values of p.The inset indicates the slope of the linear portions of the curves shown in the main figure as a function of p.

另一方面,基于理論公式(3),在占據(jù)概率分別為p=0.250365 ,0.250367 ,0.250368 ,0.250369,以 及 0.250371 時 S Q-1,2 格子鍵滲流sτ-2P≥s隨s-Ω的變化曲線如圖3 所示.由圖3 可以很明顯地看出,當(dāng)占據(jù)概率p偏離滲流閾值pc時,曲線表現(xiàn)出明顯的偏離線性行為;而當(dāng)占據(jù)概率p逐漸趨于滲流閾值pc時,曲線表現(xiàn)出越來越好的線性行為.于是,通過圖3 可以確定鍵滲流閾值的范圍為0.250368＜pc＜0.250369.

圖3 SQ-1,2 格子鍵滲流在不同占據(jù)概率p 下sτ-2P≥s隨 s -Ω 的變化曲線,其中 τ =187/91 ,Ω=72/91Fig.3.Plot of s τ-2P≥s vs.s -Ω with τ =187/91 and Ω=72/91 for the bond percolation of the S Q-1,2 lattice under different values of p.

綜上兩種方法,可以確定 S Q-1,2 格子鍵滲流閾值為pc=0.2503683(7),其中括號里面的數(shù)字表示閾值末位的誤差.對于本文所模擬的其他格點(diǎn)模型,其分析過程與上述類似;因此不再贅述每一個格點(diǎn)模型具體的結(jié)果圖形及分析過程,而是直接給出鍵滲流閾值的數(shù)值處理結(jié)果,并將結(jié)果總結(jié)在表1 中.

表1 含復(fù)雜近鄰格點(diǎn)的二維正方格子的鍵滲流閾值Table 1.Bond percolation thresholds of square lattice with complex neighborhoods.

表1 第2 列給出了每一個格點(diǎn)模型的配位數(shù)z.可以看出存在明顯的“簡并”現(xiàn)象,即許多不同的格點(diǎn)模型具有相同的配位數(shù)z,但鍵滲流閾值卻不相同.圖4 給出了這些格點(diǎn)模型鍵滲流閾值pc隨格點(diǎn)配位數(shù)z變化的對數(shù)-對數(shù)曲線,也可以明顯地看出“簡并”行為.這些結(jié)果說明,對于含復(fù)雜(非緊密)近鄰的格點(diǎn)滲流模型,滲流閾值不僅僅取決于格點(diǎn)配位數(shù),而需要對格點(diǎn)模型的其他幾何參數(shù)進(jìn)行更深入的分析.

圖4 表1 中所示格點(diǎn)模型 pc 隨z 變化的對數(shù)-對數(shù)曲線Fig.4.The log-log plot of pc vs.z for the lattices shown in Table 1.

以上分析表明需要采取有效的技術(shù)手段,引入新的參數(shù)來消除“簡并”.最近,Malarz[39]綜合考慮第i近鄰格點(diǎn)的配位數(shù)zi以及其到中心格點(diǎn)距離的平方引入了新的標(biāo)量參數(shù)

來有效消除“簡并”,并將其應(yīng)用于分析含復(fù)雜近鄰格點(diǎn)的二維三角格子的點(diǎn)滲流,取得了成功.實(shí)際上,這種消除簡并的方式在圖論中分析樹結(jié)構(gòu)的拓?fù)洳蛔冃訹36],有機(jī)化學(xué)中分析分子拓?fù)渲笖?shù)[37,38]等方面也有成功的應(yīng)用.本文嘗試采用這種方法來分析含復(fù)雜近鄰二維正方格子的鍵滲流.首先,根據(jù)(8)式計算得到每一個格點(diǎn)模型對應(yīng)的標(biāo)量指數(shù)ξ,如表1 中第3 列所示.對應(yīng)于不同近鄰格點(diǎn)的zi和ri2參見表2.從計算結(jié)果可以看出,引入新的參數(shù)后,在很大程度上消除了“簡并”現(xiàn)象.然后,作出滲流閾值pc隨參數(shù)ξ變化的對數(shù)-對數(shù)曲線,如圖5 所示.結(jié)果顯示,滲流閾值pc隨ξ的變化呈現(xiàn)出很好的冪律關(guān)系pc∝ξ-γ,并且線性擬合得到冪律指數(shù)γ≈1 .可見,(8)式消除簡并的方法也能夠成功地應(yīng)用于含復(fù)雜近鄰格點(diǎn)模型的鍵滲流.

表2 正方格子不同近鄰格點(diǎn)的相關(guān)參數(shù)Table 2.Parameters of different nearest neighbors on square lattice.

圖5 表1 中所示格點(diǎn)模型 pc 隨ξ 變化的對數(shù)-對數(shù)曲線Fig.5.The log-log plot of pc vs.ξ for the lattices shown in Table 1.

4 結(jié)論

綜上,為了深入分析滲流閾值與格點(diǎn)模型結(jié)構(gòu)之間的關(guān)系,本文基于高效的單團(tuán)簇生長算法,對含復(fù)雜近鄰的二維正方格子鍵滲流進(jìn)行了大量的蒙特卡羅模擬研究.得出主要結(jié)論如下:

1) 單團(tuán)簇生長算法能夠有效地模擬含復(fù)雜近鄰的格點(diǎn)模型的鍵滲流,并獲得了高精度的鍵滲流閾值.對于這類格點(diǎn)模型,由于模型中存在交叉(非相交)鍵,這在很大程度上限制了之前許多算法或方法的有效應(yīng)用,而單團(tuán)簇生長算法則提供了模擬此類格點(diǎn)模型的有效途徑.