物理学中的蒙特卡洛方法
The following article is from 中国科学院理论物理研究所 Author 夏晨
作者简介 /Profile/
夏晨，中国科学院理论物理研究所19级博士研究生，研究方向是宇宙线与暗物质直接探测，导师是周宇峰研究员。
什么是蒙特卡洛
1946年，在研究原子弹的“曼哈顿计划”中，数学家斯塔尼斯拉夫·乌拉姆在一次生病后的恢复期间玩纸牌游戏。他开始想用排列组合计算一下赢牌的概率，但是转念一想，如果“无脑”地反复玩很多次，最后数一数赢了多少次，也可以近似得到答案。当时正值第一台通用电子计算机 ENIAC 发明出来，乌拉姆马上联想到核武器研究中关于中子扩散的问题，也可以通过计算机模拟一个个中子的随机运动来研究。他将这个想法告诉了冯·诺伊曼，随后两人开始了研究[1]。为了保密，乌拉姆和冯·诺伊曼的工作需要一个代号。他们的同事 Metropolis 建议了蒙特卡洛 (Monte Carlo) 这个名字，来源于摩纳哥的一座城市蒙特卡洛，因为乌拉姆的一位叔叔喜欢向亲戚借钱去那里赌博，而赌博暗含了概率和随机性。后来蒙特卡洛逐渐从一个神秘代号演变成了一个术语，用来代指各种利用随机性来解决问题的方法[2]。
本文通过三个例子来介绍蒙特卡洛方法的典型应用方式，第一个例子是利用随机撒点计算图形面积，它经常作为蒙特卡洛方法的入门介绍，后两个例子是在物理学中的应用，分别关于统计物理和粒子物理领域。三个例子互相独立，读者可以选读感兴趣的内容。
单位圆的面积
我们从一个经典的例子开始：计算单位圆的面积。单位圆即半径为 1 的圆，它的面积我们都知道等于 。现在我们假装不知道 的值是多少，然后通过蒙特卡洛方法来求得。
单位圆可以被嵌入到边长为 2 的外接正方形中，正方形面积等于 4。如果在整个正方形中均匀地撒点，那么落在单位圆中点的个数应该正比于圆的面积。利用圆与正方形的面积之比等于落在圆与正方形内部的点数之比，我们就得到了计算单位圆面积的蒙特卡洛方法。
为了得到随机点，我们可以在纸上画一个圆和外接正方形，然后闭上眼睛用笔随机戳点。但这种做法不仅劳累，而且还不容易保证点分布的均匀性。这种工作显然适合交给计算机来做，计算机可以利用算法产生均匀分布的伪随机数来模拟撒点。
正方形中的随机点可以用坐标表示为 ，其中 和 都是从 到 之间均匀抽取的随机数。生成 个随机点后，满足条件 的点即在单位圆内，如下图所示：
1000 个随机点
假设圆内的点有 个，那么圆的面积，也就是 的估计值为
可以想象撒的点越多， 的估计值应该更准确。我们随机抽取不同数量 的点做了多次实验，结果展示在下表中：
能够看出随着 增大， 的估计值有接近真实值 的趋势，但似乎也不是 越大结果就一定越好，比如表中一万个点的结果 (3.1372) 反而比十万个点的结果 (3.14732) 更接近 。这是由于蒙特卡洛方法的本质是使用随机性，所以结果总会存在涨落，如果再进行一组实验将会得到一张不同的表。为了确认蒙特卡洛方法的结果有多可靠，需要估计结果的误差。所以我们再对每一个 都重复实验 1000 次，算出结果的平均值和标准偏差画在下图中：
蒙特卡洛模拟结果的误差随 的变化
我们看到误差随着撒点个数 的增加而减小，并且减小的规律是正比于函数 。稍作思考我们可以推导出这个规律，实际上每个随机点要么落在圆内，要么不落在圆内，落在圆内的概率等于 ，所以落在圆内的个数 满足二项分布。使用概率论的语言，随机变量 满足二项分布 ，其中 为落到圆内的概率， 为总点数。由二项分布的性质， 的平均值，或者说期望值 ，方差 。我们的蒙特卡洛模拟结果其实是对随机变量 的采样，它的期望值 ，而方差 ，所以标准差 ，正比于 。由此我们可以估计，如果要正确计算到 的小数点后第五位，至少要撒 也就是一千亿个点才有较大的把握。
这一节的方法当然不限于计算单位圆面积，任意图形的面积都可以这样计算。实际上这是蒙特卡洛数值积分最简单的应用形式，尽管计算精度随 的增大提升得比较慢，但是蒙特卡洛积分不会受到积分维度的显著影响。对于多变量的高维数值积分，蒙特卡洛方法几乎是唯一的选择。
统计物理中的伊辛模型
第二个例子是蒙特卡洛方法在统计物理中的应用。统计物理研究的是由大量微观粒子组成的宏观物体的统计规律，主要通过概率的语言来描述，所以蒙特卡洛方法应用到统计物理相当直接。统计物理的原理告诉我们，任何一个温度为 的复杂物理系统，它处在某一微观状态 的概率 满足玻尔兹曼分布
其中 表示系统的能量， 表示归一化系数，也称为配分函数，系统的宏观热力学量都可以通过它导出。 是玻尔兹曼常数，表示温度和能量之间的单位转换，理论物理工作者习惯取 ，也就是用能量单位来表示温度。系统在热扰动下其微观状态根据玻尔兹曼分布不断发生变化，但各种宏观热力学量的观测值可以用各种微观态下相应物理量的平均值来表示，比如系统的内能 ，统计物理就这样简单地和热力学联系了起来。然而原理虽然简单，实际上通常 和 这样的求和是很困难的，因为宏观物体的微观状态数量是天文数字，并且 的形式可能很复杂。这时便可以考虑蒙特卡洛方法，其想法非常直接，既然概率分布已经有了，直接按概率对系统状态进行采样，再计算样本的平均值就行了。根据概率分布采样也称为重要性采样，相比于均匀地遍历所有可能的微观态求和，重要性采样只需要抽取出少量的样本就能反应出系统的主要特征。
下面通过非常经典的伊辛 (Ising) 模型来演示。伊辛模型是关于物质磁性的简化模型，假设物质由格点上的一个个小磁针构成，每个格点 上的小磁针 只能有两个朝向，用 和 表示它的磁矩指向“上”和“下”。每个磁针只与相邻格点上的磁针存在相互作用，它们如果同向则能量为 ，反向则能量为 ，这可以用 来概括。在不存在外磁场的情况下，系统的总能量等于所有相邻格点之间能量的和 这里符号 表示 和 是相邻格点。对于伊辛模型，系统的一个微观状态 ，就是所有格点上磁针的某一种排列方式， ， 表示格点个数。另外我们关心单位格点上的平均磁矩 相应的总磁矩即为 。系统处于平衡状态时的单位格点磁矩则计算为 。
一维格点伊辛模型的严格解在 1925 年已由伊辛本人求出，二维无外磁场伊辛模型由昂萨格于 1944 年求解，杨振宁先生在 1952 年发表了二维伊辛模型磁化强度的推导，而三维伊辛模型的精确解到目前为止还没有得到。二维伊辛模型的精确解表明系统存在铁磁相变，对应的相变临界温度为 当温度高于临界温度 时，系统处于顺磁相，平均磁矩为零，而温度低于 时处于铁磁相，平均磁矩不为零。类似平均磁矩这种在一个相中为零，另一个相中不为零的量可以叫作序参量。像伊辛模型这种序参量在临界温度处连续变化的相变称为连续相变。下面我们介绍如何用蒙特卡洛方法模拟二维伊辛模型。
在这里，蒙特卡洛方法的核心问题是：如何通过玻尔兹曼分布对系统微观态采样？答案之一是采用马尔科夫链蒙特卡洛 (Markov Chain Monte Carlo, MCMC) 方法。所谓马尔科夫链是指一个离散的随机过程，系统从初始状态开始，由一个状态按照一定的概率跳转到另一个状态，每一次跳转的概率只与当前的状态有关，与此前的历史无关，多次跳转后形成一条链。MCMC 的要点在于巧妙地设置跳转概率，使得链条上的样本满足目标分布。
构造马尔科夫链的一个传统算法是 Metropolis 算法，对于伊辛模型其基本流程为：
基于当前状态 ，随机选择一个格点将其磁矩翻转，得到一个新的状态 。
如果能量 ，则接受新状态 ；如果 则按概率 决定是否接受，若拒绝 则保留当前状态 作为新的状态。
上面两个步骤反复迭代就能得到一系列的系统构型。可以看到如果新状态的能量更低，则一定会被接受，系统可以很快地向能量更低的状态移动，同时也有机会跳转到能量更高的状态，最后得到的样本可以证明满足玻尔兹曼分布。此外 MCMC 方法不需确切地知道概率分布的归一化系数，从而避免了配分函数的计算。
下图以 的二维正方格点为例，展示伊辛模型的 MCMC 模拟结果[3]：
二维伊辛模型 MCMC 模拟
和 的格点分别显示为白色和绿色， 仅表示马尔科夫链的迭代过程，并不代表真实的时间，并且每两帧之间间隔了 100 次迭代。我们看到在高温下， 时，格点上的磁矩排布是杂乱无章的，但温度降到临界温度 后，系统自动地出现磁矩指向一致的磁性团块。注意到系统并没有施加外磁场，这是一种自发磁化行为。
在不同温度下进行模拟，计算样本的单位格点磁矩，可以探究序参量随温度的变化关系，研究系统的相变特性，如下图所示：
二维伊辛模型单位格点磁矩随温度的变化
这里我们为了简便，取了 。MCMC 在每一个温度下的结果都由 1 万个样本点计算得到。我们通过 的格点模拟得到了与精确解符合得还不错的结果，而精确解是在热力学极限，也就是格点数量 的条件下得到的。三维伊辛模型虽然没有求出精确解，但同样可以通过蒙特卡洛方法来研究。
本节主要介绍的是如何通过给定概率分布来进行采样的蒙特卡洛方法，因为按概率采样本质上需要去寻找概率分布的峰的位置，所以这种类型的应用可以拓展到求函数极值的最优化问题，同样也是蒙特卡洛方法最重要的应用场景之一。
暗物质在地球内部的运动
第三个例子我们进入粒子物理领域，以暗物质粒子在地球内部的运动为例，介绍随机游走过程的蒙特卡洛模拟。
许多天文观测发现，宇宙中我们熟悉的可见物质，如恒星、行星、星云等等，不足以提供足够的引力来解释观测到的物质运动方式，例如星系的旋转速度太大，星系团内的星系运动太快，光线在引力场附近的弯曲过强等等，因此提出可能存在看不见的暗物质，来弥补缺失的质量。并且现代宇宙学根据宇宙微波背景辐射的观测数据，推测出暗物质应占宇宙物质总量的 左右，这意味着人类对于宇宙的认识可能还只在冰山一角，探索暗物质的本质是当前物理学的前沿课题。
我们已经知道普通物质由原子构成，原子又由基本粒子构成，那么暗物质是否也是某种未知的基本粒子呢？粒子物理学家们提出了众多的粒子模型，为了能够在宇宙中产生，这些模型或多或少都要求暗物质与普通物质之间存在除引力之外的相互作用，这就为暗物质粒子的实验探测带来了可能。目前世界各地建立起了大量的暗物质直接探测实验，清华大学主导的 CDEX 实验和上海交通大学主导的 PandaX 实验就是其中的佼佼者，它们位于四川锦屏山隧道中的锦屏地下实验室，垂直埋深达到 2.4 千米，是世界上最深的地下实验室。之所以建造在地下，是因为暗物质直接探测实验的目标是寻找暗物质粒子与靶材料之间的碰撞事件，需要利用厚厚的土层和岩石来屏蔽高能宇宙线的干扰。到目前为止，还没有探测到暗物质的明确信号，实验技术仍在不断发展之中。
实验室建造在地下能够屏蔽背景的同时，如果暗物质与物质相互作用不太弱的话，也有可能屏蔽掉我们想要探测的暗物质粒子，这个问题就可以使用蒙特卡洛模拟方法来研究。暗物质粒子从地表进入到地球内部后的运动可以看作随机游走的过程，我们只要模拟大量粒子的运动轨迹，就能重建出在地下实验室中的暗物质分布。
银河系可能被一个巨大的暗物质晕包围，其中暗物质粒子的速度满足麦克斯韦分布，平均速度大致和银河系中星体的运动速度相当，约为 ，称为标准暗晕模型 (Standard Halo Model, SHM)。暗物质粒子在地表处的初始速度将通过这一速度分布抽样得到，随后的随机游走则由两个步骤反复迭代进行：
自由传播：粒子在发生碰撞之前沿直线自由传播一段距离，距离的长度称为自由程。自由程满足特定的概率分布，其平均值即平均自由程，由粒子与地球内部元素相互作用强度和地球的密度确定。模拟中首先计算平均自由程，然后自由程根据相应的概率分布抽样得到。
碰撞：自由运动结束后暗物质粒子与地球内部元素发生碰撞，碰撞将导致暗物质粒子损失一部分能量而减速，并且运动方向改变。新的速度大小和方向由相互作用的具体形式按概率抽样确定，随后重复进行下一段自由传播。
通过这样一步一步的随机过程，可以模拟出暗物质粒子折线形式的运动轨迹，如下图所示[4]：
暗物质粒子轨迹模拟，由坐标原点处垂直向下出发
模拟大量轨迹之后，收集每个粒子经过实验室深度时的速度，就可以统计出地下实验室中的暗物质速度分布，作为直接探测实验信号分析的重要输入条件。下图展示了质量约等于质子质量的暗物质粒子，在地下 2.4 km 深度的速度分布模拟结果：
地下暗物质粒子速度分布
图中的速度分布归一化到暗物质粒子数密度，即曲线下的面积代表暗物质粒子数量的相对大小。SHM 标记的虚线表示地球外部的标准暗晕速度分布。不同颜色的实线代表不同的相互作用强度，通过暗物质粒子与质子的散射截面 来刻画。散射截面可以理解为如果把暗物质粒子看作小球的话，它的横截面积的大小。截面越大，暗物质粒子在地球内部碰撞的机会就更多，从而损失更多能量，使得高速运动的暗物质粒子数量变少，并在低速部分堆积。由于探测器需要暗物质粒子具有一定的动能才能触发，这就使得相互作用较强的暗物质粒子可能反而探测不到[5]。
本节介绍的通过随机游走模拟暗物质粒子运动的方法，本质上等同于曼哈顿计划中代号为蒙特卡洛的核裂变反应的中子扩散模拟，此外高能粒子对撞机中的探测器模拟也是采用类似的方法，只是在这些应用中，碰撞过程中会产生的大量次级粒子需要记录和追踪。
结语
蒙特卡洛方法并不特指某种特定的算法，而是对利用随机性来解决问题的方法的统称，我们通过几个例子展示了蒙特卡洛方法的典型应用。物理学，以及所有的科学，都是以实验为基础的，使用蒙特卡洛方法相当于在计算机中进行模拟实验。尽管新物理只可能在真实的实验中发现，模拟实验只能输入已知的物理定律，但正如我们看到的，蒙特卡洛方法可以为理论推导提供佐证，可以用简单的方式解决困难的问题。在开始真实的实验之前，蒙特卡洛模拟也是检验实验方案可行性的重要手段，特别是对于实验成本非常高昂的大科学装置，如暗物质探测器，高能粒子对撞机等等。同样的，在各种工程和应用领域如航空航天工业，军事武器研发等等，蒙特卡洛模拟都是必要的环节。对于我们个人来说，利用一台小小的电脑，就能研究磁性系统的相变，观察暗物质粒子的运动，甚至模拟宇宙的演化、星系的形成这些不可能在现实中直接观测的过程，这本身就是非常有趣的事情。
参考文献：
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.itp.ac.cn/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.
微信号｜ITP-CAS
开放 交融 求真 创新
· 中科院理论物理研究所 ·