刘 岩，刘献龙，罗志环，李 海，杨 意

（华南农业大学 理学院应用物理系，广东 广州 510642）

格点规范理论课程中的数值模拟方法的应用

刘 岩，刘献龙，罗志环，李 海，杨 意

（华南农业大学 理学院应用物理系，广东 广州 510642）

本文在格点规范理论课程中，采用数值模拟中的Metropolis方法，组织学生解决了一个可以使用理论方法计算的问题，即规范场作用量，并把数值结果与理论结果做了比较，加深了学生对离散化格点这一模型的理解，并进而熟悉和掌握在离散化的格点下研究问题的原理和方法，取得了很好的学习效果.

格点规范理论；数值方法；教学改革

自然界存在的四种基本相互作用——强相互作用，弱相互作用，电磁相互作用以及引力相互作用，它们都有各自的理论体系进行描述，其中对于描述强相互作用的可重整化规范理论，即量子色动力学（QCD），是在上个世纪70年代才被提出.不过，QCD是微扰的理论，仅在耦合常数较小，或者说强相互作用高能的部分有效，而在强相互作用的低能部分，微扰的QCD失效.为了研究强相互作用的低能非微扰过程，1974年K.Wilson首先提出了格点规范理论（LGT），其基本思想是将连续的时空用离散的晶格代替，这样，连续理论无法计算（发散的）物理量就可在格点上计算，最后还要用重整化群方法把结果外推到连续极限.如图1所示，其中格点间的距离即格距为a，费米场定义在格点上，而规范场则定义在两个格点间的链上.

图1 格点规范理论的基本元素.其中费米场和规范场分别定义在a）,b）处

在格点规范理论的教学中，会教导学生首先掌握规范场和费米场的作用量，而格点下的规范场和费米场作用量是由连续理论下的作用量在格点下离散化后得到的，为了让学生理解连续理论离散化到格点的，以及格点理论下的结果在格距a→0时可以回归到连续理论的过程，一般在课堂上会给出以上过程的理论推导，但是，为了帮助学生进一步的理解连续性和离散化的相互转化，并进而熟悉和掌握在离散化的格点下研究各种问题的原理和方法，可以在教导理论推导的同时，组织学生采用数值计算的方法研究格点问题.下面以格点下的规范场作用量为例，让学生使用Metropolis方法研究1+1维U（1）格点规范场作用量与β的函数关系，同时计算对应平均作用量的理论值，将理论值与数值结果进行比较，并作出相应的对比图形.

由Wilson提出的方块算符（plaquette）构成的规范场作用量如下：

其中Pμv代表格点的方块变量，它由四条规范链的有序乘积所组成，如图2所示：

图2 方块算符Pμv

下面组织学生用Metropolis方法数值计算该规范场作用量.首先对Metropolis方法进行说明：

（1）通常情况下，组态U是被随机定义在格点上的，为了求解相关的物理量，我们先引入物理量对应算符F的期望值[1]：

关于应用到格点上的组态，当组态数目较大时，我们可以以格点离散的形式得到近似的表达式[2]：

上式可以看到，要求得相应的表达式必须对所有组态进行求和.不管怎么说，这是不太现实的，在统计物理方面，人们选择“重要抽样”的思想，只对贡献较大的组态进行求和，对于贡献较少的部分则省略掉.这样，物理量可以通过求选取的组态平均值得到.

（2）组态的选择取决于玻尔兹曼因子（Boltzmann factor），即e-△S（U），Metropolis方法中通过Update实现组态的选择：

假设更新前的组态和更新后的组态分别为U和U'，二者满足：

其中，△S（U）是组态从U变化到U'所产生的作用量的变化量，并且（4）中的玻尔兹曼因子 e-△S（U）满足[3]：

组态更新与否并非完全确定的，组态U以一定的概率确定是否更新为U'，用P（U→U'）表示这一概率，则P（U→U'）满足：

通常，由于组态的初始化是通过随机数实现的，所以组态的更新需要进行一定的时间才能达到合理的状态，即组态的分布需要一定次数的更新之后才能满足重要抽样条件.Metropolis方法的关键也在于找到合适的系统组态，之后再在该系统组态下，计算相关的物理量.

我们使用以上介绍的Metropolis方法计算规范场作用量，该模拟过程主要分为以下四个步骤：

1 确定格点的初组态

为了计算格点的方块平均值，我们需要确定格点最初的组态，包括格点的大小和格点上的规范场作用量.格点的大小可以根据实际计算确定，通常为了减少运算量，可以用4×4的格点验证程序的可行性，进而可减小格距，增大格子数目，如使用16×16、32×32的格点等.

