vegan稀释曲线 基因丰度_【2020.12】如何模拟宇宙中暗物质丰度的演化？

之前我已经写了好几篇关于暗物质的文章，都是对暗物质的观测证据及其可能的物理本质进行介绍。今天来写点专业的东西——暗物质丰度演化的模拟，内容包括：①什么叫暗物质的丰度；②玻尔兹曼方程及其数值求解方法；③一些具体的模拟结果。其中数值方法是我做相关研究的过程中总结出来的一套方法，希望可以对想要入门这方面研究的读者有所帮助。

本文的内容比较专业，所以需要读者具有一定的数学和物理基础（微积分、统计物理和一些粒子物理的基本概念）。非专业的读者可以跳过本文第三部分的模拟方法，只需要看第一、二部分的基础介绍以及第四部分的模拟结果即可。此外，不熟悉暗物质基本概念的读者建议先看以下几篇文章或回答：

王清扬：【2020.10】暗物质的基本概念​zhuanlan.zhihu.com

暗物质是由什么组成的？如何把暗能量整合进当今的粒子物理学？​www.zhihu.com

王清扬：至暗之处​zhuanlan.zhihu.com

---

一、暗物质及其丰度

长期以来的理论与实验探究表明，我们所熟悉的物质世界在微观上可以被粒子物理的标准模型描述。标准模型中包含61种基本粒子，这个数量看似不少，但其中的大多数都不能独立或稳定地存在，它们在产生后会迅速衰变至为数不多的几种稳定粒子。因此能参与构成宇宙物质组分的标准模型粒子都是稳定的，包括电子、光子、中微子以及由夸克组合成的质子、中子等^[1]。这些粒子都可以被人们直接探测到，所以被称为可见物质，它们通常以星体、星际物质或背景辐射的形式存在。

大量的天文观测显示，可见物质只占宇宙物质总量的一小部分，宇宙中还有不可见的“暗物质”，它们的总质量达到了可见物质的五倍以上。暗物质不是标准模型中的粒子，它们之所以不可见是因为它们不参与电磁相互作用，这使得其无法产生电磁波信号而被人们观测到。目前人们只能通过暗物质的巨大质量产生的引力效应感受到它们的存在，相关的观测证据包括星系旋转曲线的平坦结果、引力透镜效应、CMB各向异性谱以及宇宙大尺度结构的观测等。其中CMB的观测结果告诉我们暗物质在宇宙诞生初期就已存在，它们和可见物质一样都是从宇宙大爆炸中产生的。虽然至今人们对暗物质的物理本质没有任何了解，但理论家们根据观测得到的暗物质的特性，构造了许多超出标准模型的理论来给出暗物质的候选者。这些候选者大多是粒子模型，包括弱作用大质量粒子（WIMP）、惰性中微子、轴子、超大质量“哥斯拉”粒子等。更多关于暗物质的观测证据及候选者的内容可以参考我之前写的那几篇文章。

对于暗物质，目前人们能够精确测量的参数只有暗物质的密度，尤其是现阶段它在整个宇宙中的平均密度，大约是

，这相当于每立方米1.3个氢原子的质量。相比之下可见物质的平均密度只有

。我们假设暗物质是某种粒子，且这种粒子的质量为

，则可以定义暗物质粒子在整个宇宙中的

数密度

，它的意义为单位体积中暗物质粒子的数量。需要注意的是，

宇宙是在膨胀的，即宇宙的总体积在增大，这样一来暗物质的数密度就会随宇宙尺度

的增大而被稀释，有

。在实际研究中，为方便起见，我们需要再定义一个不被宇宙膨胀稀释的数密度——

丰度

，它的定义是给数密度

乘以一个体积因子：

。可以看到体积因子

会抵消掉宇宙膨胀带来的

的稀释，这相当于定义了一个随宇宙膨胀而增大的体积中粒子的数量（专业术语叫

共动体积中粒子的数密度）。这样一来，我们在研究问题时就不需要考虑宇宙膨胀对暗物质粒子密度的稀释，处理问题就会很方便。

