下载:
-
强相互作用是自然界存在的四种基本相互作用之一，其禁闭性质至今仍是未解的世界难题。夸克和胶子通过强相互作用形成色单态的强子，由于夸克禁闭，人们只能观测到强子而无法观测到单独的夸克。宇宙中99%以上的可见物质都由核子构成，由于强相互作用的复杂性，人们对核子内部结构所知甚少，核子的一些基本性质，例如它的质量和自旋的来源，也还不清楚。电子离子对撞机(EIC)是研究核子结构的最佳实验手段，目前国内外都在大力推动EIC的建设，美国已于2020年批准EIC实验[1-2]，中国电子离子对撞机(EicC)也正在推进[3-4]。从理论上，尤其是从第一性原理来研究核子结构，对理解强项互作用和夸克禁闭具有重要意义，并且将为实验研究提供重要理论支撑。
描述强相互作用的基本理论是量子色动力学(Quantum Chromodynamics, QCD)，由于强相互作用在禁闭区域是非微扰的，必须用非微扰的方法来研究，格点QCD是从第一性原理来研究非微扰强相互作用最重要的方法之一[5]。它利用大规模数值计算来对QCD理论进行数值求解。随着近年来计算机技术的高速发展和相关算法和理论的不断突破，格点QCD在强相互作用研究中发挥着越来越重要的作用。格点QCD基本思想是将QCD理论定义在有限大小的离散的欧氏时空中(如图1所示)，利用场论的路径积分形式，数值地计算强相互作用相关物理量。有限的格距和体积分别为量子场论中的紫外和红外发散提供了自然的截断。通常我们需要在不同的格距和体积下进行计算，最后取连续极限和无限体积极限，取极限的过程就是量子场论中的重整化过程。
图 1 (在线彩图)格点QCD基本方法示意
在量子场论中，一个体系的所有物理都包含在路径积分中。离散的四维欧氏空间中QCD理论的路径积分为
$$ Z = \displaystyle\int DU D\bar{\varPsi} D\varPsi \exp\big(-S[U, \bar{\varPsi}, \varPsi]\big) {\text{，}} $$ (1) 其中
$ U $ ，$ \varPsi $ 和$ \bar{\varPsi} $ 分别是规范场、夸克场和反夸克场。$ S[U, \bar{\varPsi}, \varPsi] $ 是欧氏空间中离散化的QCD作用量。任意一个物理量$ \mathcal{O} $ 的值可由下式得到：$$ \langle \mathcal{O} \rangle = \frac{1}{Z}\int DU D\bar{\varPsi} D\varPsi \mathcal{O} \exp\big(-S[U, \bar{\varPsi}, \varPsi]\big) {\text{。}} $$ (2) 为了计算这个物理量，我们首先用蒙特卡罗方法产生一组组态，并使这些组态按照密度分布函数
$ \exp\big({-S[U,\varPsi, \varPsi]}\big)/Z $ 来分布。这样一来，$ \langle \mathcal{O} \rangle $ 的值可以由对各个组态求平均来得到。产生组态是格点QCD中计算量最大的一部分，也是整个格点计算的基础。国际上，组态通常由一些大的格点合作组来完成。中国格点组还没有自己的组态，目前这项工作刚刚展开，本文的研究目的就是利用开源软件Chroma[6]产生一组组态，格点上的夸克作用量采用Wilson Clover作用量，格距约为 0.1 fm，格子大小为
$ 32^3 \times 64 $ ，$ {\rm{{\text{π}}}}$ 介子质量约为 250 MeV。具有这样的参数的组态可以用来研究核子结构和强子谱等物理问题并能得到有意义的物理结果。下面我们将在第2节介绍产生组态的算法，即杂化蒙特卡罗算法(Hybrid Monte Carlo, HMC)；在第3节中介绍具体产生组态的过程；在第4节中展示用产生的组态计算的
$ {\rm{{\text{π}}}}$ 介子和K介子质量；最后在第5节给出总结。 -
HMC是从杂化算法和分子动力学算法所衍生出来的一种普遍的产生玻尔兹曼分布的方法。
假设我们需要产生的概率分布为
$$ P[\phi]\propto {\rm e}^{-S[\phi]} {\text{，}} $$ (3) 这里
$ S\left[\phi\right] $ 表示体系的作用量，$ \phi $ 表示体系的自由度。自然界中空气的分子运动和碰撞形成了一种天然的玻尔兹曼分布。 为了模拟这个过程，我们把$ \phi_x $ 作为所研究体系的广义坐标，这样$ S\left[\phi\right] $ 就成了仅依赖于坐标的势能，再引入与$ \phi_x $ 共轭的广义动量$ {\text{π}}_x $ 和体系的动能。于是体系的哈密顿量可以定义为$$ \mathcal{H}\left[{\text{π}},\phi\right] = \sum\limits_{x}\frac{{\text{π}}_x^2}{2}+S\left[\phi\right] {\text{。}} $$ (4) 我们要求体系按照经典的哈密顿正则方程演化：
$$ \begin{split}& \dot{{\text{π}}}_{x} \equiv \frac{\mathrm{d}{\text{π}}_{x}}{\mathrm{d}\tau} = -\frac{\partial \mathcal{H}[{\text{π}},\phi]}{\partial \phi_x} = -\frac{\partial S[\phi]}{\partial \phi_x} {\text{，}}\\& \dot{\phi}_{x} \equiv\frac{\mathrm{d} \phi_x}{\mathrm{d}\tau} = \frac{\partial \mathcal{H}[{\text{π}},\phi]}{\partial {\text{π}}_x} = {\text{π}}_{x} {\text{，}} \end{split} $$ (5) 这里
$ \tau $ 是人为引入的一个额外的虚拟时间，称为蒙特卡罗时间。从式(5)可以看出哈密顿量实际上是动能部分和势能部分之和，因此体系关于坐标和动量的概率分布是完全独立的，也就是说，只要产生正则分布
$$ {\rm e}^{-\mathcal{H}} = \exp\left(-\sum\limits_{x}\frac{{\text{π}}_x^2}{2}\right)\exp\big(-S\left[\phi\right]\big) {\text{，}} $$ (6) 那么它的势能部分就是我们所要的玻尔兹曼分布，它的动能部分则是可以直接产生的高斯分布。
在详细介绍整个HMC算法过程之前，先简要介绍这个算法过程中常用到的蛙跳积分方法。 蛙跳积分方法主要有以下三个步骤，首先利用初始的场变量
$ \phi(0) $ 将动量积分半步，然后可以利用刚刚更新到半步的动量$ {\text{π}}\left(\delta\tau/2\right) $ 将场变量$ \phi $ 更新一步，最后再利用已经更新的场变量将动量更新半步。具体过程如下$$ \begin{split}& {\text{π}}_{x}\Big(\frac{\delta\tau}{2}\Big) = {\text{π}}_{x}(0)-\frac{\partial S[\phi(0)]}{\partial \phi_x}\cdot \frac{\delta\tau}{2} {\text{，}}\\& \phi_x(\delta\tau) = \phi_x(0)+{\text{π}}_{x}\Big(\frac{\delta\tau}{2}\Big)\cdot \delta\tau {\text{，}}\\ & {\text{π}}_{x}(\delta\tau) = {\text{π}}_{x}\Big(\frac{\delta\tau}{2}\Big)-\frac{\partial S[\phi(\delta\tau)]}{\partial \phi_x}\cdot \frac{\delta\tau}{2} {\text{，}} \end{split} $$ (7) 这里的
$ \delta\tau $ 是我们给定的一个步长，在实际的计算中，上述过程会重复进行$ N_{\rm step} $ 次，下面给出的是HMC算法的具体过程。假设初始的组态为
$ {\phi_x(0)} $ ，每个径迹的长度为$ \tau_0 \!=\! N_{\rm step}\delta\tau $ ，其中$ \delta\tau $ 为步长。(1) 计算每一次哈密顿量的改变量$ \varDelta\mathcal{H} = \mathcal{H}(\tau_0)-\mathcal{H}(0) $ ；(2)如果$ \varDelta\mathcal{H} < 0 $ ，接收新的组态；(3) 如果$\varDelta\mathcal{H}\geqslant 0$ (1)，产生一个在区间(0，1)内的均匀分布的随机数r；(4) 如果$ r < {\rm e}^{-\varDelta\mathcal{H}} $ ，接收新的组态；(5) 否则拒绝接收，重复第一步。HMC算法最主要的一个特点是对于任意
$ \delta\tau $ 都满足细致平衡条件，也就是说，原则上只要迭代次数足够多，最终一定可以得到所期望的玻尔兹曼分布。HMC算法的一大优点是帮助我们完全规避了对不同的$ \delta\tau $ 进行模拟然后再外推的繁琐过程，只需要在开始前选择一个合适的$ \delta\tau $ 就可以了。蛙跳积分会导致$ \varDelta\mathcal{H}\propto(\delta\tau)^2 $ ，因此接受率会当$ \delta\tau $ 过大时急剧下降。接受率虽然不会影响算法的严格性，但是过低的接受率会使在有限的计算时间内获得的等效的统计独立的组态数目大幅减少，进而会放大计算过程中的统计误差。HMC算法的另一个优点是其接受率几乎仅仅依赖于$ \delta\tau $ ，与$ \tau_0 $ 的依赖并不敏感。根据经验，实际计算时一般选择使得接受率在$ 0.5\sim0.9 $ 左右的$ \delta\tau $ 为宜。 -
图2描述的是由美国USQCD合作组织研发的Chroma软件套件[7]，该软件系统设计高度便携化，可以在多种不同的环境里搭建使用。为了实现代码的执行性和可移植性，不同的软件模块层被创造出来，最上面一层是应用层，下面三层是支撑库层，分别实现消息通讯、数据结构定义和基本算法实现的功能。
图 2 (在线彩图)软件模块结构
格点QCD的模拟和测量等计算应用必须建立在不同的软件模块之上。第一层包括格点QCD消息传递应用程序接口(QMP)和格点QCD线性代数应用程序接口(QLA)。QMP实现并行机器上处理器之间的通信，而QLA提供了一个线性代数例程的接口。第二层包括格点QCD输入输出应用程序接口(QIO)和格点QCD数据并行应用程序接口(QDP和QDP++)。 QIO使用QMP在单节点系统和并行机器上实现格点数据的输入和输出。QDP和QDP++的主要特点是可以处理不同的数据类型，因此其可以根据计算需要创建任意大小和维度的格子。此外，QMP、QLA和QIO模块可以在构建应用程序接口时使用。第三层包括几个优化包，QUDA是基于CUDA的格点QCD计算库[8]，是在GPU上的软件包，实际计算过程中主要靠它来计算庞大的费米子矩阵的逆，也是计算中最耗时间的地方，它可以和应用层的Chroma集成。
-
我们使用的规范场作用量为Symanzik改进作用量，并使用了tree-level的tadpole改进，这种作用量有限格距的误差是
$ O(a^2) $ 量级的。夸克作用量使用的是Wilson clover作用量，这种作用量在Wilson作用量上加入Sheikholeslami-Wohlert项[9]:$$ S_{\rm SW} = c_{\rm SW}\frac{\rm i}{4}\sum\limits_{x,\mu\nu}{\bar{\psi}\left(x\right)\sigma_{\mu\nu}\mathcal{F}_{\mu\nu}\left(x\right)\psi\left(x\right)} {\text{，}} $$ (8) 其中
$ c_{\rm SW} $ 是一个依赖于耦合常数的参数，通过调节这个参数获得物理量的改进。$ \mathcal{F}_{\mu\nu} $ 是格点上规范场场强张量。如果考虑tree-level的改进，那么$ c_{\rm SW}\! =\! 1 $ 。在tadpole改进下，最终我们在计算中使用的Wilson clover作用量可写为$$ \begin{split} S_f^{\left(TI\right)} =& \sum\limits_{xy}{{\bar{\psi}}_x\left({\widetilde{D}}_W+m_0\right)_{xy}\psi_y}+\\ &\frac{\rm i}{4u_0^4}\sum\limits_{x,\mu\nu}{\bar{\psi}\left(x\right)\sigma_{\mu\nu}\mathcal{F}_{\mu\nu}\left(x\right)\psi\left(x\right)}{\text{，}} \end{split}$$ (9) 这里的参数
$ u_0 $ 是tadpole改进参数，$ {\widetilde{D}}_W $ 与Wilson作用量中定义的一样，只不过这里的规范场$ U_\mu(x) $ 都换成了$ U_\mu(x)/u_0 $ 。$ u_0 $ 的取值为wilson plaquette的值开四次方，需要在产生组态的过程进行调节。需要指出的是，Wilson Clover作用在有限格距下破坏了手征对称性，这会导致一个相加性的夸克质量重整化。通常u/d夸克的裸质量
$ m_0 $ 是一个负值，从而Wilson Dirac算符$ D_W + m_0 $ 可能会出现接近于0的本征值，这种情况下在HMC算法中对Dirac算符求逆时会出现趋于无穷大的值，因此在HMC演化过程中可能会出现一些“exceptional”的组态。在这些“exceptional”组态上计算的物理量会出现明显的异常，下面我们计算$ {\rm{{\text{π}}}}$ 介子和K介子的两点关联函数时观察到了少数几个“exceptional”组态，在数据分析中去除了这些组态。 -
格点作用量中的参数包括裸耦合常数
$ \beta $ ，u/d夸克裸质量参数$ m_{\rm u} $ 和s夸克裸质量参数$ m_{\rm s} $ ，以及tadpole改进参数$ u_0 $ 。其中$ \beta $ 与格距有关，我们采用文献[10]中的办法，通过计算Wilson flow来确定格距。为了得到0.1 fm左右的格距，我们取$ \beta \!=\! 6.2 $ 。调节u/d夸克裸质量参数使$ {\rm{{\text{π}}}}$ 介子质量在250 MeV左右，s夸克质量通过K介子质量来调节，使K介子质量在物理的K质量附近。tadpole改进参数$ u_0 $ 的值取wilson plaquette的值开四次方，在组态产生过程中需要不断调整，使得输入的$ u_0 $ 值与用产生的组态计算的wilson plaquette开四次方的值吻合。我们从任意一个组态开始，用HMC算法对组态进行更新，得到一个新的组态，然后这样不断进行下去产生足够多的组态，这些组态最终会满足玻尔兹曼分布，每更新一次称作一个trajectory。由于我们是从任意一个组态开始的，首先需要预热一段时间使组态分布达到平衡，即达到我们想要的分布，在达到平衡之前的组态不能用来测量物理量。判断组态是否达到平衡通常可以通过观察某一个物理量随着蒙特卡罗时间的变化来判断，在达到平衡之后物理量会稳定在一个值。图3展示了wilson plaquette的值随着蒙特卡罗时间的演化，我们可以看到从一开始Wilson plaquette不断变化，最后在大约1000个trajectory之后稳定在一个值，图中在trajectory number为600的时候出现的变化是由于我们在这里调整了参数所致。由于新的组态都是由上一个组态得到，相邻的trajectory一定会存在一定的关联性，我们需要统计上相互独立的组态，因此，我们需要每隔一定数量的trajectory保存一个组态用来测量物理量。我们的做法是每隔5个trajectory保存一个组态，但这只是一个经验的选择，我们仍然需要通过测量物理量来判断这些组态之间是否存在自关联。下一节中我们计算
$ {\rm{{\text{π}}}}$ 介子两点关联函数，将通过关联函数的误差分析来判断自关联性。
图 3 (在线彩图)Wilson plaquette 随蒙特卡罗时间的演化
通过调节之后我们最终取的参数为：
$ m_{\rm u} \!=\! -0.279\,0 $ ,$ m_{\rm s} \!=\! -0.240\,0 $ ，$ u_0 $ =0.855 52。从产生的组态计算Wilson plaquette开四次方的值为0.855 55，只在第五位小数上与输入的$ u_0 $ 值有差异，吻合较好。计算Wilson flow得到的格距为0.105 fm。在硬件上，格点QCD产生组态的程序通常需要GPU加速卡的支持，计算规模和时间与格子体积大小、夸克质量、夸克作用量以及格距都有关系。本文的计算工作使用Nvidia V100 GPU卡，对于
$ 32^3 \times 64 $ 大小的格子，至少需要16块V100 GPU卡并行，产生一个组态需要70 min左右。 -
我们知道
$ {\text{π}}^{+} $ 的价夸克由一个u夸克和一个反d夸克构成，因此可以用算符$ \overline{d}_x\gamma_5u_x $ 来表示$ {\text{π}}^{+} $ 介子：$$ {\text{π}}^+(x) = \overline{d}_x\gamma_5u_x {\text{，}} $$ (10) 这里略写了Dirac指标和色指标。这个算符的厄密共轭为
$$ ({\text{π}}^+(x))^{\dagger} = (\overline{d}_x\gamma_5u_x)^{\dagger} = -\overline{u}_x\gamma_5d_x = -{\text{π}}^-_{x} {\text{。}} $$ (11) K介子与
$ {\rm{{\text{π}}}}$ 介子类似，将$ {\rm{{\text{π}}}}$ 介子中的反d夸克换成反s夸克即可得到K介子算符：$ K^{+}(x) = \overline{s}_x\gamma_5u_x $ 。下面我们将以$ {\rm{{\text{π}}}}$ 介子为例说明介子质量的计算过程。$ {\rm{{\text{π}}}}$ 介子的关联函数的定义为$$ \begin{split} C^{\text{π}}(t,t_0)\!=& \big\langle\widetilde{{\text{π}}}^+(t,{0}) (\widetilde{{\text{π}}}^+(t_0,{0}))^{\dagger}\big\rangle\\ =& -\frac{1}{V_3}\!\!\sum\limits_{{{x}},{{y}}}\!\big\langle \overline{d}(t,{{x}})\gamma_5u(t,{{x}})\overline{u}(t_0,{{y}})\gamma_5d(t_0,{{y}})\big\rangle _{U,F}\\ =&\frac{1}{V_3}\sum\limits_{{{x}},{{y}}}\big\langle {\rm Tr}[d_b(t_0,{{y}})\overline{d}_a(t,{{x}})\gamma_5u_a(t,{{x}})\times\\ & \overline{u}_b(t_0,{{y}})\gamma_5] \big\rangle _{U} {\text{，}}\\[-13pt] \end{split} $$ (12) 这里
$ a $ 和$ b $ 表示色指标。式(12)的第二行表示对理论的所有规范场$ U $ 和所有费米子场$ F $ 求期望值，第三行是利用了Wick定理将对费米子场的期望值写成了相应费米子场的缩并，Tr表示对Dirac旋量指标求迹。这里的费米子场缩并(又称夸克传播子)实际上就是理论中费米子场矩阵的逆。由于费米子矩阵是一个巨大的稀疏矩阵，因此要完全计算它的逆矩阵在数值上几乎是不可能的，但是费米子矩阵的逆矩阵的特定矩阵元可以通过数值方法计算。我们把夸克场的所有指标统一记为
$ A,B,C,\cdots $ , 这里的A包括了Dirac指标$ \alpha $ ，色指标$ a $ 和时空指标$ (t,{{x}}) $ ，我们要求的夸克传播子为$$ \psi_A\overline{\psi}_B = \big[\mathcal{M}[U_\mu]\big]^{-1}_{A;B} {\text{。}} $$ (13) 上述矩阵元可以通过设定点源，然后求解线性方程得到。定义除了在给定点B之外都等于零的源场
$ \phi^{(B)}_{A} $ :$$ \phi^{(B)}_{A} = \delta_{A;B} {\text{，}} $$ (14) 这里的
$ \delta_{A;B} $ 表示一系列Kronecker符号的乘积。如果以$ \phi^{(B)} $ 为源求解线性方程$$ \mathcal{M}[U_\mu]\cdot X^{(B)} = \phi^{(B)} {\text{，}} $$ (15) 这个线性方程的解
$ X^{(B)} $ 就是我们要求的传播子[11]$$ [X^{(B)}]_{A} = \big[\mathcal{M}[U_\mu]\big]^{-1}_{A;B} {\text{。}} $$ (16) 最为简单的零动量的
$ {\rm{{\text{π}}}}$ 介子关联函数经过设源求解线性方程之后可以用线性方程的解表达为式(17)的形式，因此它的数值可以由Monte Carlo计算获得。$$ C^{\text{π}}(t,t_0) = \frac{1}{V_3}\sum\limits_{\alpha,a,{x}}\sum\limits_{\beta,b}\Big\langle \left| X^{(\beta,b,t_0)}_{\alpha,a,t,{x}}\right| ^{2}\Big\rangle _{U} {\text{。}} $$ (17) 另一方面，假定
$ \tau = t-t_0>0 $ ，有$$ \begin{split}& C^{\text{π}}(\tau) = \frac{1}{Z}{\rm Tr}\big[{\rm e}^{-(T - \tau)H_{\widetilde{{\text{π}}}^+}}(0,{0}){\rm e}^{-\tau H}[\widetilde{{\text{π}}}^+(0,{0})]^{\dagger}\big]\\ =& \frac{1}{Z}\sum\limits_{n = 0}^{\infty}{\rm e}^{-(T - \tau)E_n}\Big\langle \varOmega_n\left| \widetilde{{\text{π}}}^+(0,{0}){\rm e}^{-\tau H}[\widetilde{{\text{π}}}^+(0,{0})]^{\dagger} \right| \varOmega_n\Big\rangle {\text{，}} \end{split}$$ (18) 这里的
$ \left| \varOmega_n\right\rangle,n = 0,1,\cdots $ 表示哈密顿量H的一组完备的能量本征态，其能量本征值记为$ E_n $ ，即$ H\left| \varOmega_n\right\rangle = E_n\left| \varOmega_n\right\rangle $ ，$ T $ 为格点上时间方向总长度。 我们约定$ \left| \varOmega_{n = 0}\right\rangle $ 。对应于QCD真空并且取$ E_0 = 0 $ ，同时假设其余的能量已经按照n的顺序排列好，即$ 0 = E_0<E_1<E_2\cdots $ 。对于
$ {\rm{{\text{π}}}}$ 介子的关联函数，在两个算符之间再插入一组完备的基，得到$$ C^{\text{π}}(\tau) = \frac{1}{Z}\!\!\sum\limits_{n,m = 0}^{\infty}\left| \langle \varOmega_n\left| \widetilde{{\text{π}}}^+(0,{0})\right| \varOmega_m\rangle \right| ^{2}{\rm e}^{-(T -\tau)E_{n}}{\rm e}^{-\tau E_m} {\text{。}} $$ (19) 因为所有其他
$ n\neq 0 $ 的态都具有比QCD真空更高的能量，上式说明所有这些态的贡献无论是沿着正的$ \tau $ 方向，还是沿着负的$ \tau $ 方向一定都是指数衰减的。一般而言，因为$ \left| \langle \varOmega_0\left| \widetilde{{\text{π}}}^{+}(0,{0})\right| \varOmega_0\rangle \right| = 0 $ ，所以$ n = m = 0 $ 的项是没有贡献的。所以衰减最慢的项是$ n = 0,\; m = 1 {\text{和}} $ $ n = 1,\; m = 0 $ 的项，即$$ \begin{split} C^{\text{π}}(\tau)\approx &\frac{1}{Z}\left| \langle \varOmega_0\left| \widetilde{{\text{π}}}^+(0,{0})\right| \varOmega_1\rangle \right|^{2}[{\rm e}^{-(T -\tau)E_1}+{\rm e}^{-\tau E_1}]\\ \approx &\frac{\left| \langle \varOmega_0\left| \widetilde{{\text{π}}}^+(0,{0})\right| \varOmega_1\rangle \right|^{2}}{\cosh(T E_1/2)}\cosh[(T/2-\tau)E_1] {\text{。}} \end{split}$$ (20) 对于像
$ {\rm{{\text{π}}}}$ 介子这样的强子，它的关联函数对于$ \tau $ 的依赖是一个双曲余弦的形式，对称点在时间方向的一半，即$ T/2 $ 的地方。我们可以对Monte Carlo数值模拟中计算得到的$ C^{\text{π}}(\tau) $ 进行拟合从而获得$ E_1 $ 的中心值及其误差。在格点QCD中，我们经常构造一个被称为有效质量的“物理量”，它的具体结构如下：
$$ R(\tau) = \frac{C^{\text{π}}(\tau+1)+C^{\text{π}}(\tau-1)}{{2C}^{\text{π}}(\tau)} {\text{。}} $$ (21) 可以证明，如果
$ C^{\text{π}}(\tau) $ 严格具有双曲余弦的形式，那么这样构造的$ R(\tau) $ 应当是一个不依赖于$ \tau $ 的常数$ \cosh(E_1) $ 。 因此我们可以构造有效质量$ m_{\rm eff}(\tau) $ 如下：$$ m_{\rm{eff}}(\tau) = \cosh^{-1}{\left[ \frac{C^{\text{π}}(\tau+1)+C^{\text{π}}(\tau-1)}{{2C}^{\text{π}}(\tau)}\right] } {\text{。}} $$ (22) 在上一节中我们提到组态之间可能存在自关联，可以通过对样本做bin之后的误差变化来观察是否存在自关联。假设原有
$ N $ 个样本，对这些样本每$ n $ 个取平均值，得到$ N/n $ 个新的样本，$ n $ 就是Bin size。图4中展示了$ {\rm{{\text{π}}}}$ 介子的两点关联函数的误差随Bin size的变化趋势，图中横轴是Bin size，纵轴是$ {\rm{{\text{π}}}}$ 介子两点关联函数在时间$ \tau \!=\! 10 $ 的误差和信号的比值，即$ \delta C^{\text{π}}(\tau \!=\! 10)/ C^{\text{π}}(\tau \!=\! 10) $ 。可以看到，随着Bin size的增大，误差一直在增大，直到Bin size等于10左右，趋于稳定。这说明我们产生的组态存在自关联，自关联长度约为10。因此我们取Bin size为10，对原始样本做bin之后来对$ {\rm{{\text{π}}}}$ 介子的关联函数进行分析。我们总共产生了548个组态，bin之后的样本数为54。
图 4 (在线彩图)
$\pi$ 介子两点关联函数在$\tau=10$ 的误差随着bin size的变化图5是对
$ {\rm{{\text{π}}}}$ 介子的关联函数$ C^{\text{π}}(\tau) $ 按照$ \cosh $ 函数形式，即式(20)，进行拟合得到的拟合曲线。拟合得到的能量$ aE \!=\! 0.117\pm0.002 $ ，根据上一节定的格距$ a \!=\! 0.105 $ fm，得到$ {\text{π}}$ 介子质量为$ (220\pm3) $ MeV。类似地，对K介子关联函数进行拟合，如图6所示，得到K介子质量为$ (490\pm 1) $ MeV，与物理的K介子质量接近。在图7和图8中分别给出了$ {\rm{{\text{π}}}}$ 介子和K介子的有效质量图。对有效质量做常数拟合得到的$ {\rm{{\text{π}}}}$ 介子和K介子与用$ \cosh $ 函数形式拟合关联函数得到的结果相同。
图 5 (在线彩图)
$\pi$ 介子两点关联函数$C^\pi(\tau)$ 随时间$\tau$ 的函数关系
图 6 (在线彩图)K介子两点关联函数
$C^{\rm K}(\tau)$ 随时间$\tau$ 的函数关系
图 7 (在线彩图)
$\pi$ 介子有效质量$m^K_{\rm eff}(\tau)$ 随时间$\tau$ 的函数关系
图 8 (在线彩图)K介子有效质量
$m^{\rm K}_{\rm eff}(\tau)$ 随时间$\tau$ 的函数关系 -
本文介绍了利用开源软件Chroma用HMC算法产生组态的过程。我们产生了一组具有格距0.105 fm，体积
$ 32^3 \times 64 $ ，$ {\rm{{\text{π}}}}$ 介子质量为220 MeV的组态，格点上夸克作用量采用Wilson clover作用量。格点计算中的系统误差主要来自有限格距、有限体积以及非物理的$ {\rm{{\text{π}}}}$ 介子质量。需要在不同的格距、体积和$ {\rm{{\text{π}}}}$ 介子质量下进行计算来系统地消除这些误差。我们将来计划产生多组具有不同格距、体积和$ {\rm{{\text{π}}}}$ 介子质量的组态，尤其是具有物理的$ {\rm{{\text{π}}}}$ 介子质量的组态。本文中的方法可直接用于将来的工作。致谢 感谢广州南方核科学计算中心提供的帮助，本论文使用的全部计算资源均由其提供。
-
摘要: 组态是格点QCD计算的基础，本文利用开源软件Chroma产生了一组格点QCD组态，格距为0.105 fm，体积为
$32^3\times 64$ ，π介子质量为220 MeV，格点上的夸克作用量采用Wilson clover作用量。这组组态可用于格点QCD中研究核子结构和强子谱等物理问题。Abstract: Configurations are the basis of lattice QCD calculations. In this paper, we generate a set of configurations with the Wilson clover quark action using the open source software Chroma. The parameters are tuned such that the lattice spacing is ~0.105 fm, pion mass is ~220 MeV. The volume is$32^3 \times 64$ . These configurations will be used to study nucleon structure and exotic hadronic states in lattice QCD.-
Key words:
- lattice QCD /
- strong interactions /
- chroma /
- lattice QCD configurations /
- pion
-
图 1 (在线彩图)格点QCD基本方法示意
格子具有有限大小$ L $和有限格距$ a $，将QCD拉氏量离散化之后得到格点QCD拉氏量，夸克场$ \varPsi $定义在格点上，胶子场$ U $在连接相邻格点的连接线上。
-
[1] ACCARDI A, ALBACETEET J L, ANSELMINO M,et al. Eur Phys J A, 2016, 52(9): 268. doi: 10.1140/epja/i2016-16268-9 [2] ABDUL K R, ACCARDI A, ADAM J, et al. arXiv: 2103.05419. [3] 曹须, 常雷,畅宁波, 等. 核技术, 2020, 43(02): 3. doi: 10.11889/j.0253-3219.2020.hjs.43.020001 CAO Xu, CHANG Lei, CHANG Ningbo, et al. Nuclear Technology, 2020, 43(02): 3. (in Chinese) doi: 10.11889/j.0253-3219.2020.hjs.43.020001 [4] ANDERLE D P, BERTONE V, CAO X, et al. arXiv: 2102.09222. [5] WILSON K G. Phys Rev D, 1974, 10: 2445. doi: 10.1103/PhysRevD.10.2445 [6] EDWARDS R G, JOO B. Nucl Phys B Proc Suppl, 2005, 140: 832. doi: 10.1016/j.nuclphysbps.2004.11.254 [7] 田英齐, 毕玉江, 贺雨晴, 等. 计算机系统应用, 2019, 028(009): 25. TIAN Yingqi, BI Yujiang, HE Yuqing, et al. Computer Systems and Applications, 2019, 028(009): 25. (in Chinese) [8] CLARK M A, BABICH R, BARROS K, et al. Comput Phys Commun, 2010, 181: 1517. doi: 10.1016/j.cpc.2010.05.002 [9] SHEIKHOLESLAMI B, WOHLERT R. Nucl Phys B, 1985, 259: 572. doi: 10.1016/0550-3213(85)90002-1 [10] BORSANYI S, DÜRR S, FODOR Z, et al. JHEP, 2012, 09: 010. doi: 10.1007/JHEP09(2012)010 [11] 刘川. 格点量子色动力学导论[M]. 北京: 北京大学出版社, 2017. LIU Chuan. An Introduction to Lattice Chromodynamics[M]. Beijing: Peking University Press, 2017 (in Chinese) -
点击查看大图
图(8)
计量
- 文章访问数: 1233
- HTML全文浏览量: 386
- PDF下载量: 87
- 被引次数: 0
甘公网安备 62010202000723号