|
作者簡介/Profile/ 夏晨，中國科學(xué)院理論物理研究所19級博士研究生，研究方向是宇宙線與暗物質(zhì)直接探測，導(dǎo)師是周宇峰研究員。 什么是蒙特卡洛 1946年，在研究原子彈的“曼哈頓計劃”中，數(shù)學(xué)家斯塔尼斯拉夫·烏拉姆在一次生病后的恢復(fù)期間玩紙牌游戲。他開始想用排列組合計算一下贏牌的概率，但是轉(zhuǎn)念一想，如果“無腦”地反復(fù)玩很多次，最后數(shù)一數(shù)贏了多少次，也可以近似得到答案。當(dāng)時正值第一臺通用電子計算機(jī) ENIAC 發(fā)明出來，烏拉姆馬上聯(lián)想到核武器研究中關(guān)于中子擴(kuò)散的問題，也可以通過計算機(jī)模擬一個個中子的隨機(jī)運(yùn)動來研究。他將這個想法告訴了馮·諾伊曼，隨后兩人開始了研究[1]。為了保密，烏拉姆和馮·諾伊曼的工作需要一個代號。他們的同事 Metropolis 建議了蒙特卡洛 (Monte Carlo) 這個名字，來源于摩納哥的一座城市蒙特卡洛，因為烏拉姆的一位叔叔喜歡向親戚借錢去那里賭博，而賭博暗含了概率和隨機(jī)性。后來蒙特卡洛逐漸從一個神秘代號演變成了一個術(shù)語，用來代指各種利用隨機(jī)性來解決問題的方法[2]。 本文通過三個例子來介紹蒙特卡洛方法的典型應(yīng)用方式，第一個例子是利用隨機(jī)撒點(diǎn)計算圖形面積，它經(jīng)常作為蒙特卡洛方法的入門介紹，后兩個例子是在物理學(xué)中的應(yīng)用，分別關(guān)于統(tǒng)計物理和粒子物理領(lǐng)域。三個例子互相獨(dú)立，讀者可以選讀感興趣的內(nèi)容。 單位圓的面積 我們從一個經(jīng)典的例子開始：計算單位圓的面積。單位圓即半徑為 1 的圓，它的面積我們都知道等于 。現(xiàn)在我們假裝不知道 的值是多少，然后通過蒙特卡洛方法來求得。 單位圓可以被嵌入到邊長為 2 的外接正方形中，正方形面積等于 4。如果在整個正方形中均勻地撒點(diǎn)，那么落在單位圓中點(diǎn)的個數(shù)應(yīng)該正比于圓的面積。利用圓與正方形的面積之比等于落在圓與正方形內(nèi)部的點(diǎn)數(shù)之比，我們就得到了計算單位圓面積的蒙特卡洛方法。 為了得到隨機(jī)點(diǎn)，我們可以在紙上畫一個圓和外接正方形，然后閉上眼睛用筆隨機(jī)戳點(diǎn)。但這種做法不僅勞累，而且還不容易保證點(diǎn)分布的均勻性。這種工作顯然適合交給計算機(jī)來做，計算機(jī)可以利用算法產(chǎn)生均勻分布的偽隨機(jī)數(shù)來模擬撒點(diǎn)。 正方形中的隨機(jī)點(diǎn)可以用坐標(biāo)表示為 ，其中 和 都是從 到 之間均勻抽取的隨機(jī)數(shù)。生成 個隨機(jī)點(diǎn)后，滿足條件 的點(diǎn)即在單位圓內(nèi)，如下圖所示：  1000 個隨機(jī)點(diǎn) 假設(shè)圓內(nèi)的點(diǎn)有 個，那么圓的面積，也就是 的估計值為 可以想象撒的點(diǎn)越多， 的估計值應(yīng)該更準(zhǔn)確。我們隨機(jī)抽取不同數(shù)量 的點(diǎn)做了多次實(shí)驗，結(jié)果展示在下表中： 能夠看出隨著 增大， 的估計值有接近真實(shí)值 的趨勢，但似乎也不是 越大結(jié)果就一定越好，比如表中一萬個點(diǎn)的結(jié)果 (3.1372) 反而比十萬個點(diǎn)的結(jié)果 (3.14732) 更接近 。這是由于蒙特卡洛方法的本質(zhì)是使用隨機(jī)性，所以結(jié)果總會存在漲落，如果再進(jìn)行一組實(shí)驗將會得到一張不同的表。為了確認(rèn)蒙特卡洛方法的結(jié)果有多可靠，需要估計結(jié)果的誤差。所以我們再對每一個 都重復(fù)實(shí)驗 1000 次，算出結(jié)果的平均值和標(biāo)準(zhǔn)偏差畫在下圖中：  蒙特卡洛模擬結(jié)果的誤差隨 的變化 我們看到誤差隨著撒點(diǎn)個數(shù) 的增加而減小，并且減小的規(guī)律是正比于函數(shù) 。稍作思考我們可以推導(dǎo)出這個規(guī)律，實(shí)際上每個隨機(jī)點(diǎn)要么落在圓內(nèi)，要么不落在圓內(nèi)，落在圓內(nèi)的概率等于 ，所以落在圓內(nèi)的個數(shù) 滿足二項分布。使用概率論的語言，隨機(jī)變量 滿足二項分布 ，其中 為落到圓內(nèi)的概率， 為總點(diǎn)數(shù)。由二項分布的性質(zhì)， 的平均值，或者說期望值 ，方差 。我們的蒙特卡洛模擬結(jié)果其實(shí)是對隨機(jī)變量 的采樣，它的期望值 ，而方差 ，所以標(biāo)準(zhǔn)差 ，正比于 。由此我們可以估計，如果要正確計算到 的小數(shù)點(diǎn)后第五位，至少要撒 也就是一千億個點(diǎn)才有較大的把握。 這一節(jié)的方法當(dāng)然不限于計算單位圓面積，任意圖形的面積都可以這樣計算。實(shí)際上這是蒙特卡洛數(shù)值積分最簡單的應(yīng)用形式，盡管計算精度隨 的增大提升得比較慢，但是蒙特卡洛積分不會受到積分維度的顯著影響。對于多變量的高維數(shù)值積分，蒙特卡洛方法幾乎是唯一的選擇。 統(tǒng)計物理中的伊辛模型 第二個例子是蒙特卡洛方法在統(tǒng)計物理中的應(yīng)用。統(tǒng)計物理研究的是由大量微觀粒子組成的宏觀物體的統(tǒng)計規(guī)律，主要通過概率的語言來描述，所以蒙特卡洛方法應(yīng)用到統(tǒng)計物理相當(dāng)直接。統(tǒng)計物理的原理告訴我們，任何一個溫度為 的復(fù)雜物理系統(tǒng)，它處在某一微觀狀態(tài) 的概率 滿足玻爾茲曼分布 其中 表示系統(tǒng)的能量， 表示歸一化系數(shù)，也稱為配分函數(shù)，系統(tǒng)的宏觀熱力學(xué)量都可以通過它導(dǎo)出。 是玻爾茲曼常數(shù)，表示溫度和能量之間的單位轉(zhuǎn)換，理論物理工作者習(xí)慣取 ，也就是用能量單位來表示溫度。系統(tǒng)在熱擾動下其微觀狀態(tài)根據(jù)玻爾茲曼分布不斷發(fā)生變化，但各種宏觀熱力學(xué)量的觀測值可以用各種微觀態(tài)下相應(yīng)物理量的平均值來表示，比如系統(tǒng)的內(nèi)能 ，統(tǒng)計物理就這樣簡單地和熱力學(xué)聯(lián)系了起來。然而原理雖然簡單，實(shí)際上通常 和 這樣的求和是很困難的，因為宏觀物體的微觀狀態(tài)數(shù)量是天文數(shù)字，并且 的形式可能很復(fù)雜。這時便可以考慮蒙特卡洛方法，其想法非常直接，既然概率分布已經(jīng)有了，直接按概率對系統(tǒng)狀態(tài)進(jìn)行采樣，再計算樣本的平均值就行了。根據(jù)概率分布采樣也稱為重要性采樣，相比于均勻地遍歷所有可能的微觀態(tài)求和，重要性采樣只需要抽取出少量的樣本就能反應(yīng)出系統(tǒng)的主要特征。 下面通過非常經(jīng)典的伊辛 (Ising) 模型來演示。伊辛模型是關(guān)于物質(zhì)磁性的簡化模型，假設(shè)物質(zhì)由格點(diǎn)上的一個個小磁針構(gòu)成，每個格點(diǎn) 上的小磁針 只能有兩個朝向，用 和 表示它的磁矩指向“上”和“下”。每個磁針只與相鄰格點(diǎn)上的磁針存在相互作用，它們?nèi)绻騽t能量為 ，反向則能量為 ，這可以用 來概括。在不存在外磁場的情況下，系統(tǒng)的總能量等于所有相鄰格點(diǎn)之間能量的和 這里符號 表示 和 是相鄰格點(diǎn)。對于伊辛模型，系統(tǒng)的一個微觀狀態(tài) ，就是所有格點(diǎn)上磁針的某一種排列方式， ， 表示格點(diǎn)個數(shù)。另外我們關(guān)心單位格點(diǎn)上的平均磁矩 相應(yīng)的總磁矩即為 。系統(tǒng)處于平衡狀態(tài)時的單位格點(diǎn)磁矩則計算為 。 一維格點(diǎn)伊辛模型的嚴(yán)格解在 1925 年已由伊辛本人求出，二維無外磁場伊辛模型由昂薩格于 1944 年求解，楊振寧先生在 1952 年發(fā)表了二維伊辛模型磁化強(qiáng)度的推導(dǎo)，而三維伊辛模型的精確解到目前為止還沒有得到。二維伊辛模型的精確解表明系統(tǒng)存在鐵磁相變，對應(yīng)的相變臨界溫度為 當(dāng)溫度高于臨界溫度 時，系統(tǒng)處于順磁相，平均磁矩為零，而溫度低于 時處于鐵磁相，平均磁矩不為零。類似平均磁矩這種在一個相中為零，另一個相中不為零的量可以叫作序參量。像伊辛模型這種序參量在臨界溫度處連續(xù)變化的相變稱為連續(xù)相變。下面我們介紹如何用蒙特卡洛方法模擬二維伊辛模型。 在這里，蒙特卡洛方法的核心問題是：如何通過玻爾茲曼分布對系統(tǒng)微觀態(tài)采樣？答案之一是采用馬爾科夫鏈蒙特卡洛 (Markov Chain Monte Carlo, MCMC) 方法。所謂馬爾科夫鏈?zhǔn)侵敢粋€離散的隨機(jī)過程，系統(tǒng)從初始狀態(tài)開始，由一個狀態(tài)按照一定的概率跳轉(zhuǎn)到另一個狀態(tài)，每一次跳轉(zhuǎn)的概率只與當(dāng)前的狀態(tài)有關(guān)，與此前的歷史無關(guān)，多次跳轉(zhuǎn)后形成一條鏈。MCMC 的要點(diǎn)在于巧妙地設(shè)置跳轉(zhuǎn)概率，使得鏈條上的樣本滿足目標(biāo)分布。 構(gòu)造馬爾科夫鏈的一個傳統(tǒng)算法是 Metropolis 算法，對于伊辛模型其基本流程為： 基于當(dāng)前狀態(tài) ，隨機(jī)選擇一個格點(diǎn)將其磁矩翻轉(zhuǎn)，得到一個新的狀態(tài) 。 如果能量 ，則接受新狀態(tài) ；如果 則按概率 決定是否接受，若拒絕 則保留當(dāng)前狀態(tài) 作為新的狀態(tài)。 上面兩個步驟反復(fù)迭代就能得到一系列的系統(tǒng)構(gòu)型。可以看到如果新狀態(tài)的能量更低，則一定會被接受，系統(tǒng)可以很快地向能量更低的狀態(tài)移動，同時也有機(jī)會跳轉(zhuǎn)到能量更高的狀態(tài)，最后得到的樣本可以證明滿足玻爾茲曼分布。此外 MCMC 方法不需確切地知道概率分布的歸一化系數(shù)，從而避免了配分函數(shù)的計算。 下圖以 的二維正方格點(diǎn)為例，展示伊辛模型的 MCMC 模擬結(jié)果[3]：   二維伊辛模型 MCMC 模擬 和 的格點(diǎn)分別顯示為白色和綠色， 僅表示馬爾科夫鏈的迭代過程，并不代表真實(shí)的時間，并且每兩幀之間間隔了 100 次迭代。我們看到在高溫下， 時，格點(diǎn)上的磁矩排布是雜亂無章的，但溫度降到臨界溫度 后，系統(tǒng)自動地出現(xiàn)磁矩指向一致的磁性團(tuán)塊。注意到系統(tǒng)并沒有施加外磁場，這是一種自發(fā)磁化行為。 在不同溫度下進(jìn)行模擬，計算樣本的單位格點(diǎn)磁矩，可以探究序參量隨溫度的 變化關(guān)系，研究系統(tǒng)的相變特性，如下圖所示：  二維伊辛模型單位格點(diǎn)磁矩隨溫度的變化 這里我們?yōu)榱撕啽?，取?/div> 。MCMC 在每一個溫度下的結(jié)果都由 1 萬個樣本點(diǎn)計算得到。我們通過 的格點(diǎn)模擬得到了與精確解符合得還不錯的結(jié)果，而精確解是在熱力學(xué)極限，也就是格點(diǎn)數(shù)量 的條件下得到的。三維伊辛模型雖然沒有求出精確解，但同樣可以通過蒙特卡洛方法來研究。 本節(jié)主要介紹的是如何通過給定概率分布來進(jìn)行采樣的蒙特卡洛方法，因為按概率采樣本質(zhì)上需要去尋找概率分布的峰的位置，所以這種類型的應(yīng)用可以拓展到求函數(shù)極值的最優(yōu)化問題，同樣也是蒙特卡洛方法最重要的應(yīng)用場景之一。 暗物質(zhì)在地球內(nèi)部的運(yùn)動 第三個例子我們進(jìn)入粒子物理領(lǐng)域，以暗物質(zhì)粒子在地球內(nèi)部的運(yùn)動為例，介紹隨機(jī)游走過程的蒙特卡洛模擬。 許多天文觀測發(fā)現(xiàn)，宇宙中我們熟悉的可見物質(zhì)，如恒星、行星、星云等等，不足以提供足夠的引力來解釋觀測到的物質(zhì)運(yùn)動方式，例如星系的旋轉(zhuǎn)速度太大，星系團(tuán)內(nèi)的星系運(yùn)動太快，光線在引力場附近的彎曲過強(qiáng)等等，因此提出可能存在看不見的暗物質(zhì)，來彌補(bǔ)缺失的質(zhì)量。并且現(xiàn)代宇宙學(xué)根據(jù)宇宙微波背景輻射的觀測數(shù)據(jù)，推測出暗物質(zhì)應(yīng)占宇宙物質(zhì)總量的 左右，這意味著人類對于宇宙的認(rèn)識可能還只在冰山一角，探索暗物質(zhì)的本質(zhì)是當(dāng)前物理學(xué)的前沿課題。 我們已經(jīng)知道普通物質(zhì)由原子構(gòu)成，原子又由基本粒子構(gòu)成，那么暗物質(zhì)是否也是某種未知的基本粒子呢？粒子物理學(xué)家們提出了眾多的粒子模型，為了能夠在宇宙中產(chǎn)生，這些模型或多或少都要求暗物質(zhì)與普通物質(zhì)之間存在除引力之外的相互作用，這就為暗物質(zhì)粒子的實(shí)驗探測帶來了可能。目前世界各地建立起了大量的暗物質(zhì)直接探測實(shí)驗，清華大學(xué)主導(dǎo)的 CDEX 實(shí)驗和上海交通大學(xué)主導(dǎo)的 PandaX 實(shí)驗就是其中的佼佼者，它們位于四川錦屏山隧道中的錦屏地下實(shí)驗室，垂直埋深達(dá)到 2.4 千米，是世界上最深的地下實(shí)驗室。之所以建造在地下，是因為暗物質(zhì)直接探測實(shí)驗的目標(biāo)是尋找暗物質(zhì)粒子與靶材料之間的碰撞事件，需要利用厚厚的土層和巖石來屏蔽高能宇宙線的干擾。到目前為止，還沒有探測到暗物質(zhì)的明確信號，實(shí)驗技術(shù)仍在不斷發(fā)展之中。 實(shí)驗室建造在地下能夠屏蔽背景的同時，如果暗物質(zhì)與物質(zhì)相互作用不太弱的話，也有可能屏蔽掉我們想要探測的暗物質(zhì)粒子，這個問題就可以使用蒙特卡洛模擬方法來研究。暗物質(zhì)粒子從地表進(jìn)入到地球內(nèi)部后的運(yùn)動可以看作隨機(jī)游走的過程，我們只要模擬大量粒子的運(yùn)動軌跡，就能重建出在地下實(shí)驗室中的暗物質(zhì)分布。 銀河系可能被一個巨大的暗物質(zhì)暈包圍，其中暗物質(zhì)粒子的速度滿足麥克斯韋分布，平均速度大致和銀河系中星體的運(yùn)動速度相當(dāng)，約為 ，稱為標(biāo)準(zhǔn)暗暈?zāi)Ｐ?/span>(Standard Halo Model, SHM)。暗物質(zhì)粒子在地表處的初始速度將通過這一速度分布抽樣得到，隨后的隨機(jī)游走則由兩個步驟反復(fù)迭代進(jìn)行： 自由傳播：粒子在發(fā)生碰撞之前沿直線自由傳播一段距離，距離的長度稱為自由程。自由程滿足特定的概率分布，其平均值即平均自由程，由粒子與地球內(nèi)部元素相互作用強(qiáng)度和地球的密度確定。模擬中首先計算平均自由程，然后自由程根據(jù)相應(yīng)的概率分布抽樣得到。 碰撞：自由運(yùn)動結(jié)束后暗物質(zhì)粒子與地球內(nèi)部元素發(fā)生碰撞，碰撞將導(dǎo)致暗物質(zhì)粒子損失一部分能量而減速，并且運(yùn)動方向改變。新的速度大小和方向由相互作用的具體形式按概率抽樣確定，隨后重復(fù)進(jìn)行下一段自由傳播。 通過這樣一步一步的隨機(jī)過程，可以模擬出暗物質(zhì)粒子折線形式的運(yùn)動軌跡，如下圖所示[4]：  暗物質(zhì)粒子軌跡模擬，由坐標(biāo)原點(diǎn)處垂直向下出發(fā) 模擬大量軌跡之后，收集每個粒子經(jīng)過實(shí)驗室深度時的速度，就可以統(tǒng)計出地下實(shí)驗室中的暗物質(zhì)速度分布，作為直接探測實(shí)驗信號分析的重要輸入條件。下圖展示了質(zhì)量約等于質(zhì)子質(zhì)量的暗物質(zhì)粒子，在地下 2.4 km 深度的速度分布模擬結(jié)果：  地下暗物質(zhì)粒子速度分布 圖中的速度分布?xì)w一化到暗物質(zhì)粒子數(shù)密度，即曲線下的面積代表暗物質(zhì)粒子數(shù)量的相對大小。SHM 標(biāo)記的虛線表示地球外部的標(biāo)準(zhǔn)暗暈速度分布。不同顏色的實(shí)線代表不同的相互作用強(qiáng)度，通過暗物質(zhì)粒子與質(zhì)子的散射截面 來刻畫。散射截面可以理解為如果把暗物質(zhì)粒子看作小球的話，它的橫截面積的大小。截面越大，暗物質(zhì)粒子在地球內(nèi)部碰撞的機(jī)會就更多，從而損失更多能量，使得高速運(yùn)動的暗物質(zhì)粒子數(shù)量變少，并在低速部分堆積。由于探測器需要暗物質(zhì)粒子具有一定的動能才能觸發(fā)，這就使得相互作用較強(qiáng)的暗物質(zhì)粒子可能反而探測不到[5]。 本節(jié)介紹的通過隨機(jī)游走模擬暗物質(zhì)粒子運(yùn)動的方法，本質(zhì)上等同于曼哈頓計劃中代號為蒙特卡洛的核裂變反應(yīng)的中子擴(kuò)散模擬，此外高能粒子對撞機(jī)中的探測器模擬也是采用類似的方法，只是在這些應(yīng)用中，碰撞過程中會產(chǎn)生的大量次級粒子需要記錄和追蹤。 結(jié)語 蒙特卡洛方法并不特指某種特定的算法，而是對利用隨機(jī)性來解決問題的方法的統(tǒng)稱，我們通過幾個例子展示了蒙特卡洛方法的典型應(yīng)用。物理學(xué)，以及所有的科學(xué)，都是以實(shí)驗為基礎(chǔ)的，使用蒙特卡洛方法相當(dāng)于在計算機(jī)中進(jìn)行模擬實(shí)驗。盡管新物理只可能在真實(shí)的實(shí)驗中發(fā)現(xiàn)，模擬實(shí)驗只能輸入已知的物理定律，但正如我們看到的，蒙特卡洛方法可以為理論推導(dǎo)提供佐證，可以用簡單的方式解決困難的問題。在開始真實(shí)的實(shí)驗之前，蒙特卡洛模擬也是檢驗實(shí)驗方案可行性的重要手段，特別是對于實(shí)驗成本非常高昂的大科學(xué)裝置，如暗物質(zhì)探測器，高能粒子對撞機(jī)等等。同樣的，在各種工程和應(yīng)用領(lǐng)域如航空航天工業(yè)，軍事武器研發(fā)等等，蒙特卡洛模擬都是必要的環(huán)節(jié)。對于我們個人來說，利用一臺小小的電腦，就能研究磁性系統(tǒng)的相變，觀察暗物質(zhì)粒子的運(yùn)動，甚至模擬宇宙的演化、星系的形成這些不可能在現(xiàn)實(shí)中直接觀測的過程，這本身就是非常有趣的事情。 參考文獻(xiàn)： 1. Eckhardt, Roger (1987). “Stan Ulam, John von Neumann, and the Monte Carlo method” . Los Alamos Science (15): 131–137. 2. Metropolis, N. (1987). “The beginning of the Monte Carlo method”. Los Alamos Science (1987 Special Issue dedicated to Stanislaw Ulam): 125–130. 3. 使用 Julia 語言 Ising2D.jl 程序包制作。 4. 使用作者編寫的 darkprop 程序包計算。http://yfzhou./darkprop. 5. Emken T, Kouvaris C. “How blind are underground and surface detectors to strongly interacting Dark Matter?” Phys. Rev. D, 2018, 97(11): 115047. arXiv:1802.04764.
|
|
|
