Perturbation Calculation of Charge Radii in Nuclear Lattice Effective Field Theory

从现实核力出发的第一性原理计算(ab initio calculations)已成为当今核理论领域的重要前沿研究方向。它不仅能够提供高精度的计算结果与可靠的理论预言，同时由于没有额外的参数依赖，能够在有效场论的框架下建立起与量子色动力学(QCD)联系的桥梁。因此第一性原理计算对核理论的发展、核实验的解释及预言、核天体环境下元素合成过程的理解以及中子星结构的研究等都非常重要。近些年来随着手征有效核力的逐步完善^[1−6]与量子多体方法的不断发展，第一性原理计算取得了许多令人瞩目的成就。较早出现的无芯壳模型(NCSM)^[7]与量子蒙特卡罗(QMC)方法^[8]能够对较轻的原子核进行严格计算，但由于其计算量随核子数呈指数增加，难以触及更重的原子核体系。现实核力(Realistic Nuclear Forces)壳模型(SM)^[9−10]通过多体微扰的方法构建出模型空间哈密顿量, 在保证精度的同时有效地降低了计算量, 拓展了其理论描述的能力。非微扰的Coupled Cluster^[11]和介质中相似重整化群IM-SRG^[12−13]将满壳或半满壳核作为参考态, 以此能够对中重核进行理论计算与预言。而价空间的IMSRG^[14]，则能够对开壳原子核进行可靠的描述。通过采用Gamow-Berggren基，现实核力Gamow壳模型^[15−25]与Gamow IM-SRG^[26]能够将连续谱效应考虑进来对弱束缚区原子核进行有效的理论计算。

原子核格点有效场论(Nuclear Lattice Effective Field Theory, NLEFT)是近些年发展起来的一种第一性原理计算方法^[27−37]。它将有效场论与格点蒙特卡罗方法有机地结合在一起。相较于其他第一性原理方法，NLEFT最大的优势是在格点的框架下能够同时构建起核子间相互作用与求解量子多体问题。得益于蒙特卡罗方法的特性，它的计算复杂度随核子数成多项式增长，同时能够提取诸如核集团结构^{[34, 37]}、α-α散射^[32]等其它第一性原理方法难以计算的强关联信息。

与其它量子蒙特卡罗方法一样，NLEFT也面临着“符号问题”的困难。在NLEFT中，我们使用辅助场(Auxiliary field)方法将核子间相互作用转变为核子与“背景场”间的相互作用。但核力的排斥部分会使得核子与“背景场”间的相互作用出现振荡的虚部相位，从而降低蒙卡计算中的统计信号，导致计算结果出现很大的统计误差。为此，原子核格点有效场论合作组发展出波函数匹配(Wave Function Matching, WFM)方法^[38]与微扰量子蒙特卡罗方法(perturbative Quantum Monte Carlo method, pQMC)^[36]，其核心是构建出一个幺正变换，它能将高精度的手征哈密顿量H变换到$ H' $，其中$ H' $靠近一个没有符号问题的、容易处理的哈密顿量$H^0$。因此，较小的$ H' - H^0 $部分则可以通过微扰论的方法进行计算。基于以上方法，我们与合作者发展了Rank-one Operator (RO)方法^[39]对算符进行微扰计算。

在本文中，我们进一步发展格点上的微扰计算方法，通过引入两点关联函数进行原子核电荷半径的计算。在第1节，我们将对格点有效场论中的微扰计算进行简要介绍；在第2节，给出在格点上、半径算符的计算方法；在第3节，将讨论微扰计算的收敛性与有限体积效应，并给出使用完整N³LO核力的计算结果；最后作出总结与展望。

格点有效场论是一种基于时空格点离散化的量子蒙特卡罗方法^{[27, 31]}。它使用虚时投影得到体系的基态波函数$ |\varPsi\rangle=\lim _{L_t \rightarrow \infty} M^{L_t / 2}|\varPsi_T\rangle $，其中$ |\varPsi_T\rangle $为试探波函数，而$ M $是每一个虚时时间隔$ a_t $所对应的转移矩阵$ M=: e^{-a_tH}: $。冒号表示算符正规乘积。在得到基态波函数后，算符的期望通过振幅$ {\cal{M}} $的比值求得，$ \langle {\cal{O}} \rangle = {{\cal{M}}({\cal{O}})}/{{\cal{M}}} $，其中振幅$ {{\cal{M}}({\cal{O}})}=\langle \varPsi \vert {{\cal{O}}} \vert \varPsi \rangle $，代表了包含有算符的关联矩阵行列式(determinant of correlation matrix)。通过辅助场(Auxiliary Field)的引入，我们可以严格地将哈密顿量H中两体和三体相互作用转化为单个粒子与辅助场间的相互作用。但经过辅助场变换后，核力中的排斥部分会造成计算中出现复的相因子。这个相因子会导致在蒙卡统计中出现强烈的正负相消，也就是著名的"符号问题"。迄今为止还没有解决"符号问题"的一般方案。所有解决方法都仅适用于特定系统。

