关于模拟粒子在物质中的输运的一些思考
这篇文章
目录
展开/折叠
前言
高能物理领域和核物理领域经常需要模拟大量微观粒子（如带电轻子、\(\gamma\) 射线、重子、介子等）在物质中的输运过程。模拟这些中高能粒子的输运有助于开展高能物理实验设计、核物理实验设计、评估射线屏蔽效果、以及模拟各种粒子探测器的响应和性能指标等等 [1, 2, 3, 4] 。在放射生物学、核医学领域，模拟高能粒子在生物组织中的输运是评估射线对生物组织的影响的重要一环 [5] 。
大量中高能粒子在物质中输运的模拟方法属于一种蒙特卡洛方法。这一计算方法已被大量的实践充分地证明了其有效性，算法的基本框架已经非常成熟 [6, 7] 。在这种方法中，粒子被当作具有确定的动力学量的经典粒子处理，单个粒子的初态从指定的能谱抽样，动力学量按经典力学演化。特定物理过程的相互作用点以及相互作用末态从相应的截面确定的分布抽样。粒子的存活时间从对应的指数分布抽样，分布参数由宽度决定。这种模拟方法已经应用于成熟的中高能粒子输运计算软件中，如 Geant4 [8, 9] 和 FLUKA [10, 11] 。其中 Geant4 是开源软件，FLUKA 是专有软件。
微观粒子输运的蒙特卡洛方法的基础在于什么？这个问题值得从头思考。在模拟大量粒子的输运时，我们简单地将被模拟的粒子当作一群经典粒子，但微观粒子本该是用概率幅描述的，服从量子力学规律，为什么可以把它当成经典的来计算？这个问题的答案不是显而易见的。除了计算结果符合实验结论这一基础性的保障之外，更仔细地考虑一下这种算法的物理基础是值得的。
要准确模拟粒子在物质中的输运，关键在于正确计算粒子在相互作用点之间的运动、正确选取相互作用点、以及根据相互作用截面确定粒子的末态 [7] 。
粒子运动的处理
与低能粒子相比，中高能粒子与物质的相互作用截面通常很小。在宏观尺度上看，相互作用点离散地分布在粒子的径迹上，在相互作用点之间粒子的运动可视为仅受宏观外场作用的运动。对于微观粒子而言，其状态由各动力学量的振幅描述，粒子状态的演化由对应的运动方程给出。但准确模拟大量微观粒子态的计算开销极其恐怖，在实际应用中没有价值。
因此为了简化问题，需要一种近似手段。一种朴素的想法是把这些微观粒子当作准经典粒子来处理，它们的动力学量是普通的c数。这样的操作明显和量子力学是不相容的，这种近似的合理性在哪里？这个问题涉及到粒子运动的物理基础。我们从动力学的角度来看看这样操作的合理大概在哪。
物理基础
从这种角度来看，在这种方法中，模拟的粒子和现实世界中的粒子是一一对应的关系。被计算的粒子动力学量应当视为实际动力学量的期望值，例如计算中的位矢应当视为粒子位矢的期望 \(\langle\vec{x}\rangle\)，计算中的动量应当视为粒子动量的期望 \(\langle\vec{p}\rangle\)，等等。
要保证这种经典近似操作的合理性，需要粒子动力学量的均值演化服从经典力学，并且粒子动力学量的展宽始终足够小。作为阐释，前者可由Ehrenfest定理说明，后者可从波包宽度的演化规律中略窥一二。
Ehrenfest定理
有质量粒子的Ehrenfest定理写作 [12] \[ \frac{\mathrm{d}\langle\vec{p}\rangle}{\mathrm{d}t}=\left\langle\vec{F}(\vec{x})\right\rangle~. \] 其中 \(\vec{F}(\vec{x})\) 是粒子所处的经典背景诱导的经典力场。例如由电磁场产生的电磁力，由引力场产生的引力之类。从Ehrenfest定理已经可以看出经典运动方程的影子，但它们还略有不同。不同之处在于Ehrenfest定理的右手边是 \(\left\langle\vec{F}(\vec{x})\right\rangle\)，而经典运动方程的右手边是 \(\vec{F}(\langle\vec{x}\rangle)\) 。事实上，对于不太陡峭的背景，即变化足够平缓的 \(\vec{F}(\vec{x})\)，Ehrenfest 定理和经典动力学方程是可以近似等同的。为了看出这点，将 \(\vec{F}(\vec{x})\) 在粒子位置的期望 \(\langle\vec{x}\rangle\) 处展开： \[ \vec{F}(\vec{x})=\vec{F}(\langle\vec{x}\rangle)+\frac{\partial\vec{F}}{\partial x^i}\bigg|_{\langle\vec{x}\rangle}(x^i-\langle x^i\rangle)+\frac{\partial^2\vec{F}}{\partial x^i\partial x^j}\bigg|_{\langle\vec{x}\rangle}(x^i-\langle x^i\rangle)(x^j-\langle x^j\rangle)+\cdots~, \] 并求期望， \[ \left\langle\vec{F}(\vec{x})\right\rangle=\vec{F}(\langle\vec{x}\rangle)+\frac{\partial\vec{F}}{\partial x^i}\bigg|_{\langle\vec{x}\rangle}\left\langle x^i-\langle x^i\rangle\right\rangle+\frac{\partial^2\vec{F}}{\partial x^i\partial x^j}\bigg|_{\langle\vec{x}\rangle}\left\langle(x^i-\langle x^i\rangle)(x^j-\langle x^j\rangle)\right\rangle+\cdots~. \] 注意期望具有的线性性，一阶项为零。二阶项是位矢分布的协方差： \[\left\langle x^i-\langle x^i\rangle\right\rangle = \langle x^i\rangle-\langle x^i\rangle=0~,\] \[\left\langle(x^i-\langle x^i\rangle)(x^j-\langle x^j\rangle)\right\rangle = \left\langle\Delta\vec{x}^2\right\rangle^{ij}~,\] 则 \[ \left\langle\vec{F}(\vec{x})\right\rangle=\vec{F}(\langle\vec{x}\rangle)+\frac{\partial^2\vec{F}}{\partial x^i\partial x^j}\bigg|_{\langle\vec{x}\rangle}\left\langle\Delta\vec{x}^2\right\rangle^{ij}+\cdots~. \] 由此可看出 \(\left\langle\vec{F}(\vec{x})\right\rangle\) 和 \(\vec{F}(\langle\vec{x}\rangle)\) 之间的差别只在于 \(\partial_i\partial_j\vec{F}|_{\langle\vec{x}\rangle}\left\langle\Delta\vec{x}^2\right\rangle^{ij}+\cdots\)。即它们之间的误差主要由 \(\vec{F}(\vec{x})\) 在 \(\langle\vec{x}\rangle\) 处的曲率以及波包的展宽贡献。高能粒子的特性决定了波包的展宽本就非常小，因此只要 \(\vec{F}(\vec{x})\) 足够平缓，平缓到其二阶导数在 \(\sqrt{\det\left\langle\Delta\vec{x}^2\right\rangle}\) 的尺度内可以看作常数，那么 \(\left\langle\vec{F}(\vec{x})\right\rangle\) 就很好地近似于 \(\vec{F}(\langle\vec{x}\rangle)\) 。对于常见的背景场，如重力场、加速电场等，这一条件可以被非常好地满足。因此粒子的平均动力学量自动满足经典动力学方程 \[ \frac{\mathrm{d}\langle\vec{p}\rangle}{\mathrm{d}t}=\vec{F}(\langle\vec{x}\rangle)~. \]
粒子性
但这样还不够，除了期望值按经典力学演化之外，还需要粒子的动力学量分布的展宽足够小，远小于我们关心的精度，从而允许其被当做经典粒子处理。
众所周知 [12, 13, 14]，波包是随时间缓慢扩散的，其扩散速度大致随粒子能量的升高而减缓。因此高能粒子动力学量分布散开的速度本身是相当缓慢的。并且其传播的时间通常也相当短，因此其动力学量的展宽一直是相对小的。
用于阐述，我们考虑一个非相对论性自由高斯波包的扩散。对于一个在笛卡尔坐标系下初始动量为 \(\vec{p}_0\) ，初始展宽为 \(\sqrt{\left\langle\Delta\vec{x}^2\right\rangle}=\mathrm{diag}(\xi,\eta,\zeta)\) ，初始位于原点的高斯波包 \[\begin{equation}\label{gausswavepacket} \psi(\vec{x},0)=\frac{1}{(2\pi)^{3/4}(\xi\eta\zeta)^{1/4}}\exp\left(-\frac{1}{4}\left(\frac{x^2}{\xi^2}+\frac{y^2}{\eta^2}+\frac{z^2}{\zeta^2}\right)\right)\exp\left(i\vec{p}_0\cdot\vec{x}\right)~, \end{equation}\] 其概率幅
附注
附录
非相对论性高斯波包的演化
考虑波包 \(\eqref{gausswavepacket}\) 的自由运动，非相对论性 Green 函数为 \[ iG(\vec{x},t;\vec{x}_0,t_0)=\left(\frac{m}{2\pi i(t-t_0)}\right)^{3/2}\exp\left(i\frac{m(\vec{x}-\vec{x}_0)^2}{2(t-t_0)}\right)~, \] 将其作用到波包上， \[ \begin{aligned} \psi(\vec{x},t)&=\int\mathrm{d}\vec{x}_0iG(\vec{x},t;\vec{x}_0,0)\psi(\vec{x}_0,0)\\ &=\left(\frac{m}{2\pi it}\right)^{3/2}\frac{1}{(2\pi)^{3/4}(\xi\eta\zeta)^{1/4}}\\ &\qquad\int\mathrm{d}\vec{x}_0\exp\left(-\frac{1}{4}\left(\frac{x_0^2}{\xi^2}+\frac{y_0^2}{\eta^2}+\frac{z_0^2}{\zeta^2}\right)+i\left(\frac{m(\vec{x}-\vec{x}_0)^2}{2(t-t_0)}+\vec{p}_0\cdot\vec{x}_0\right)\right)~, \end{aligned} \] 接下来要算这个积分。首先要把它拆开到三个方向，三个方向的积分应当是同一个形式， \[ \begin{aligned} \psi(\vec{x},t)&=(\cdots)\int\mathrm{d}\vec{x}_0\exp\left(-\frac{1}{4}\left(\frac{x_0^2}{\xi^2}+\frac{y_0^2}{\eta^2}+\frac{z_0^2}{\zeta^2}\right)+i\left(\frac{m(\vec{x}-\vec{x}_0)^2}{2(t-t_0)}+\vec{p}_0\cdot\vec{x}_0\right)\right)\\ &=(\cdots)\int\mathrm{d}x_0\exp\left(-\frac{x_0^2}{4\xi^2}+i\left(\frac{m(x-x_0)^2}{2(t-t_0)}+{p_x}_0x_0\right)\right)\\ &\qquad\int\mathrm{d}y_0\exp\left(-\frac{y_0^2}{4\eta^2}+i\left(\frac{m(y-y_0)^2}{2(t-t_0)}+{p_y}_0y_0\right)\right)\\ &\qquad\quad\int\mathrm{d}z_0\exp\left(-\frac{z_0^2}{4\zeta^2}+i\left(\frac{m(z-z_0)^2}{2(t-t_0)}+{p_z}_0z_0\right)\right)\\ &:=(\cdots)I_xI_yI_z~, \end{aligned} \] 只用算一个方向的积分，剩下两个方向换个参数代进去就好。以 \(I_x\) 为例，先将指数整理成关于 \(x_0\) 的二次多项式， \[ \begin{aligned} I_x&=\int\mathrm{d}x_0\exp\left(-\frac{x_0^2}{4\xi^2}+i\left(\frac{m(x-x_0)^2}{2(t-t_0)}+{p_x}_0x_0\right)\right)\\ &=\int\mathrm{d}x_0\exp\left(\right)\\ \end{aligned} \]
参考文献
[1] Ai-Yu Bai et al. (MACE Working Group) Snowmass2021 Whitepaper: Muonium to antimuonium conversion. e-Print: 2203.11406 [hep-ph].
[2] R. Abramishvili et al. (COMET Collaboration). COMET Phase-I Technical Design Report. e-Print: 1812.09018 [physics.ins-det].
[3] L. Bartoszek et al. (Mu2e Collaboration). Mu2e Technical Design Report. e-Print: 1501.05241 [physics.ins-det].
[4] K. Arndt et al. (Mu3e Collaboration). Technical design of the phase I Mu3e experiment. Nucl. Instrum. Meth. A, 1014:165679, 2021. e-Print: 2009.11690 [physics.ins-det].
[5] Ioanna Kyriakou et al. Review of the Geant4-DNA Simulation Toolkit for Radiobiological Applications at the Cellular and DNA Level. Cancers, 14(1):35, 2022.
[6] 马文淦. 计算物理学. 北京:科学出版社, 2005.
[7] Geant4 Collaboration. Geant4 Physics Reference Manual. Geant4 documentation.
[8] S. Agostinelli et al. (Geant4 Collaboration). GEANT4--a simulation toolkit. Nucl. Instrum. Meth. A, 506:250-303, 2003.
[10] Giuseppe Battistoni et al. Overview of the FLUKA code. Annals Nucl. Energy, 82:10-18, 2015.
[11] https://fluka.cern/
[12] 曾谨言. 量子力学教程. 北京:科学出版社, 2008.
[13] 喀兴林. 高等量子力学, 第2版. 北京:高等教育出版社, 2001.
[14] 林琼桂. 高等量子力学讲义. 中山大学物理学院课程讲义, 未出版.