上一段绕来绕去说了这么多，读者可能都被绕糊涂了。到底什么叫数密度？什么叫丰度？我来举个例子吧，这个例子可能对有密集恐惧症的同学不太友好……考虑一个爬满蚂蚁的气球，在一开始气球表面每平方米有100只蚂蚁，这是它此时的数密度。把这个气球吹大，让气球的表面积增大一倍，此时每平方米就只有50只蚂蚁了，这就叫膨胀对数密度的稀释。而丰度是什么呢？我在气球表面画一个圈，圈内有100只蚂蚁，不管这个气球将来是充气还是漏气，圈里面蚂蚁的数量总是不变的，这就叫共动体积中的数密度——丰度。这样一来我们在研究气球表面蚂蚁密度的演化时就不需要再考虑气球膨胀对蚂蚁数密度的稀释，就可以专注于其他因素对数密度的影响。

图1 宇宙的膨胀或收缩会影响到宇宙中单位体积内粒子的数量，即数密度，但不会影响共动体积（红圈）内粒子的数量，即丰度。（当然，前提是粒子不会产生湮灭）

二、玻尔兹曼方程

虽然我们定义的暗物质的丰度不被宇宙膨胀所稀释，但这并不意味着丰度这个物理量不在宇宙膨胀的过程中改变。这是因为在宇宙膨胀的过程中宇宙的温度（粒子平均能量）会降低，而温度的变化会影响到暗物质粒子的产生与湮灭，这样一来暗物质的丰度就会在宇宙膨胀的过程中发生演化。用前面的气球蚂蚁的例子来比喻就是，气球膨胀的过程中虽然蚂蚁的丰度没有稀释，但蚂蚁之间互相远离对方会感到孤独然后抑郁而死，这使得气球膨胀的过程中会有很多蚂蚁死掉，从而导致蚂蚁的丰度发生演化。

那么蚂蚁……啊呸，暗物质的丰度如何随宇宙的温度演化呢？这就是玻尔兹曼方程所研究的事情了。在介绍这个方程之前我需要再做两件事，首先是给大家描述一下早期宇宙中暗物质粒子的产生湮灭过程。通常来说，在微观上我们会假定暗物质的产生与湮灭是通过“2到2过程”发生的：

，其中1 2为两个暗物质粒子，3 4为两个其他粒子。也就是说，我们假定

暗物质粒子由两个其他粒子的碰撞产生，两个暗物质粒子相遇也可以湮灭为两个其他粒子。这样一来，其他粒子较多的时候，它们就会频繁碰撞产生暗物质粒子（如果可以的话），导致暗物质丰度上升。当暗物质粒子较多的时候，它们就会频繁相遇湮灭为两个其他粒子，导致暗物质丰度下降。当然，正向过程与逆向过程也可以达成 热平衡，此时暗物质的丰度不变。

要做的第二件事是引入一个物理量——质温比

，即暗物质的质量与宇宙温度的比值。之所以定义这个物理量是因为，宇宙的温度是通过“操纵”暗物质的产生湮灭来影响暗物质的丰度。而一般情况下温度只有接近暗物质的质量时

^[2]，才可以对暗物质的产生湮灭进行有效的操控。比如质量为100GeV的WIMP型暗物质粒子，宇宙温度是1000GeV还是10000GeV对它来说没有什么区别，只有宇宙温度（粒子平均能量）下降到它质量以下时，它再难以通过其他粒子的碰撞产生出来，它的丰度才会发生显著的变化。因此对暗物质来讲，宇宙温度的绝对数值没有意义，质温比才是有意义的。

现在我们就可以借助玻尔兹曼方程去讨论暗物质丰度的演化了，它是暗物质的丰度

关于质温比

的一个非线性微分方程

^[3]：

方程右边的

是自旋为

的暗物质粒子在某质温比

下处于热平衡时的丰度，具体表达式为：

玻尔兹曼方程中的参数

定义为：

其中

为系统的自由度数量，

为暗物质粒子的质量，

为普朗克质量。

为前述2到2过程的平均散射截面，它是

的函数。

综上所述可以看到玻尔兹曼方程中有两个待定参数，暗物质粒子的自旋

、

质量

。以及一个待定函数——暗物质与其他粒子的

平均散射截面

