物理学中的蒙特卡洛方法

luyuanhong

物理学中的蒙特卡洛方法

撰文 | 夏晨

什么是蒙特卡洛

1946 年，在研究原子弹的“曼哈顿计划”中，数学家斯塔尼斯拉夫·乌拉姆在一次生病后的恢复期间玩纸牌游戏。他开始想用排列组合计算一下赢牌的概率，但是转念一想，如果“无脑”地反复玩很多次，最后数一数赢了多少次，也可以近似得到答案。当时正值第一台通用电子计算机 ENIAC 发明出来，乌拉姆马上联想到核武器研究中关于中子扩散的问题，也可以通过计算机模拟一个个中子的随机运动来研究。他将这个想法告诉了冯·诺伊曼，随后两人开始了研究[1]。为了保密，乌拉姆和冯·诺伊曼的工作需要一个代号。他们的同事 Metropolis 建议了蒙特卡洛 (Monte Carlo) 这个名字，来源于摩纳哥的一座城市蒙特卡洛，因为乌拉姆的一位叔叔喜欢向亲戚借钱去那里赌博，而赌博暗含了概率和随机性。后来蒙特卡洛逐渐从一个神秘代号演变成了一个术语，用来代指各种利用随机性来解决问题的方法[2]。

本文通过三个例子来介绍蒙特卡洛方法的典型应用方式，第一个例子是利用随机撒点计算图形面积，它经常作为蒙特卡洛方法的入门介绍，后两个例子是在物理学中的应用，分别关于统计物理和粒子物理领域。三个例子互相独立，读者可以选读感兴趣的内容。

单位圆的面积

我们从一个经典的例子开始：计算单位圆的面积。单位圆即半径为 1 的圆，它的面积我们都知道等于 π 。现在我们假装不知道 π 的值是多少，然后通过蒙特卡洛方法来求得。

单位圆可以被嵌入到边长为 2 的外接正方形中，正方形面积等于 4 。如果在整个正方形中均匀地撒点，那么落在单位圆中点的个数应该正比于圆的面积。利用圆与正方形的面积之比等于落在圆与正方形内部的点数之比，我们就得到了计算单位圆面积的蒙特卡洛方法。

为了得到随机点，我们可以在纸上画一个圆和外接正方形，然后闭上眼睛用笔随机戳点。但这种做法不仅劳累，而且还不容易保证点分布的均匀性。这种工作显然适合交给计算机来做，计算机可以利用算法产生均匀分布的伪随机数来模拟撒点。

正方形中的随机点可以用坐标表示为 (x,y) ，其中 x 和 y 都是从 -1 到 1 之间均匀抽取的随机数。生成 N 个随机点后，满足条件 x^2+y^2＜1 的点即在单位圆内，如下图所示：


1000 个随机点

假设圆内的点有 M 个，那么圆的面积，也就是 π 的估计值为 π≈4×M/N 。

可以想象撒的点越多，π 的估计值应该更准确。我们随机抽取不同数量 N 的点做了多次实验，结果展示在下表中：



能够看出随着 N 增大，π 的估计值有接近真实值 3.1415926… 的趋势，但似乎也不是 N 越大结果就一定越好，比如表中一万个点的结果 (3.1372) 反而比十万个点的结果 (3.14732) 更接近 π 。这是由于蒙特卡洛方法的本质是使用随机性，所以结果总会存在涨落，如果再进行一组实验将会得到一张不同的表。为了确认蒙特卡洛方法的结果有多可靠，需要估计结果的误差。所以我们再对每一个 N 都重复实验 1000 次，算出结果的平均值和标准偏差画在下图中：


蒙特卡洛模拟结果的误差随 N 的变化



这一节的方法当然不限于计算单位圆面积，任意图形的面积都可以这样计算。实际上这是蒙特卡洛数值积分最简单的应用形式，尽管计算精度随 N 的增大提升得比较慢，但是蒙特卡洛积分不会受到积分维度的显著影响。对于多变量的高维数值积分，蒙特卡洛方法几乎是唯一的选择。

统计物理中的伊辛模型



下图以 100×100 的二维正方格点为例，展示伊辛模型的 MCMC 模拟结果[3]：


二维伊辛模型 MCMC 模拟

si=+1 和 -1 的格点分别显示为白色和绿色，t 仅表示马尔科夫链的迭代过程，并不代表真实的时间，并且每两帧之间间隔了 100 次迭代。我们看到在高温下，T=2Tc 时，格点上的磁矩排布是杂乱无章的，但温度降到临界温度 Tc 后，系统自动地出现磁矩指向一致的磁性团块。注意到系统并没有施加外磁场，这是一种自发磁化行为。

在不同温度下进行模拟，计算样本的单位格点磁矩，可以探究序参量随温度的变化关系，研究系统的相变特性，如下图所示：


二维伊辛模型单位格点磁矩随温度的变化

这里我们为了简便，取了 J=1 ，k=1 。MCMC 在每一个温度下的结果都由 1 万个样本点计算得到。我们通过  100×100 的格点模拟得到了与精确解符合得还不错的结果，而精确解是在热力学极限，也就是格点数量 N→∞ 的条件下得到的。三维伊辛模型虽然没有求出精确解，但同样可以通过蒙特卡洛方法来研究。

本节主要介绍的是如何通过给定概率分布来进行采样的蒙特卡洛方法，因为按概率采样本质上需要去寻找概率分布的峰的位置，所以这种类型的应用可以拓展到求函数极值的最优化问题，同样也是蒙特卡洛方法最重要的应用场景之一。