为了克服“符号问题”所带来的计算困难，我们发展了WMF^[38]方法和pQMC^[36]方法。从高精度的手征哈密顿量$ H_{{\rm{N3LO}}} $出发，WMF方法构建出一组幺正变换将$ H_{{\rm{N3LO}}} $变为$ H'_{{\rm{N3LO}}} $。这组幺正变换使得在两体系统下，$ H'_{{\rm{N3LO}}} $与一个易处理的$ H^S $的波函数在短程部分相匹配。此幺正变换确保了高精度哈密顿量与易处理的哈密顿量之间的差异$ H'_{{\rm{N3LO}}}-H^S $较小，并且能够通过微扰的方法进行计算。通常我们选取近似SU(4)对称性(Leading-order, LO)手征哈密顿量作为$ H^S $。更多相关细节，请查看文献[36, 38−39]。

为了计算格点上的微扰波函数与算符，我们将哈密顿量$ H $拆分成非微扰与微扰两个部分$ H=H_0+H_1 $，与之相对应的转移矩阵为$ M=M_0+M_1 $。对于基态波函数$ |\varPsi\rangle=\lim _{L_t \rightarrow \infty} (M_0+M_1)^{L_t / 2}|\varPsi_T\rangle $，我们有0阶非微扰部分：

与1阶微扰修正：

对于能量期望我们有$ E = E^{(0)}+ E^{(1)}+\cdots $，其中

以及

而算符期望则相对复杂，包含0阶非微扰部分：

与1阶修正：

其中共有四种不同类型的振幅：$ {\cal{M}}^{(0)} $、$ {\cal{M}}^{(1)} $、$ {\cal{M}}_o^{(0)} $与$ {\cal{M}}_o^{(1)} $。前三种振幅的计算相对容易处理，这里$ {\cal{M}}^{(0)} $为非微扰波函数的关联矩阵行列式$ {\mathcal{M}}^{(0)}=\det {\mathcal{A}}^{(0)} $, 关联矩阵$ {\cal{A}}^{(0)} $的矩阵元由单粒子波函数积分求得$ {\cal{A}}^{(0)}_{ij} = \langle \phi_i \vert \phi_j \rangle $。类似地$ {\cal{M}}^{(1)} $为微扰波函数的关联矩阵行列式。对于包含有算符$ {\cal{O}} $的振幅$ {\cal{M}}_o^{(0)} $，我们可以通过Jacobi方法求得：

更多细节请查看文献[27, 31]。真正的困难出现在微扰的转移矩阵$ M_1 $与算符$ {\cal{O}} $同时存在的$ {\cal{M}}_o^{(1)} $。最直接的原因来自于多于三体的算符插入将会使得Jacobi公式变得非常复杂，而微扰转移矩阵$ M_1 $中已经包含了手征两体与三体算符，与此同时我们还需要考虑可观测量的算符插入(例如两点关联函数算符)，这将使计算变得异常困难。

