撰文 | 夏晨

什么是蒙特卡洛

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

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

单位圆的面积

1000 个随机点

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

统计物理中的伊辛模型

一维格点伊辛模型的严格解在 1925 年已由伊辛本人求出，二维无外磁场伊辛模型由昂萨格于 1944 年求解，杨振宁先生在 1952 年发表了二维伊辛模型磁化强度的推导，而三维伊辛模型的精确解到目前为止还没有得到。二维伊辛模型的精确解表明系统存在铁磁相变，对应的相变临界温度为

上面两个步骤反复迭代就能得到一系列的系统构型。可以看到如果新状态的能量更低，则一定会被接受，系统可以很快地向能量更低的状态移动，同时也有机会跳转到能量更高的状态，最后得到的样本可以证明满足玻尔兹曼分布。此外 MCMC 方法不需确切地知道概率分布的归一化系数，从而避免了配分函数的计算。

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

二维伊辛模型 MCMC 模拟

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

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

暗物质在地球内部的运动

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

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

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

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

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

自由传播：粒子在发生碰撞之前沿直线自由传播一段距离，距离的长度称为自由程。自由程满足特定的概率分布，其平均值即平均自由程，由粒子与地球内部元素相互作用强度和地球的密度确定。模拟中首先计算平均自由程，然后自由程根据相应的概率分布抽样得到。碰撞：自由运动结束后暗物质粒子与地球内部元素发生碰撞，碰撞将导致暗物质粒子损失一部分能量而减速，并且运动方向改变。新的速度大小和方向由相互作用的具体形式按概率抽样确定，随后重复进行下一段自由传播。

通过这样一步一步的随机过程，可以模拟出暗物质粒子折线形式的运动轨迹，如下图所示[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级博士研究生，研究方向是宇宙线与暗物质直接探测，导师是周宇峰研究员。

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

特 别 提 示

1. 进入『返朴』微信公众号底部菜单“精品专栏“，可查阅不同主题系列科普文章。

2. 『返朴』提供按月检索文章功能。关注公众号，回复四位数组成的年份+月份，如“1903”，可获取2019年3月的文章索引，以此类推。