^[4]，这个需要根据具体的粒子物理模型去算。因此我们模拟暗物质丰度演化要做的就是：①给定一个包含暗物质候选者的粒子物理模型；②写出此模型中暗物质粒子的自旋、质量，并通过量子场论去计算它的产生湮灭过程的散射截面；③把这些代入玻尔兹曼方程并求解得到函数

。

这个函数描述了随着时间演化/宇宙膨胀/温度降低（质温比增大），暗物质的丰度如何发生变化。那么如何求解玻尔兹曼方程呢？很遗憾，这个方程太过复杂没有通解。不过可以证明在设定明确的初始条件的前提下，此方程有唯一的解 ^[5]。因此我们可以通过数值计算的方法去找到这个数值解，从而模拟暗物质丰度的演化规律。

三、数值模拟方法

我习惯用MATLAB写程序，所以采用以下的思路进行数值模拟：
①设置暗物质粒子的自旋

，数值积分得到

的数值表达式。根据积分的渐近性质，使用MATLAB的拟合工具包来将

拟合为一个解析函数。

②根据具体的粒子物理模型手动计算暗物质产生湮灭过程的散射截面

。之后设定一个温度区间，根据粒子在各温度下的动量统计分布，对截面中的动量进行数值积分，得到各温度或质温比下的平均截面作为一个数值函数

，将其拟合为一个解析函数后输入玻尔兹曼方程。

③根据粒子物理模型设定方程的其他参数以及初始条件。
④用MATLAB的ode23s指令解方程。
接下来讨论这些的具体实现。

3.1

的数值积分

• 我们假定暗物质是某种自旋为1的粒子，则：
  
  数值积分的思路是将
  和
  离散化，对于每一个
  将
  数值积分掉（可以用最简单的欧拉法），这样就得到了
  的数值函数。这里有两点需要注意。
• 第一，从积分表达式可以看出
  是随
  的增大而指数压低的，当
  左右时
  的值就会低于
  ，这是双精度浮点运算的极限，我们无法计算比这更小的数。因此如果使用双精度浮点运算，函数
  的定义域就只能是大约
  ，如果想解锁
  的更大的取值就必须要使用更高精度的浮点运算。
• 第二，关于
  的积分范围，显然我们不可能在数值计算中从0一直积分到∞，因此需要对积分上限进行截断。我采用的积分范围是
  ，这个上限可以保证积分的精度也兼顾了程序运行的速度。

3.2

的拟合

• 前面的数值积分给出了
  的数值函数，现在需要将它拟合为一个近似的解析表达式。为了确定拟合函数的形式，我们根据前面的积分表达式来计算
  的渐近性质。
• 首先看
  较小时
  的性质。将其在
  附近做泰勒展开：
  
  可以计算得到奇数阶导数都为0，所以：
  
  其中
  ，
  ，拟合函数需要满足这样的泰勒展开式。另外，拟合函数还要满足
  远离0点时
  会迅速指数压低向一个定值收敛。同时满足以上两条性质的函数可以是
  或者
  ，这就是
  较小时
  的近似表达式。
• 再来看
  较大时
  的性质。我们使用广义二项式定理把前面积分表达式的分母中的
  在
  做展开：
  
  取展开式前两项代入前面的积分表达式得：
  
  由于
  较大，所以可以忽略分母中的
  ，此时积分有解析表达式：
  。这就是
  较大时
  的近似表达式，为了将它与前面的
  较小时的表达式融合，我们还需要给它做一点小小的改动：
  ，这样一来就不会在
  时有
  了。
• 最终，我们将两种情况下的近似表达式融为一体得到：
  这就是
  的拟合函数的形式
  。将这个拟合函数与3.1小节得到的数值函数一同输入MATLAB的拟合工具包中，我们就可以得到四个待定参数的数值，分别为
  ，
  ，
  ，
  。拟合结果如下图所示，可以看到拟合效果是非常好的。最后要提醒的是，
  这四个参数的数值是针对自旋1粒子拟合的，对于自旋0和自旋1/2的粒子需要重新进行3.1中的数值积分，然后对新的数值函数进行拟合。