程序模拟中通过函数rand（）产生一组0～1之间的随机数来确定各个格点上的规范场作用量U，具体的表达式满足：

2 格点组态的更新

格点组态的更新是为了得到一个相对平衡的格点系统，主要通过Metropolis方法条件来进行更新组态：更新前的组态U与另外3个组态U1U2U3共同组成一个方块UP，即更新前有：

之后用另一组随机数产生一个新的组态U'，可以组成一个新的方块：

将（8）和（9）分别代入（1），可以计算得到两个方块的作用量，其中β的值已知.

然后，根据Metropolis方法原理，比较两个方块作用量的变化：

（1）如果新的组态产生的方块作用量比原组态构成的方块作用量小，则用新的组态U'代替原来的组态U；

（2）如果新的组态产生的方块作用量比原组态构成的方块作用量大，则可能用新组态U'代替原组态U，也可能保持原组态U不变，其更新概率按照（6）式进行.

以上即完成了一个组态的更新，运用同样的方法对格点系统中的每一组态进行更新，可完成一次整个格点系统的组态更新.通常情况下，格点的更新需要反复进行，才能得到尽可能准确的规范场分布，从而进行作用量计算.

3 格点作用量平均值的计算

（1）格点系统的预热过程：在计算作用量平均值之前，先根据步骤2对格点系统更新200次；

（2）计算一个格点系统的总作用量：完成200次格点系统更新后，在最新的规范场分布下利用（1）式计算每个方块的作用量，求和得到一个总的作用量；

（3）持续计算一组格点系统的总作用量：计算出一个格点系统总的作用量后，按步骤2继续更新格点组态，每完成5次格点系统的更新计算一个系统总的作用量，总共计算500个格点系统总的作用量；

（4）计算一个格点方块的作用量平均值：对所求的500组格点系统总作用量取算术平均可得到格点系统的总作用量平均值，再除以格点系统的格点方块数量，即可得到单个格点方块的作用量平均值.

4 均匀改变公式（1）式中β的值，重复步骤1-3

程序中β的取值范围为[0,6]，为了得到有效的拟合曲线，β的取值间隔越小越好，这里以0.2为取值间隔，得到多组β值下对应的格点方块作用量平均值.

5 求作用量的理论值

为了验证格点QCD方法的正确性，我们在得到格点QCD作用量的平均值后，需要计算出它的理论值，从而能进行结果比较.

作用量当中的Re（UP）的理论值为：

其中，I0（β）和 I1（β）满足[4]：

由此得到：

由此，我们分别得到了规范场作用量的数值计算和理论计算结果，列于表1：

表1 不同β下的方块作用量理论计算以及数值计算结果

表中的理论计算结果可通过第一类的虚宗量贝塞尔函数得到，实际操作中，可以通过Matlab、Excel等工具计算得到.

以β为横坐标、方块作用量为纵坐标，可以得到如图3所示的对比曲线：

图3 方块作用量数值计算与理论计算结果的对比曲线

从图形可以看出，计算机模拟的结果与理论值达到了有效的吻合，二者可分辨的误差极小.这正好验证了格点方法在非微扰QCD上运用的正确性.将数据结果进行曲线拟合后，可以近似给出格点QCD经过连续极限外推下的真实值.

通过这一问题的解决，可以看到学生能够使用数值计算方法比较容易的得到与理论解一致的结果，从而熟悉和掌握在离散化的格点下研究问题的原理和方法，并进一步理解连续理论离散化到格点，以及格点理论下的结果在格距a→0时可以回归到连续理论的过程.因此，在格点规范理论课程中，基于离散化这一特殊背景，对一些问题由浅入深的逐步引入数值计算，是可以有效的增加学生的基本功以及对格点理论的理解和掌握，同时也培养了学生数值计算的能力，对其今后进一步的科研工作也有着积极的作用.

〔1〕（美）黄卓然.高能重离子碰撞导论[M].哈尔滨:哈尔滨工业大学出版社,2002.

〔2〕罗向前,刘岩.格点QCD对混杂态介子基态及第一激发态质量谱的研究[D].广州:中山大学,2006.

〔3〕Kerson H.Quarks,leptons and gauge fields[M].2nd ed.Beijing:World Scientific Publishing,1992.

〔4〕吴崇试.数学物理方法[M].北京:北京大学出版社,2003.

G420

A

1673-260X（2012）09-0211-03

华南农业大学2010年度教育教学改革与研究自筹项目（JG10108）；华南农业大学2012年度教育教学改革与研究重点项目（JG12008）