为此我们与合作者发展出RO (Rank-one Operator)方法^[39]。在这个方法中，单体算符可以被拆分成秩为1的组分，$ {\cal{O}} = \sum_{\alpha\alpha',\, is} F_{\alpha',\, is}^{\dagger}F_{\alpha,\, is} $，其中$ \alpha $为轨道标记，$ F_\alpha= \sum_{{\boldsymbol{n}}} a_{is}({\boldsymbol{n}}) f_{\alpha,\, is}({\boldsymbol{n}}) $，$ F_\alpha^{\dagger}=\sum_{{\boldsymbol{n}}} a_{\alpha, is}^{\dagger}({\boldsymbol{n}}) f_{\alpha, is}^*({\boldsymbol{n}}) $以及代表同位旋与自旋的$i,\; s$。可以看出$ F_{\alpha',\, is}^{\dagger}F_{\alpha,\, is} $的秩为1，任何指数高于1的$ F_{\alpha',\, is}^{\dagger}F_{\alpha,\, is} $组合在normal ordering后作用在单体波函数上将为0。因此我们有关$ F_{\alpha',\, is}^{\dagger}F_{\alpha,\, is} $线性关系式：

并可以应用到求解$ {\cal{M}}_o^{(0)} $中：

对于求解$ {\cal{M}}_o^{(1)} $，我们能够通过RO方法与Jacobi方法分别处理算符$ {\cal{O}} $与微扰转移矩阵$ M^{(1)} $的插入，最终对其进行计算。值得一提的是，Rank-one Operator方法拓展到两体或多体算符非常直接，并且原则上能够应用到任何依赖于使用行列式做反对称化的量子多体方法中。更多有关细节请参考文献[39]。

原子核电荷半径的理论计算通常包含以下几项的贡献^[40]：

质子相对于原子核质心的半径、质子半径、中子半径修正以及相对论修正项。本文中质子半径$ R_{\mathrm{p}} = 0.840^{+0.003+0.002}_{-0.002-0.002}\; \mathrm{fm} $^[41]，中子半径$ R_{\mathrm{n}}^2= -0.105^{+0.005}_{-0.006}\; \mathrm{fm}^2 $^[42], 以及$\frac{3 \hbar^2}{4 m_{\mathrm{p}}^2 c^2} \approx 0.033~172 \;\mathrm{fm}^2$。除此之外，自旋轨道耦合也会带来很小的相对论修正，这部分的贡献将在未来的工作中加入。电荷半径的主要贡献来自于质子相对于原子核质心的半径$ r_{\rm pp} $，这也是计算中的难点。

在NLEFT中目前有两种方案来计算$ r_{\mathrm{pp}} $：Pinhole方法与两体关联函数。Pinhole方法中^[34]，我们通过引入Pinhole来确定每个核子的位置坐标，进而计算质心位置与半径。但是由于Pinhole方法通常会引入更大的“符号问题”，在本文中我们选择使用两体关联函数来计算质子相对于质心的半径：

其中：$ \rho(\boldsymbol{r}) $为单体密度算符；$ |\boldsymbol{r}-\boldsymbol{r}^{\prime}| $则代表了两个密度测量点间的距离。

在具有$ L^3 $个格点的立方盒子中，我们采取周期边界条件，并有$ -{L}/{2}< r_x, r_y, r_z \leqslant{L}/{2} $。在格点上计算两点距离时我们也通常要求$ -{L}/{2}< r_x - r'_x \leqslant{L}/{2} $。可以看到，当格点$ L $不足够大时由于$ | r_x - r'_x |\leqslant{L}/{2} $会使测量半径偏小，引入额外的有限体积效应(Finite Volume Effects)。为此，我们发展出一个新的方法来测量两点间距离，包含有两个步骤：1) 通过计算在$ x,\, y,\, z $方向上的单体密度并找到大致的质心位置$ (x_{cm},\, y_{cm},\, z_{cm}) $；2) 两点间的距离由$ | {\rm{Lmod}}(r_x - x_{cm}) - {\rm{Lmod}}(x_{cm} - r'_x) | $给出，其中$ {\rm{Lmod}} $会将数值约化到$ (-{L}/{2},\, {L}/{2}] $区间。在下一个小节中我们将看到，使用新定义的格点上的两点距离将有效地减小有限体积效应。

为了在格点上进行半径计算，我们首先定义点状粒子半径算符：

其中：$ {{\boldsymbol{x}}} $和$ {{\boldsymbol{y}}} $代表三维格点坐标；指标i代表同位旋：质子、中子或不加指标代表质子加中子。因此我们可以写出：

与

我们结合RO方法与式(7)第一行中的求导方法，可以得到关于算符$ \hat{{\cal{R}}}^2_i $的振幅计算公式：

其中$ \hat{\rho}_{is} = a^{\dagger}_{is}a_{is} $是秩为1的算符，为了公式简洁，我们新定义了算符$ \hat{{\cal{P}}} = \sum_{{{\boldsymbol{x}}}} \hat{\rho}[{{\boldsymbol{x}}}+{{\boldsymbol{y}}}] {{\boldsymbol{r}}}^2[{{\boldsymbol{x}}}] $。

本文中我们将首次应用上述方案计算原子核半径，同时结合第1节中所概述的微扰蒙特卡罗方法，使用高精度的手征N³LO核力对原子核半径进行微扰计算。在以下计算中，我们设置格点的间距$ a=1/150 \; {\rm{MeV}}^{-1}\approx 1.32 \; {\mathrm{fm}}$以及虚时间隔$ a_t=1/1000 \;{\rm{MeV}}^{-1} $。

为了验证半径微扰计算的可靠性，我们首先使用一个简单模型进行验证。在这个模型中，哈密顿量$ H $包含两部分$ H=T+V $，其中$ T $是动能项，$ V $是两体接触项，耦合强度为$ -4.2\times 10^{-7} $MeV$ ^{-2} $。我们将$ H $拆分成非微扰的$ H_0=T + (1-x)V $与微扰的$ H_1=xV $，同时引入了参数$ x $控制微扰部分的大小。

这里我们选取$ L^3=6^3 $的格点以及$ L_t = 100 $对¹⁶O进行微扰验证计算。当$ x=0 $时，非微扰计算将给出完整哈密顿量下的$ r_{\rm pp} $，而微扰计算中的微扰贡献为零。当$ x=1 $时，非微扰计算中将只有动能项，相互作用项将全部由微扰修正给出。在表1中，我们列出$ x=0.0,\, 0.2,\, 0.4 $与1.0的能量1阶微扰计算。可以看出当$ x=0.2 $时，一阶微扰已经能给出约40%的能量修正。即使当$ x=0.4 $，微扰修正后的能量与准确值间的差别也仅为4%。

在图1中，我们给出与之相对应的半径微扰计算结果。其中当$x=0$时，计算得到的半径为使用完整哈密顿量的"准确值"，我们将其作为参考并在图中以灰色虚线表示。可以看出，当$ x=0.2 $时，非微扰计算给出的$ r_{\rm pp} $显著偏离正确值，而微扰计算则能将其很好的修正。随着$ x $增加，非微扰的$ r_{\rm pp} $呈线性增长并逐渐偏离准确值。通过算符的微扰计算，$ r_{\rm pp} $则会被极大地修正到准确值附近。即便当体系中的相互作用项全部为$ H_1 $时，$ E_1 $超过两倍的体系总能量，一阶微扰修正后的$ r_{\rm pp} $与准确值间的差别也不超过2%。可以看出，此方法中，微扰对半径的修正非常准确。

在上一小节中，为了消除有限体积效应，我们引入一个新的计算格点上距离的方法。这里我们采用了领头阶(Leading Order)的手征核力，并设定$ L_t=100 $，对¹⁶O进行非微扰的计算验证。在图2中可以看到，使用常规方法计算的质子相对于质心的半径$ r_{\rm pp} $在格点$ L $较小时，明显偏小，直到$ L=10 $才足够接近无穷大空间中的极限值。而使用新方法计算的结果则明显地减少了有限体积效应。

最后，我们使用两体N³LO核力与三体N²LO核力，对⁴He、¹²C与¹⁶O进行了完整的微扰计算，并将结果列在表2中。这里我们采用$ L=10 $的格点，以及$ L_t=200 $作为参数进行计算。值得一提的是表2中的误差仅为蒙特卡罗采样的统计不确定度，更多关于手征核力的截断误差、$ L_t $外推以及有限体积效应等系统误差将在未来的工作中讨论。

从表2中我们可以看出，我们的NLEFT计算能够给出与实验符合非常好的原子核结合能。其中最大的差别出现在¹²C，我们计算的结合能比实验大3.5%。同时，通过新发展的电荷半径计算方法，我们能够非常好地复现⁴He与¹²C的实验测量结果，给出的¹⁶O的计算结果仅稍大于实验值约5%。以上结果均验证了本文中新发展的格点上微扰计算方法的可靠性。

本文中，我们发展了在格点上微扰计算原子核半径的新方法。通过使用WMF和pQMC方法，避免了严重的“符号问题”，并结合两点关联函数，能够有效地测量半径算符。为了减小格点的有限体积效应，我们用了新的方法计算格点上的两点距离。验证计算表明了格点上半径微扰计算的可靠性。接着，我们使用完整的N³LO手征核力对⁴He、¹²C与¹⁶O的结合能与半径进行了计算，结果很好地复现了实验数据。该方法将成为新的第一性原理计算原子核半径的有效工具。未来，我们将使用该方法对更多的原子核体系进行高精度的格点计算。