图2 自旋1粒子的Yeq(x)，蓝线为数值积分得到的数值函数，红线为拟合函数，可以看到二者几乎完全重合，拟合效果非常好。

3.3 假设模型 计算散射截面

现在我们假设一个包含自旋1的暗物质的粒子物理模型去计算散射截面（之所以用“假设”这个词，是因为至今我们对暗物质的物理本质一无所知……）目前人们已构建出了许多看起来很合理的暗物质的粒子物理模型，我的这篇回答的第二部分对它们进行了总结，可以点这里回顾一下：暗物质是由什么组成的？如何把暗能量整合进当今的粒子物理学？当然，这些模型都是相当复杂的模型，本文的讨论不需要那么复杂，假设一个toy model就够了，比如这样：

。看过我粒子物理专栏1.3节（1.3 强相互作用与弱相互作用）的朋友应该能看出这个公式是什么意思。式中

为某个自旋0的标量粒子，

为自旋1的暗物质粒子，

这个公式表示两个

粒子与两个

粒子直接耦合在一起发生相互作用

，

的大小表示它们耦合的强度。也就是说，在这个模型当中最简单且最普遍的相互作用过程是：

暗物质粒子由两个标量粒子碰撞产生，以及其逆向过程——两个暗物质粒子相遇湮灭为两个标量粒子。

让我们计算上述产生湮灭过程的散射截面，如下图3所示，结果为：

其中

为暗物质粒子的动量的模，

是它的质量。

是两粒子的总能量，由粒子的动量和质量唯一决定。

图3 两标量粒子Φ与两矢量粒子W直接耦合的树图阶散射截面计算。

在宇宙中，大量粒子的动量是各不相同的，因此大量粒子之间两两散射的截面也各不相同。但统计物理告诉我们大量粒子的运动是有统计规律，其动量满足费米分布或玻色分布。所以我们可以将散射截面对所有粒子的动量分布积分来得到平均散射截面

，公式如下

^[6]：

显然，

是温度

或者说质温比

的函数。通过数值积分，我们可以设定一个离散的温度区间来计算得到这个

关于

的数值函数。最后，将数值函数拟合为一个解析函数，便可代入玻尔兹曼方程开始求解工作。对我们这个模型，拟合函数的形式可以设为：

其中

，

，

，

的数值与质量的单位有关，如果质量单位是GeV，则

。

3.4 玻尔兹曼方程的数值求解

到此为止我们已将玻尔兹曼方程中的

与

处理好，接下来进行数值模拟的最后也是最重要的一步——解方程！我直接用了MATLAB的ode23s指令，把

与

直接体现在fun函数中了，因此代码很简单，只有几行。运行结果就是玻尔兹曼方程在我们所设求解范围内的数值解

，即我们的模型中暗物质的丰度随宇宙膨胀/质温比增大的演化规律。我把代码放在下面，感兴趣的读者可以拿去运行一下。

%单位GeV
gs=100;Mp=1.22e19;mw=;%暗物质粒子的质量，需要设定
c1=0.1738;c2=1.495;c3=-4.840;c4=1.594;
k1=1.403;k2=-0.5993;k3=-3.111;k4=0.8555;

i=;g=10^i;%模型中的耦合常数，需要设定
fun=@(x,Y) -(sqrt(45/(4*pi^3*gs))*mw*Mp*(exp((k1*tanh(log10(x)+k2)+k3)*log10(x)+k4)*g^2/mw^2)) ...
*(Y^2-((c1+exp(-c2*x^2+c3))*exp(-x)*(x+c4)^1.5)^2)/x^2;%玻尔兹曼方程
x_range=[0.01 1000];%求解范围
%对热产生粒子求解范围建议设为[0.01 1000]。对非热产生粒子建议设为[Tp 1000],Tp为粒子产生时的温度。
Y0=; %初始条件，需要设定（热产生粒子设定为Y0=Yeq(0)，非热产生粒子设为Y0=0）
options=odeset('InitialStep',1e-6,'MaxStep',0.1);%求解精度
[x,Y]=ode23s(fun,x_range,Y0,options);%解方程
loglog(x,Y,'LineWidth',1.5);hold on %把求解结果画图

四、模拟结果