暗物质在地球内部的运动

第三个例子我们进入粒子物理领域，以暗物质粒子在地球内部的运动为例，介绍随机游走过程的蒙特卡洛模拟。

许多天文观测发现，宇宙中我们熟悉的可见物质，如恒星、行星、星云等等，不足以提供足够的引力来解释观测到的物质运动方式，例如星系的旋转速度太大，星系团内的星系运动太快，光线在引力场附近的弯曲过强等等，因此提出可能存在看不见的暗物质，来弥补缺失的质量。并且现代宇宙学根据宇宙微波背景辐射的观测数据，推测出暗物质应占宇宙物质总量的  左右，这意味着人类对于宇宙的认识可能还只在冰山一角，探索暗物质的本质是当前物理学的前沿课题。

我们已经知道普通物质由原子构成，原子又由基本粒子构成，那么暗物质是否也是某种未知的基本粒子呢？粒子物理学家们提出了众多的粒子模型，为了能够在宇宙中产生，这些模型或多或少都要求暗物质与普通物质之间存在除引力之外的相互作用，这就为暗物质粒子的实验探测带来了可能。目前世界各地建立起了大量的暗物质直接探测实验，清华大学主导的 CDEX 实验和上海交通大学主导的 PandaX 实验就是其中的佼佼者，它们位于四川锦屏山隧道中的锦屏地下实验室，垂直埋深达到 2.4 千米，是世界上最深的地下实验室。之所以建造在地下，是因为暗物质直接探测实验的目标是寻找暗物质粒子与靶材料之间的碰撞事件，需要利用厚厚的土层和岩石来屏蔽高能宇宙线的干扰。到目前为止，还没有探测到暗物质的明确信号，实验技术仍在不断发展之中。

实验室建造在地下能够屏蔽背景的同时，如果暗物质与物质相互作用不太弱的话，也有可能屏蔽掉我们想要探测的暗物质粒子，这个问题就可以使用蒙特卡洛模拟方法来研究。暗物质粒子从地表进入到地球内部后的运动可以看作随机游走的过程，我们只要模拟大量粒子的运动轨迹，就能重建出在地下实验室中的暗物质分布。

银河系可能被一个巨大的暗物质晕包围，其中暗物质粒子的速度满足麦克斯韦分布，平均速度大致和银河系中星体的运动速度相当，约为 240 km/s ，称为标准暗晕模型 (Standard Halo Model, SHM)。暗物质粒子在地表处的初始速度将通过这一速度分布抽样得到，随后的随机游走则由两个步骤反复迭代进行：

1. 自由传播：粒子在发生碰撞之前沿直线自由传播一段距离，距离的长度称为自由程。自由程满足特定的概率分布，其平均值即平均自由程，由粒子与地球内部元素相互作用强度和地球的密度确定。模拟中首先计算平均自由程，然后自由程根据相应的概率分布抽样得到。

2. 碰撞：自由运动结束后暗物质粒子与地球内部元素发生碰撞，碰撞将导致暗物质粒子损失一部分能量而减速，并且运动方向改变。新的速度大小和方向由相互作用的具体形式按概率抽样确定，随后重复进行下一段自由传播。

通过这样一步一步的随机过程，可以模拟出暗物质粒子折线形式的运动轨迹，如下图所示[4]：


暗物质粒子轨迹模拟，由坐标原点处垂直向下出发

模拟大量轨迹之后，收集每个粒子经过实验室深度时的速度，就可以统计出地下实验室中的暗物质速度分布，作为直接探测实验信号分析的重要输入条件。下图展示了质量约等于质子质量的暗物质粒子，在地下 2.4 km 深度的速度分布模拟结果：


地下暗物质粒子速度分布



本节介绍的通过随机游走模拟暗物质粒子运动的方法，本质上等同于曼哈顿计划中代号为蒙特卡洛的核裂变反应的中子扩散模拟，此外高能粒子对撞机中的探测器模拟也是采用类似的方法，只是在这些应用中，碰撞过程中会产生的大量次级粒子需要记录和追踪。

结语

蒙特卡洛方法并不特指某种特定的算法，而是对利用随机性来解决问题的方法的统称，我们通过几个例子展示了蒙特卡洛方法的典型应用。物理学，以及所有的科学，都是以实验为基础的，使用蒙特卡洛方法相当于在计算机中进行模拟实验。尽管新物理只可能在真实的实验中发现，模拟实验只能输入已知的物理定律，但正如我们看到的，蒙特卡洛方法可以为理论推导提供佐证，可以用简单的方式解决困难的问题。在开始真实的实验之前，蒙特卡洛模拟也是检验实验方案可行性的重要手段，特别是对于实验成本非常高昂的大科学装置，如暗物质探测器，高能粒子对撞机等等。同样的，在各种工程和应用领域如航空航天工业，军事武器研发等等，蒙特卡洛模拟都是必要的环节。对于我们个人来说，利用一台小小的电脑，就能研究磁性系统的相变，观察暗物质粒子的运动，甚至模拟宇宙的演化、星系的形成这些不可能在现实中直接观测的过程，这本身就是非常有趣的事情。

参考文献

[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

作者简介 /Profile/

夏晨，中国科学院理论物理研究所 19 级博士研究生，研究方向是宇宙线与暗物质直接探测，导师是周宇峰研究员。

本文转载自微信公众号“中国科学院理论物理研究所”。