第三节介绍了具体的模拟方法，这一节我们来实践一下。3.4小节的模拟程序中只有三个地方需要设定：①暗物质粒子的质量；②相互作用的耦合常数；③求解范围与初始条件（由粒子的产生方式决定）。现在分别考虑我们模型中的

粒子作为两种不同的暗物质候选者：质量在100GeV左右的弱作用大质量暗物质粒子WIMP，以及质量在

GeV左右的超大质量“哥斯拉”粒子。我们将以它们两个为例给出宇宙中暗物质丰度演化的模拟结果。

4.1 WIMP型暗物质

先来说WIMP，WIMP的产生方式通常是“热产生”，即产生于温度远高于其质量能标的极早期宇宙当中（比如在温度

的宇宙重加热过程中产生）。在产生过程中WIMP的丰度会急剧增加，之后它会很快与产生它的粒子达到热平衡（比如我们模型中的

粒子），即其产生过程与湮灭过程相当。这使得

高温时WIMP的丰度保持在一个定值，这个定值就是前面我们说的

在

时的数值。因此，

模拟WIMP的演化时我们采用的初始条件是

。将这一初始条件输入程序，同时将WIMP质量设为

。为了让读者看到耦合常数

如何影响暗物质丰度的演化，我们将

设定为0.1 0.03 0.01 0.003 0.001这五个数值，依次进行模拟，结果如图4所示。

• 可以看到在
  ，即宇宙温度高于暗物质质量能标100GeV时，暗物质的丰度就如前面所说——保持在
  基本不变。
• 之后随着宇宙的膨胀，宇宙的温度在不断降低。当温度降低至100GeV以下，即
  时，暗物质的产生过程
  由于能量不够而再难以发生。此时
  产生过程与湮灭过程的平衡被打破，暗物质的湮灭过程开始成为主导，使得暗物质的丰度急剧下降，图上可以清楚地看到这一过程。
• 在宇宙温度降低至1GeV左右时，即
  ，暗物质的丰度已下降了上亿倍。此时宇宙中的暗物质已足够稀薄，各个暗物质粒子之间相距甚远以至于再难以碰到一起进行湮灭。也就是说，
  宇宙温度足够低时暗物质粒子之间的湮灭过程也基本不再发生。这意味着，随着宇宙的进一步膨胀，暗物质的丰度基本“冻结”在一个定值，这个定值被称为“剩余丰度”，即为我们通过天文观测确定的今天宇宙中的暗物质的丰度。
• 还可以看到，耦合常数
  越大，暗物质的剩余丰度越小。这一点很好理解，因为耦合常数表征的是相互作用的强度，相互作用越强，演化中湮灭的暗物质粒子就越多，剩余下来的自然就少了。

图4 我们的模型中质量为100GeV的WIMP型暗物质的丰度随宇宙温度演化的模拟结果。

目前人们通过天文观测获得的今天宇宙中的暗物质丰度大约是

，对于

的WIMP，有

。因此，

如果暗物质的物理本质被我们的模型所描述的话，那么从图4可以看到，目前的观测数据会限制我们模型中的耦合常数

。这是一个非常不平凡的结论，因为0.03这个数值除以

后

^[7]，恰好与粒子物理标准模型中弱作用的耦合常数在一个数量级。这意味着WIMP型暗物质相关的物理模型与已知的标准模型相比并不显得突兀，二者很可能有非常大的联系，这就是我们认为WIMP是非常合适的暗物质候选者的原因。

4.2 “哥斯拉”型暗物质

现在考虑我们模型中的

粒子作为质量在

GeV左右的“哥斯拉”型暗物质。由于哥斯拉粒子的质量非常大，所以它产生时宇宙的温度一般约等于或小于它的质量能标，这使得哥斯拉粒子从未与其他粒子处于热平衡，因此这种产生方式被称为“

非热产生”。对于非热产生的粒子，在求解玻尔兹曼方程时需要设定一个粒子产生时的质温比

，让方程从

开始求解。且初始条件应当设为

，即一开始宇宙中没有暗物质粒子，它们在宇宙温度下降至

时才产生。

现在开始对哥斯拉型暗物质的丰度演化进行模拟。将上述的初始条件与

输入程序，由于非热产生中耦合常数不是影响剩余丰度的主要因素，我们将耦合常数设为一个定值

，设置多个

来看哥斯拉粒子产生时的温度对它丰度演化的影响。模拟结果如图5所示。

• 可以看到在宇宙温度高于我们设定的阈值时，宇宙中是没有暗物质的。当温度降至阈值后，作为暗物质的哥斯拉粒子会迅速产生，并且其丰度会迅速冻结为一个定值。这个冻结丰度就是今天我们观测到的宇宙中暗物质的丰度。
• 同时也可以看到，哥斯拉粒子产生时的温度对它的冻结丰度影响巨大，温度稍一偏离，冻结丰度就会偏离好几个数量级。因此，如果认为暗物质是哥斯拉粒子的话，理论上就会存在一个精细调节的问题。也就是说，我们观测到的今天宇宙中的暗物质总量和可见物质总量相差不到一个数量级，如果认为这件事不是巧合的话，那么就需要将暗物质产生时宇宙的温度调节得足够精细（在我们的设定下，大约
  ）才能将暗物质与可见物质总量控制在一个数量级。这样一来，理论就显得非常不自然。不过一些研究表明这个精细调节问题与理论的可重整性有关，是可以被解决的，具体可以参考arXiv: 0911.1120这篇文章。

图5 我们的模型中质量为10^15GeV的“哥斯拉”型暗物质的丰度随宇宙温度演化的模拟结果，其中耦合常数设为g=0.03。

---

最后进行一个总结：

1. 粒子“丰度”的定义是共动体积中粒子的数密度，这个物理量不会随着宇宙膨胀而稀释。
2. 宇宙的持续膨胀会导致宇宙温度持续降低，从而影响宇宙中暗物质粒子的产生与湮灭。因此暗物质的丰度会随宇宙温度的降低而发生变化，这一过程被玻尔兹曼方程描述。
3. 在设定一个具体的暗物质模型之后，我们可以计算并拟合玻尔兹曼方程中的待定参数/函数，从而使用数值方法去求解玻尔兹曼方程，模拟暗物质的丰度随宇宙温度的演化规律。
4. WIMP型暗物质的产生方式是“热产生”。它会在宇宙温度远高于其质量能标时产生，并很快与其他粒子处于热平衡，丰度保持不变，直到宇宙温度降低至其质量能标以下，WIMP的产生通道被“切断”，热平衡才被打破。此时WIMP的湮灭过程开始占主导，它的丰度急剧下降。当丰度下降至足够小，WIMP足够稀薄以至于难以相遇发生湮灭，它的丰度就会“冻结”在一个定值，这个值就是今天我们观测到的宇宙中暗物质的丰度。
5. 超大质量的“哥斯拉”型暗物质的产生方式是“非热产生”。通常它产生时宇宙的温度比它的质量能标低一些，因此存在一个温度阈值使得当宇宙的温度降低至这个阈值时哥斯拉粒子会迅速产生，丰度急剧上升。不过很快产生过程会停止，哥斯拉粒子的丰度会如同WIMP一样“冻结”在一个定值，即今天观测到的暗物质的丰度。

参考

1. ^注：中子只有处于原子核内或中子星内才是稳定的，自由中子的平均寿命只有896秒。
2. ^注：在自然单位制下，质量、能量、温度的量纲是一样的，可以进行比较。
3. ^玻尔兹曼方程的推导参考杨炳麟《暗物质及相关宇宙学》第11章。点击此链接可以跳到此书英文版下载地址。 https://link.springer.com/article/10.1007/s11467-016-0583-4
4. ^注：这里其实应该是平均散射截面与粒子间相对速度之积，但是相对速度可以与截面中的相对速度因子抵消掉，所以这里只称平均截面。
5. ^解的唯一性的证明见杨炳麟《暗物质及相关宇宙学》第11.3节。
6. ^这个公式可由杨炳麟《暗物质及相关宇宙学》第11.2节(11.2.9)式推导得出。
7. ^注：之所以要除以质量的平方，是因为我们模型中的耦合常数是无量纲的，而弱作用耦合常数有[M]^-2的量纲，要比较的话需要将量纲统一。