Improving the Kinetic Couplings in Lattice Non-Relativistic QCD

TL;DR

This paper enhances lattice NRQCD by calculating higher-order kinetic corrections, determining one-loop matching coefficients, and testing the impact on bottomonium mass calculations using multiple lattice spacings.

Contribution

The work introduces improved kinetic operators in lattice NRQCD with one-loop corrections and assesses their effects on bottomonium mass determinations.

Findings

Improved dispersion relation accuracy up to (p^6)

One-loop matching coefficients for kinetic operators determined

Enhanced bottomonium mass precision with new corrections

Abstract

We improve the non-relativistic QCD (NRQCD) action by comparing the dispersion relation to that of the continuum through $\mathcal{O}(p^6)$ in perturbation theory. The one-loop matching coefficients of the $\mathcal{O}(p^4)$ kinetic operators are determined, as well as the scale at which to evaluate $\alpha_s$ in the $V$-scheme for each quantity. We utilise automated lattice perturbation theory using twisted boundary conditions as an infrared regulator. The one-loop radiative corrections to the mass renormalisation, zero-point energy and overall energy-shift of an NRQCD $b$-quark are also found. We also explore how a Fat$3$-smeared NRQCD action and changes of the stability parameter $n$ affect the coefficients. Finally, we use gluon field ensembles at multiple lattice spacing values, all of which include $u$, $d$, $s$ and $c$ quark vacuum polarisation, to test how the improvements affect the non-perturbatively determined $\Upsilon(1S)$ and $\eta_b(1S)$ kinetic masses, and the tuning of the $b$ quark mass.

Tables

Table 1. Table 1: The one-loop mean Landau gauge link Horgan:Tads and its log moments for either an unsmeared or smeared gauge-link definition.

Gauge-Link	$u_0^{\left(2\right)}=⟨f^{tad}⟩$	$⟨f^{tad}\log \left(q^2\right)⟩$	$⟨f^{tad}{\log }^2\left(q^2\right)⟩$
Unsmeared	$0.750275\left(5\right)$	$1.45755\left(2\right)$	$3.6022\left(1\right)$
Fat$3$	$0.231784\left(5\right)$	$0.26101\left(7\right)$	$0.6429\left(2\right)$
Fat$7$	$0.108244\left(5\right)$	$0.10271\left(4\right)$	$0.4117\left(2\right)$

Table 2. Table 2: The different little groups relevant for each momentum type in a finite-volume cubic lattice with PBC. The momenta are in units of 2 ​ π / L 2 𝜋 𝐿 2\pi/L and n , m , p 𝑛 𝑚 𝑝 n,m,p are non-zero integers with n ≠ m ≠ p 𝑛 𝑚 𝑝 n\neq m\neq p Moore:2005 ; Thomas:Helicity . The single cover irreps describe integer spin states.

Momentum Type	Little Group	Irreps
	(Double Cover)	(Single Cover)
$(0,0,0)$	$O_h^D$	$A_1^{\pm },A_2^{\pm },E^{\pm },T_1^{\pm },T_2^{\pm }$
$(n,0,0)$	$Di{c}_4$	$A_1,A_2,E_2,B_1,B_2$
$(n,n,0)$	$Di{c}_2$	$A_1,A_2,B_1,B_2$
$(n,n,n)$	$Di{c}_3$	$A_1,A_2,E_2$
$(n,m,0)$	$C_4$	$A_1,A_2$
$(n,n,m)$	$C_4$	$A_1,A_2$
$(n,m,p)$	$C_2$	$A$

Table 3. Table 3: Details of the gauge ensembles used in this study. β 𝛽 \beta is the gauge coupling. a Υ subscript 𝑎 Υ a_{\Upsilon} is the lattice spacing determined from the Υ ​ ( 2 ​ S − 1 ​ S ) Υ 2 𝑆 1 𝑆 \Upsilon(2S-1S) splitting Dowdall:Upsilon , where the error combines statistics, experiment and the dominant NRQCD systematic error. a ​ m q 𝑎 subscript 𝑚 𝑞 am_{q} are the sea quark masses, N s × N T subscript 𝑁 𝑠 subscript 𝑁 𝑇 N_{s}\times N_{T} gives the spatial and temporal extent of the lattices in lattice units and n cfg subscript 𝑛 cfg n_{\text{cfg}} is the number of configurations in each ensemble. In column 1 we use the numbering convention for the ensembles from Dowdall:Upsilon . Ensemble 1 1 1 is referred to as “very coarse”, 3 as “coarse,” and 5 as “fine”.

Set	$\beta$	$a_{{\rm Y}}$(fm)	$a{m}_l$	$a{m}_s$	$a{m}_c$	$N_s\times N_T$	$n_{\text{cfg}}$
$1$	$5.8$	$0.1474\left(15\right)$	$0.013$	$0.065$	$0.838$	$16\times 48$	$1020$
$3$	$6.0$	$0.1219\left(9\right)$	$0.0102$	$0.0509$	$0.635$	$24\times 64$	$1052$
$5$	$6.3$	$0.0884\left(6\right)$	$0.0074$	$0.037$	$0.440$	$32\times 96$	$1008$

Table 4. Table 4: Parameters used for the valence quarks. a ​ m b 𝑎 subscript 𝑚 𝑏 am_{b} is the bare b 𝑏 b -quark mass in lattice units, u 0 ​ L subscript 𝑢 0 𝐿 u_{0L} is the tadpole parameter Hughes:Hindered . T p subscript 𝑇 𝑝 T_{p} is the total propagation time for the b 𝑏 b -quark propagator and n t subscript 𝑛 𝑡 n_{t} is the number of time sources used per configuration. The 𝒪 ​ ( α s ) 𝒪 subscript 𝛼 𝑠 {\mathcal{O}(\alpha_{s})} matching coefficients for c 1 , c 6 subscript 𝑐 1 subscript 𝑐 6 c_{1},c_{6} and c 5 subscript 𝑐 5 c_{5} are taken from Tables 8 and 9 . As explained in Sec. III the 𝒪 ​ ( α s ) 𝒪 subscript 𝛼 𝑠 \mathcal{O}(\alpha_{s}) coefficients are functions of a ​ m b 𝑎 subscript 𝑚 𝑏 am_{b} ; the α s subscript 𝛼 𝑠 \alpha_{s} value they are combined with to give c 1 subscript 𝑐 1 c_{1} , c 5 subscript 𝑐 5 c_{5} and c 6 subscript 𝑐 6 c_{6} depends on the lattice spacing. The values are different for each version of the NRQCD action tested. As we focus on the improvements made in this study, c 2 , c 3 subscript 𝑐 2 subscript 𝑐 3 c_{2},c_{3} and c 4 subscript 𝑐 4 c_{4} are taken to be their tree-level values of 1.0 1.0 1.0 .

Set	$a{m}_b$	$u_{0L}$	$T_p$	$n_t$
$1$	$3.40$	$0.8195$	$40$	$16$
$3$	$2.80$	$0.8349$	$40$	$16$
$5$	$1.90$	$0.8525$	$48$	$16$

Table 5. Table 5: The local bilinear operators used in this study. Note the i ​ γ 5 𝑖 superscript 𝛾 5 i\gamma^{5} is needed to make the overlaps real Dudek:CharmRad . The second column gives the J P ​ C superscript 𝐽 𝑃 𝐶 J^{PC} states that these operators create at rest in an infinite volume continuum. The third column gives the helicity eigenvalues λ 𝜆 \lambda that these operators create at nonzero momentum in an infinite volume continuum which is only rotationally invariant, while the J 𝐽 J in brackets are the states which contribute to that helicity (c.f. Hughes:Hindered ; Thomas:Helicity .)

${𝒪}^{\Gamma }\left(x\right)$	$J^{PC}$	$\lambda \left(\leftarrow J^P\right)$
$\overline{\psi }i{\gamma }^5\psi$	$0^{-+}$	$0^{-}(\leftarrow J^P=0^{-},1^{+},2^{-},\dots )$
$\overline{\psi }{\gamma }^i\psi$	$1^{--}$	$0^{+}(\leftarrow J^P=0^{+},1^{-},2^{+},\dots )$
		$\left|1\right|(\leftarrow J=1,2,3,\dots )$

Table 6. Table 6: Numerical values of the one-loop coefficients. Note that the unsmeared results ( 𝒪 ​ ( p 4 ) 𝒪 superscript 𝑝 4 \mathcal{O}(p^{4}) and 𝒪 ​ ( p 6 ) 𝒪 superscript 𝑝 6 \mathcal{O}(p^{6}) ) are mean-field improved given the formulae in Appendix B and the mean field parameter in Table 1 . The smeared results (Fat3) are not mean-field improved.

$a{m}_b$	$7.0$	$4.0$	$3.4$	$2.8$	$1.9$	$1.1$
${\tilde{c}}_1^{\left(1\right)}$
$𝒪\left(p^4\right)$	$0.93051\left(24\right)$	$0.68266\left(17\right)$	$0.61998\left(16\right)$	$0.55180\left(13\right)$	$0.43443\left(15\right)$	$0.261360\left(99\right)$
$𝒪\left(p^6\right)$	$0.69607\left(24\right)$	$0.52875\left(17\right)$	$0.50348\left(16\right)$	$0.48171\left(13\right)$	$0.37427\left(10\right)$	$-1.5525\left(50\right)$
Fat$3$	$0.3991\left(37\right)$	$0.2504\left(27\right)$	$0.1914\left(24\right)$	$0.1012\left(21\right)$	$-0.1970\left(15\right)$	$-1.25492\left(97\right)$
$c_5^{\left(1\right)}$
$𝒪\left(p^4\right)$	$0.568\left(11\right)$	$0.5341\left(62\right)$	$0.5220\left(53\right)$	$0.5056\left(44\right)$	$0.4628\left(30\right)$	$0.2946\left(17\right)$
$𝒪\left(p^6\right)$	$0.07863\left(55\right)$	$0.11875\left(32\right)$	$0.13945\left(27\right)$	$0.17215\left(22\right)$	$0.26557\left(15\right)$	$0.2212\left(58\right)$
Fat$3$	$0.4813\left(87\right)$	$0.4338\left(50\right)$	$0.4100\left(42\right)$	$0.3727\left(35\right)$	$0.2569\left(24\right)$	$-0.1199\left(69\right)$
$Z_m^{\left(1\right)}$
$𝒪\left(p^4\right)$	$-0.08945\left(52\right)$	$0.08892\left(52\right)$	$0.15488\left(52\right)$	$0.24259\left(52\right)$	$0.44572\left(52\right)$	$0.77983\left(52\right)$
$𝒪\left(p^6\right)$	$-0.08988\left(52\right)$	$0.08413\left(52\right)$	$0.14947\left(52\right)$	$0.24053\left(52\right)$	$0.50827\left(52\right)$	$1.85962\left(52\right)$
Fat$3$	$0.08866\left(52\right)$	$0.22842\left(52\right)$	$0.28608\left(52\right)$	$0.36957\left(52\right)$	$0.60455\left(52\right)$	$1.23018\left(52\right)$
$W_0$
$𝒪\left(p^4\right)$	$-0.94595\left(52\right)$	$-0.84252\left(52\right)$	$-0.80232\left(52\right)$	$-0.74890\left(52\right)$	$-0.64040\left(52\right)$	$-0.72316\left(52\right)$
$𝒪\left(p^6\right)$	$-0.93235\left(52\right)$	$-0.82391\left(52\right)$	$-0.78118\left(52\right)$	$-0.72050\left(52\right)$	$-0.52587\left(52\right)$	$1.16267\left(52\right)$
Fat$3$	$-0.97946\left(52\right)$	$-0.98997\left(52\right)$	$-0.99508\left(52\right)$	$-1.00270\left(52\right)$	$-1.02215\left(52\right)$	$-1.03670\left(52\right)$
$\delta C$
$𝒪\left(p^4\right)$	$-0.22458\left(53\right)$	$-0.12171\left(54\right)$	$-0.08110\left(54\right)$	$-0.02487\left(55\right)$	$0.10867\left(59\right)$	$0.12242\left(70\right)$
$𝒪\left(p^6\right)$	$-0.22307\left(53\right)$	$-0.12185\left(54\right)$	$-0.08029\left(54\right)$	$-0.01679\left(55\right)$	$0.23149\left(59\right)$	$2.91658\left(70\right)$
Fat$3$	$-0.05126\left(53\right)$	$-0.01907\left(54\right)$	$-0.00659\left(54\right)$	$0.01146\left(55\right)$	$0.06658\left(59\right)$	$0.30555\left(70\right)$

Table 7. Table 7: Numerical values for the scale a ​ q ∗ 𝑎 superscript 𝑞 aq^{*} used to evaluate α s subscript 𝛼 𝑠 \alpha_{s} in the V-scheme when computing the one-loop contributions. Note that the unsmeared results ( 𝒪 ​ ( p 4 ) 𝒪 superscript 𝑝 4 \mathcal{O}(p^{4}) and 𝒪 ​ ( p 6 ) 𝒪 superscript 𝑝 6 \mathcal{O}(p^{6}) ) are mean-field improved given the formulae in Appendix B and the mean field parameter in Table 1 . The smeared results (Fat3) are not mean-field improved. The a ​ q ∗ 𝑎 superscript 𝑞 aq^{*} are determined as described in Sec. II.3 .

$a{m}_b$	$7.0$	$4.0$	$3.4$	$2.8$	$1.9$	$1.1$
$a{q}^{\ast }\left({\tilde{c}}_1^{\left(1\right)}\right)$
$𝒪\left(p^4\right)$	$2.4882\left(13\right)$	$2.17692\left(51\right)$	$2.06963\left(46\right)$	$1.93677\left(66\right)$	$1.6893\left(15\right)$	$1.56045\left(59\right)$
$𝒪\left(p^6\right)$	$2.6799\left(12\right)$	$2.22771\left(80\right)$	$2.10459\left(93\right)$	$1.9845\left(15\right)$	$1.7394\left(23\right)$	$2.978\left(11\right)$
Fat$3$	$2.459\left(22\right)$	$2.083\left(20\right)$	$1.961\left(21\right)$	$1.737\left(27\right)$	$2.056\left(14\right)$	$1.7764\left(19\right)$
$a{q}^{\ast }\left(c_5^{\left(1\right)}\right)$
$𝒪\left(p^4\right)$	$2.651\left(56\right)$	$2.585\left(30\right)$	$2.571\left(26\right)$	$2.560\left(22\right)$	$2.570\left(16\right)$	$2.510\left(15\right)$
$𝒪\left(p^6\right)$	$1.922\left(21\right)$	$2.035\left(30\right)$	$2.310\left(22\right)$	$2.651\left(34\right)$	$2.9942\left(87\right)$	$0.625\left(19\right)$
Fat$3$	$2.034\left(28\right)$	$1.937\left(16\right)$	$1.902\left(14\right)$	$1.854\left(12\right)$	$1.744\left(12\right)$	$1.92\left(21\right)$
$a{q}^{\ast }\left(Z_m^{\left(1\right)}\right)$
$𝒪\left(p^4\right)$	$1.375\left(35\right)$	$1.169\left(11\right)$	$1.0532\left(88\right)$	$1.1508\left(62\right)$	$1.1920\left(35\right)$	$1.0931\left(36\right)$
$𝒪\left(p^6\right)$	$1.370\left(21\right)$	$1.170\left(12\right)$	$0.9945\left(86\right)$	$1.1342\left(61\right)$	$1.3047\left(34\right)$	$1.7829\left(25\right)$
Fat$3$	$0.9598\left(80\right)$	$1.0007\left(35\right)$	$1.0416\left(29\right)$	$1.0784\left(23\right)$	$1.1316\left(24\right)$	$1.2024\left(25\right)$
$a{q}^{\ast }\left(W_0\right)$
$𝒪\left(p^4\right)$	$1.0833\left(30\right)$	$1.0370\left(32\right)$	$1.0163\left(33\right)$	$0.9874\left(34\right)$	$0.9378\left(38\right)$	$1.3173\left(47\right)$
$𝒪\left(p^6\right)$	$1.0752\left(30\right)$	$1.0248\left(32\right)$	$1.0014\left(33\right)$	$0.9636\left(35\right)$	$0.9169\left(37\right)$	$1.4451\left(84\right)$
Fat$3$	$1.0118\left(27\right)$	$1.0264\left(27\right)$	$1.0337\left(27\right)$	$1.0449\left(27\right)$	$1.0761\left(27\right)$	$1.1309\left(28\right)$
$a{q}^{\ast }\left(\delta C\right)$
$𝒪\left(p^4\right)$	$1.4998\left(91\right)$	$1.527\left(91\right)$	$1.589\left(74\right)$	$1.63\left(12\right)$	$2.508\left(45\right)$	$2.56\left(11\right)$
$𝒪\left(p^6\right)$	$1.5868\left(48\right)$	$1.63\left(12\right)$	$1.69\left(11\right)$	$1.79\left(11\right)$	$2.428\left(21\right)$	$2.3549\left(29\right)$
Fat$3$	$1.866\left(26\right)$	$1.390\left(76\right)$	$0.7903\left(92\right)$	$2.89\left(34\right)$	$1.699\left(49\right)$	$1.470\left(18\right)$

Table 8. Table 8: Numerical values of the one-loop radiative shift relevant for non-perturbative calculations. Note that the unsmeared results ( 𝒪 ​ ( p 4 ) 𝒪 superscript 𝑝 4 \mathcal{O}(p^{4}) and 𝒪 ​ ( p 6 ) 𝒪 superscript 𝑝 6 \mathcal{O}(p^{6}) ) are mean-field improved given the formulae in Appendix B and the mean field parameter in Table 1 . The smeared results (Fat3) are not mean-field improved. This data is plotted in Figure 6 . To determine the physical scale, q ∗ superscript 𝑞 q^{*} , we use a − 1 = 1.3 superscript 𝑎 1 1.3 a^{-1}=1.3 , 1.6 1.6 1.6 , 2.2 2.2 2.2 and 3.3 3.3 3.3 GeV corresponding to very coarse, coarse, fine and superfine MILC ensembles used by the HPQCD collaboration Dowdall:Upsilon . The a ​ m b 𝑎 subscript 𝑚 𝑏 am_{b} values of 3.4, 2.8, 1.9 and 1.1 are the appropriate ones (approximately) for the b 𝑏 b quark on very coarse, coarse, fine and superfine ensembles respectively. However, here we give results for all 4 lattice masses at each lattice spacing for completeness.

${\alpha }_s\left(q^{\ast }\right){\tilde{c}}_1^{\left(1\right)}$
$a{m}_b$	$3.4$	$2.8$	$1.9$	$1.1$	$3.4$	$2.8$	$1.9$	$1.1$
Very Coarse					Coarse
$𝒪\left(p^4\right)$	$0.2138\left(39\right)$	$0.1983\left(38\right)$	$0.1713\left(37\right)$	$0.1094\left(26\right)$	$0.1905\left(29\right)$	$0.1756\left(28\right)$	$0.1494\left(27\right)$	$0.0944\left(18\right)$
$𝒪\left(p^6\right)$	$0.1719\left(31\right)$	$0.1704\left(32\right)$	$0.1445\left(31\right)$	$-0.4417\left(64\right)$	$0.1533\left(23\right)$	$0.1513\left(24\right)$	$0.1265\left(22\right)$	$-0.4031\left(52\right)$
Fat$3$	$0.0682\left(16\right)$	$0.0391\left(12\right)$	$-0.0682\left(14\right)$	$-0.4775\left(99\right)$	$0.0605\left(13\right)$	$0.03424\left(97\right)$	$-0.0607\left(11\right)$	$-0.4190\left(73\right)$
Fine					Superfine
$𝒪\left(p^4\right)$	$0.1641\left(21\right)$	$0.1503\left(20\right)$	$0.1260\left(18\right)$	$0.0788\left(12\right)$	$0.1403\left(15\right)$	$0.1278\left(14\right)$	$0.1059\left(12\right)$	$0.06573\left(79\right)$
$𝒪\left(p^6\right)$	$0.1323\left(17\right)$	$0.1298\left(17\right)$	$0.1071\left(15\right)$	$-0.3564\left(40\right)$	$0.1132\left(12\right)$	$0.1106\left(12\right)$	$0.0902\left(10\right)$	$-0.3117\left(31\right)$
Fat$3$	$0.05184\left(97\right)$	$0.02897\left(75\right)$	$-0.05227\left(80\right)$	$-0.3555\left(50\right)$	$0.04413\left(75\right)$	$0.02441\left(59\right)$	$-0.04466\left(60\right)$	$-0.3001\left(34\right)$

Table 9. Table 9: The same as in Table 8 but for α s ​ c 5 ( 1 ) subscript 𝛼 𝑠 superscript subscript 𝑐 5 1 \alpha_{s}c_{5}^{(1)} .

${\alpha }_s\left(q^{\ast }\right)c_5^{\left(1\right)}$
$a{m}_b$	$3.4$	$2.8$	$1.9$	$1.1$	$3.4$	$2.8$	$1.9$	$1.1$
Very Coarse					Coarse
$𝒪\left(p^4\right)$	$0.1596\left(31\right)$	$0.1549\left(28\right)$	$0.1415\left(24\right)$	$0.0912\left(15\right)$	$0.1444\left(25\right)$	$0.1402\left(23\right)$	$0.1281\left(20\right)$	$0.0824\left(13\right)$
$𝒪\left(p^6\right)$	$0.04513\left(79\right)$	$0.05182\left(85\right)$	$0.0754\left(11\right)$	$-$	$0.04055\left(62\right)$	$0.04698\left(68\right)$	$0.06879\left(86\right)$	$-$
Fat$3$	$0.1491\left(34\right)$	$0.1377\left(31\right)$	$0.0990\left(23\right)$	$-0.0434\left(40\right)$	$0.1318\left(26\right)$	$0.1214\left(24\right)$	$0.0867\left(18\right)$	$-0.0384\left(32\right)$
Fine					Superfine
$𝒪\left(p^4\right)$	$0.1265\left(20\right)$	$0.1228\left(18\right)$	$0.1122\left(15\right)$	$0.07210\left(95\right)$	$0.1098\left(16\right)$	$0.1065\left(14\right)$	$0.0973\left(11\right)$	$0.06244\left(72\right)$
$𝒪\left(p^6\right)$	$0.03526\left(45\right)$	$0.04125\left(51\right)$	$0.06085\left(66\right)$	$-$	$0.03037\left(33\right)$	$0.03584\left(38\right)$	$0.05324\left(49\right)$	$0.0915\left(39\right)$
Fat$3$	$0.1126\left(19\right)$	$0.1035\left(17\right)$	$0.0734\left(13\right)$	$-0.0328\left(25\right)$	$0.0956\left(15\right)$	$0.0877\left(13\right)$	$0.06187\left(92\right)$	$-0.0279\left(20\right)$

Equations

$$G(\mathbf{x},t+1)$$

$$G(\mathbf{x},t_{\text{src}})$$

$$e^{-aH}$$

$$\times {\left(1-\frac{a{H}_0{|}_t}{2n}\right)}^n\left(1-\frac{a\delta H{|}_t}{2}\right),$$

$$a{H}_0$$

$$a\delta H$$

$$a\delta H_{v^4}$$

$$-c_3\frac{1}{8(a{m}_b{)}^2}\mathbf{\sigma }\cdot \left(\tilde{\nabla }\times \tilde{\mathbf{E}}-\tilde{\mathbf{E}}\times \tilde{\nabla }\right)$$

$$-c_4\frac{1}{2a{m}_b}\mathbf{\sigma }\cdot \tilde{\mathbf{B}}+c_5\frac{{\Delta }^{(4)}}{24a{m}_b}-c_6\frac{({\Delta }^{(2)}{)}^2}{16n(a{m}_b{)}^2},$$

$$\delta H_{p^6}$$

$$-\frac{c_{p^6}}{180a{m}_b}{\Delta }^{(6)}$$

$$+\frac{c_{p^2{p}^4}}{48(a{m}_b{)}^3}({\Delta }^{(2)}{\Delta }^{(4)}).$$

$$c_i$$

$$G(p{)}^{-1}$$

$$G^{(0)}(p{)}^{-1}$$

$$F(p)$$

$$F_1(p)$$

$$+\frac{{\tilde{c}}_1}{m_b^3}\left[1+\frac{m_b}{2n}\right]{\left[\sum _j{\sin }^2(p_j/2)\right]}^2$$

$$-\frac{2c_{(p^2{)}^3}}{m_b^5}\left[1-\frac{m_b^2}{6n^2}\right]{\left[\sum _j{\sin }^2(p_j/2)\right]}^3$$

$$-\frac{8c_{p^6}}{45m_b}\sum _j{\sin }^6(p_j/2)$$

$$+\frac{2c_{p^2{p}^4}}{3m_b^3}\sum _{j,k}{\sin }^2(p_j/2){\sin }^4(p_k/2)$$

$$\omega (\mathbf{p})$$

$${\omega }_0(\mathbf{p})$$

$$+{\alpha }_s\left\{c_5^{(1)}\frac{{\mathbf{p}}^4}{24m_b}-{\tilde{c}}_1^{(1)}\left(\frac{1}{m_b}+\frac{1}{2n}\right)\frac{({\mathbf{p}}^2{)}^2}{8m_b^2}\right\}.$$

$$\Sigma (p)$$

$${\Sigma }_0(\omega )=\Sigma (\mathbf{p}=0),$$

$${\Sigma }_1(\omega )=m_b{\frac{{\partial }^2\Sigma (p)}{\partial p_z^2}|}_{\mathbf{p}=0},$$

$${\Sigma }_2(\omega )=m_b^2{\frac{{\partial }^4\Sigma (p)}{\partial p_z^2\partial p_y^2}|}_{\mathbf{p}=0},$$

$${\Sigma }_3(\omega )=\frac{1}{24}{\left(\frac{{\partial }^4\Sigma (p)}{\partial p_z^4}-3\frac{{\partial }^4\Sigma (p)}{\partial p_z^2\partial p_y^2}\right)|}_{\mathbf{p}=0},$$

$${\Sigma }_m^{(l)}={(-i{)}^l\frac{1}{l!}\frac{{\partial }^l{\Sigma }_m(p_4)}{\partial p_4^l}|}_{p_4=0}.$$

$$\Sigma ({\omega }_0,\mathbf{p})=W_0+\frac{{\mathbf{p}}^2}{2m_b}{Z}_m^{(1)}$$

$$+\frac{({\mathbf{p}}^2{)}^2}{8m_b^2}\left\{W_1-\frac{3Z_m^{(1)}}{m_b}\right\}+W_2{\mathbf{p}}^4,$$

$$m_b^r=Z_m{m}_b=m_b\left(1+{\alpha }_s{Z}_m^{(1)}+\mathcal{O}({\alpha }_s^2)\right),$$

$$Z_m^{(1)}={\Sigma }_0^{(1)}+{\Sigma }_1^{(0)},$$

$$W_0={\Sigma }_0^{(0)},$$

$$W_1=2{\Sigma }_0^{(2)}+2{\Sigma }_1^{(1)}+{\Sigma }_2^{(0)}+\frac{2{\Sigma }_0^{(1)}}{m_b}+\frac{3{\Sigma }_1^{(0)}}{m_b},$$

$$W_2={\Sigma }_3^{(0)}$$

$$\omega (\mathbf{p})=\frac{{\mathbf{p}}^2}{2m_b^r}-\frac{({\mathbf{p}}^2{)}^2}{8(m_b^r{)}^3}+\frac{({\mathbf{p}}^2{)}^3}{16(m_b^r{)}^5}$$

$$\displaystyle\hskip 28.45274pt-\alpha_{s}\bigg{\{}W_{0}+{{\bf\bm{p}}}^{4}\left[W_{2}-\frac{c_{5}^{(1)}}{24m_{b}}\right]$$

$$\displaystyle\hskip 56.9055pt+\frac{({{\bf\bm{p}}}^{2})^{2}}{8m_{b}^{2}}\left[\left(\frac{1}{m_{b}}+\frac{1}{2n}\right)\tilde{c}_{1}^{(1)}+W_{1}\right]\bigg{\}}.$$

$${\tilde{c}}_1^{(1)}=-{\left(\frac{1}{m_b}+\frac{1}{2n}\right)}^{-1}{W}_1,$$

$$c_5^{(1)}=24m_b{W}_2,$$

$$C={\omega }^{(QCD)}-\omega =m_b^r+{\alpha }_s{W}_0=m_b(1+{\alpha }_s\delta C),$$

$$\delta C=Z_m^{(1)}+\frac{W_0}{m_b}.$$

$$\log (q^{\ast 2})$$

Full text

††preprint: FERMILAB-PUB-18-698-T

I Introduction

The Standard Model (SM) of particle physics has been incredibly successful at describing experimental data to date PDG:2016 ; QWG:2010 . However, in many ways, this success has been a double-edged sword; while SM predictions have overwhelmingly agreed with experimental measurements within errors, this has left little room for large new-physics effects to be observed. Consequently, to illuminate any new-physics phenomena high-precision tests of the SM must be performed. In the $b$-quark sector, the LHCb and BELLE II experiments will generate increasingly precise measurements. In response to this, we make the next level of improvement to the HPQCD collaboration’s formulation of the NRQCD action Dowdall:Upsilon which has been used for a number of state-of-the-art $b$-physics calculations Daldrop:Dwave ; Dowdall:Hl ; Dowdall:BMeson ; Brian:LeptonicWidth ; Hughes:Hindered ; Dowdall:Hyperfine ; Dowdall:ErratumHF ; Hughes:tetra .

In this study we will include, for the first time, operators in the NRQCD action which reproduce the correct quark dispersion relation to ${\mathcal{O}(p^{6})}$. Then, with different values of the NRQCD stability parameter $n$, we use lattice perturbation theory to compute the kinetic matching coefficients to ${\mathcal{O}(\alpha_{s}p^{4})}$. We remove the unphysical tadpole contributions Lepage:Pert from the lattice action and give perturbative results for two different tadpole improvement programs: the first by using a mean-field improvement parameter in Landau gauge $u_{0}$ Lepage:Pert , and the second via Fat$3$ smearing Orginos:Fat3 . Additionally, we determine the one-loop (bare-to-pole) mass renormalisation and zero-point energy of the $b$-quark. These can be combined to give the one-loop energy shift of the NRQCD heavy quark, and added to non-perturbatively obtained static masses to give numerical results which, after converting from lattice units to GeV, can be compared to experimental data. Further, for each of these quantities the scale $\mu=q^{*}$ at which to evaluate the strong coupling constant defined in the V-scheme is determined using the Brodsky-Lepage-Mackenzie (BLM) procedure Lepage:Pert ; Hornbostel:Scale .

After perturbatively determining the full one-loop radiative corrections to the kinetic couplings, we non-perturbatively determine the $\Upsilon(1S)$ and $\eta_{b}(1S)$ energies in order to examine how improving the NRQCD action, both with additional ${\mathcal{O}(p^{6})}$ operators and with the ${\mathcal{O}(\alpha_{s}p^{4})}$ couplings, reduces the effect of lattice artefacts.

This paper is organised as follows. In Section II we describe the improved NRQCD action. In Section III we match the ${\mathcal{O}(\alpha_{s}p^{4},p^{6})}$ NRQCD dispersion relation to the continuum, describe our tadpole improvement procedures and how the scale at which to evaluate $\alpha_{V}$ is found. Section III.1 describes the computational setup of the automated lattice perturbation theory, while Section III.2 gives an analysis of the perturbative results. Section IV gives details of the non-perturbative calculation and Section IV.2 presents the non-perturbative results. We summarise our findings in Section V.

II $b$-Quarks Using NRQCD

Information about processes involving heavy quarks can be computed on the lattice using correlation functions constructed from combinations of heavy-quark propagators. Current lattice ensembles have small enough lattice spacings and large enough volumes so that accurate relativistic $c$-quark formalisms (e.g., Highly Improved Staggered Quarks (HISQ) HISQAction ) are available. Since the $b$-quark has a Compton wavelength of $\mathcal{O}(0.04)$ fm, most current lattice ensembles cannot resolve relativistic $b$-quarks since $am_{b}>1$ PhysRevD.93.094510 111Combining results at multiple lattice spacing values and multiple heavy quark masses with a highly improved relativistic action does allow results to be obtained at the physical $b$ quark mass PhysRevD.86.074503 .. However, it is well known that $b$-quarks are very nonrelativistic inside their bound states (with $v_{\text{rel}}^{2}\approx 0.1$ for low-lying bottomonium states) and thus using a nonrelativistic effective field theory, which has a formal expansion in $p/m_{b}=v_{\text{rel}}$ Lepage:ImprovedNRQCD , is very appropriate. This effective field theory is then discretised as lattice NRQCD Lepage:ImprovedNRQCD .

HPQCD’s formulation of lattice NRQCD has already proven successful in producing accurate $b$-physics results in the literature. For example, the NRQCD formalism that gave a quark dispersion relation correct to ${\mathcal{O}(\alpha_{s}p^{4})}$ has already been used to study bottomonium $S$, $P$ and $D$ wave mass splittings Dowdall:Upsilon ; Daldrop:Dwave , $B$ meson mass splittings Dowdall:Hl , $B$ meson decay constants Dowdall:BMeson ; Hughes:2017spc , $\Upsilon$ and $\Upsilon^{\prime}$ leptonic widths Brian:LeptonicWidth . Subsequently, the spin-dependent ${\mathcal{O}(v^{6})}$ operators were added to that NRQCD action in order to compute hindered M$1$ radiative decays Hughes:Hindered , precise bottomonium hyperfine splittings Dowdall:Hyperfine ; Dowdall:ErratumHF 222Four-quark operators were also used in this study. and to aid in the search for $bb\bar{b}\bar{b}$-type bound tetraquarks Hughes:tetra .

Given the increasingly important emphasis being put on high-precision calculations needed to keep pace with measurements from the LHCb and BELLE II experiments, we take the next steps in improving the lattice NRQCD action to reduce the systematic uncertainties in future theoretical calculations using it. The first part of this improvement is to add the necessary operators to the aforementioned NRQCD action that reproduce the correct ${\mathcal{O}(p^{6})}$ quark dispersion relation at tree level.

The NRQCD action that gives rise to a ${\mathcal{O}(p^{6})}$ correct quark dispersion relation, including ${\mathcal{O}(v^{4})}$ interaction operators Dowdall:Upsilon , produces a heavy-quark propagator which can be found through the evolution equation

[TABLE]

where $\phi(x)$ is a source function and

[TABLE]

[TABLE]

Here, $am_{b}$ is the bare $b$-quark mass, $\nabla$ is the symmetric lattice derivative, with $\tilde{\nabla}$ the improved version, and $\Delta^{(2)}$, $\Delta^{(4)}$, $\Delta^{(6)}$ are the lattice discretisations of $\sum\limits_{i}D_{i}^{2}$, $\sum\limits_{i}D_{i}^{4}$ and $\sum\limits_{i}D_{i}^{6}$ respectively, with our conventions given in Appendix A. ${{\bf\bm{\tilde{E}}}}$, ${{\bf\bm{\tilde{B}}}}$ are the improved chromoelectric and chromomagnetic fields, details of which can be found in Dowdall:Upsilon . Each of these fields, as well as the covariant derivatives, must be tadpole-improved using the same improvement procedure as in the perturbative calculation of the matching coefficients Lepage:Pert . This will be discussed further in Sec. II.2. The parameter $n$ is used to prevent instabilities at large momentum from the kinetic energy operator, and needs to satisfy the constraint $p^{2}<4nam_{b}$. A choice of $n=4$ was suitable for values of $am_{b}$ used in previous non-perturbative studies. We choose to put the $a\delta H_{p^{6}}$ corrections into $a\delta H$ rather than alter $aH_{0}$ so that the kinetic operator remains unchanged. This formulation is also consistent with previous HPQCD NRQCD actions, is symmetric with respect to time reversal and has smaller renormalisations than other formulations Lepage:ImprovedNRQCD . The rotationally-symmetry breaking operators (which vanish as $a\to 0$) with coefficients $c_{p^{6}}$ and $c_{p^{2}p^{4}}$ in $\delta H_{p^{6}}$ remove higher-order discretisation effects from using finite-difference derivatives. The operator with coefficient $c_{(p^{2})^{3}}$ correctly adds the term proportional to $({{\bf\bm{p}}}^{2})^{3}$ into the heavy-quark dispersion relation.

The matching coefficients $c_{i}$ in the above Hamiltonian take into account the high-energy UV modes from QCD processes that are not present in NRQCD. Each $c_{i}$ can be fixed by matching a particular lattice NRQCD formalism to full continuum QCD. Each $c_{i}$ can be expanded perturbatively as

[TABLE]

and, after tadpole improvement Lepage:Pert , we expect $c_{i}^{(1)}$ to be $\mathcal{O}(1)$. In Sec. II.1, we will match the on-shell NRQCD dispersion relation to that of the continuum, and determine $c_{1}^{(1)}$, $c_{6}^{(1)}$ and $c_{5}^{(1)}$. Each of these coefficients should exhibit benign behaviour as a function of $am_{b}$ in the regime where the NRQCD effective field theory is well-behaved. In contrast, the coefficient may diverge as the effective field theory breaks down as $p\sim\pi/a$ gets too large or $am_{b}$ gets too small. We take tree level values, $c_{i}=1$, for the coefficients appearing in $\delta H_{p^{6}}$.

We call the NRQCD Hamiltonian presented in Eq. (3) the ${\mathcal{O}(p^{6})}$ Hamiltonian, while choosing $a\delta H_{p^{6}}=0$ produces the ${\mathcal{O}(p^{4})}$ Hamiltonian. When including the one-loop corrections to $c_{1}$, $c_{6}$ and $c_{5}$, we denote the ${\mathcal{O}(p^{6})}$ NRQCD action as being ${\mathcal{O}(\alpha_{s}p^{4},p^{6})}$, while if $\delta H_{p^{6}}=0$ then the action is ${\mathcal{O}(\alpha_{s}p^{4})}$.

II.1 One-loop Matching to ${\mathcal{O}(\alpha_{s}p^{4},p^{6})}$

A high-precision non-perturbative calculation of mass splittings will require knowledge of at least the $\mathcal{O}(\alpha_{s})$ corrections to the matching coefficients in order to improve upon existing few percent errors. For example, when tuning the bare quark mass $am_{b}$ fully nonperturbatively in NRQCD, one computes the kinetic mass of a hadron333The static mass (the energy corresponding to zero-spatial momentum) in lattice NRQCD Dowdall:Upsilon is shifted due to the removal of the mass term from the Hamiltonian and so one can only tune static mass differences fully nonperturbatively. Dowdall:Upsilon . This kinetic mass depends on the internal kinematics of the hadron, and hence on (at least) the terms $c_{1}$, $c_{5}$, and $c_{6}$ in the Hamiltonian. These matching coefficients are known as the kinetic couplings Morningstar:KinCouplings .

The kinetic couplings can be found perturbatively by matching the NRQCD on-shell energy (which corresponds to the location of the pole of the quark propagator in the interacting theory) to the continuum QCD dispersion relation. From now on, to avoid superfluous notation, we will implicitly work in lattice units unless otherwise stated. To ${\mathcal{O}(\alpha_{s})}$ the inverse quark propagator may be written in momentum space as

[TABLE]

with $G^{(0)}(p)^{-1}$ the quark propagator obtained at tree-level from the non-interacting part of the NRQCD action, $\Sigma(p)$ the one-loop quark self-energy, $p=(p_{4},{{\bf\bm{p}}})$ a four-vector in Euclidean space and $\omega=-ip_{4}$ the energy in Minkowski space. The free quark propagator can be explicitly found as

[TABLE]

[TABLE]

where we have defined $\tilde{c}_{1}=(c_{1}+c_{6}m_{b}/2n)/(1+m_{b}/2n)$ for computational ease, and $c_{1}=c_{6}=\tilde{c}_{1}$. The term $F(p)$ arises from the non-interacting momentum space part of $(1-H_{0}/2n)$ in (2), while $F_{1}(p)$ comes from the $(1-\delta H/2)$ piece.

To find the NRQCD dispersion relation we determine the on-shell energy $\omega({{\bf\bm{p}}})$ which causes a pole in the full heavy-quark propagator. The one-loop $\omega({{\bf\bm{p}}})$ can be found from Eq. (5) and (6) as

[TABLE]

with $\omega_{0}({{\bf\bm{p}}})=-\log(F^{2n}({{\bf\bm{p}}})F_{1}^{2}({{\bf\bm{p}}}))$, with tree-level coefficients in $F$ and $F_{1}$, being the tree-level on-shell energy found by setting the tree-level inverse propagator in Eq. (6) to zero. We have constructed the Hamiltonian in Eq. (2) to produce a non-relativistic dispersion relation correct to ${\mathcal{O}(p^{6})}$, and we now include the ${\mathcal{O}(\alpha_{s}p^{4})}$ correction. This yields444We correct a typographical error in Appendix B of Dowdall:Upsilon .

[TABLE]

When matching the dispersion relation to ${\mathcal{O}(\alpha_{s}p^{4},p^{6})}$, it is necessary to decompose the self-energy $\Sigma(p)$ using the small-${{\bf\bm{p}}}$ expansion Morningstar:KinCouplings as

[TABLE]

Further, when $\omega$ is small, each function has a well-defined series expansion $\Sigma_{m}(\omega)=\sum_{l=0}^{\infty}\Sigma_{m}^{(l)}\omega^{l}$. The $\Sigma_{m}^{(l)}$ can be found from derivatives of the quark self-energy as

[TABLE]

Then, by using the tree-level $\omega_{0}$ from (10) in Eq. (11) we find

[TABLE]

where the superscript ’r’ denotes renormalised quantities and $Z_{m}^{(1)}$ is the $\mathcal{O}(\alpha_{s})$ coefficient of the bare-to-pole mass renormalisation. Substituting (8) and (17) into (9) gives the one-loop NRQCD dispersion relation to ${\mathcal{O}(\alpha_{s}p^{4},p^{6})}$ as

[TABLE]

Matching Eq. (23) to the continuum QCD dispersion relation MovingNRQCD ; Dowdall:Upsilon gives the matching coefficients for $\tilde{c}_{1}^{(1)},c_{5}^{(1)}$ as well as the energy shift of a heavy quark (to this order) as

[TABLE]

The shift $C$ is the perturbative shift of the zero of energy. For each heavy quark in a non-perturbative calculation, the shift can be added to the simulation energy and, after being converted from lattice units to GeV, this can then be compared to experimental masses MovingNRQCD ; Morningstar:KinCouplings . In practice hadron masses can be more precisely determined fully non-perturbatively through their kinetic mass in lattice QCD.

The aim of this study is to determine the one-loop coefficients $\tilde{c}^{(1)}_{1}$, $c^{(1)}_{5}$, $\delta C$ (and thus also $Z^{(1)}_{m}$ and $W_{0}$) for different improved NRQCD actions to find the best way forward for increasingly accurate non-perturbative calculations in the future. Before these coefficients can be used, it is first necessary to remove unphysical contributions from tadpole diagrams which can cause the coefficients to be rather large Lepage:Pert .

II.2 Tadpole Improvement

The authors of Ref. Lepage:Pert show that using Lie group elements when constructing the lattice field theory introduces unphysical tadpole diagrams which do not contribute to continuum schemes. These unphysical tadpole diagrams cause large, process independent renormalisations and produce a poor convergence of the perturbative series. Ref. Lepage:Pert also suggests a solution to this: a gauge-invariant mean-field improvement program (tadpole-improvement) where each lattice link, $U_{\mu}(x)$, is scaled to $U_{\mu}(x)/u_{0}$. We choose $u_{0}$ to be the mean link in Landau gauge, i.e., $u_{0}=\langle\frac{1}{3}\text{Tr}U_{\mu}(x)\rangle$. This mean-field parameter has been calculated for the Symanzik-improved gluon Dowdall:Upsilon ; Hammant:2013 action both perturbatively to one-loop (with $u_{0}=1-\alpha_{s}u_{0}^{(2)}$ giving $u_{0}^{(2)}=0.750$) Nobes:2001 and non-perturbatively Dowdall:Upsilon ; Dowdall:BMeson (where the value of $u_{0}$ depends on the ensemble used, e.g., see Table 4). After a tadpole improvement procedure has been implemented, the one-loop coefficients are expected to be ${\mathcal{O}(1)}$. The same tadpole-improvement program must of course be implemented in the nonperturbative calculations as has been used for the perturbative calculations.

Before the mean-field improvement procedure is performed, care must be taken to ensure that any link-pair cancellations $U^{\dagger}_{\mu}(x)U_{\mu}(x)=1$ occur in the lattice action used in both non-perturbative and perturbative calculations. Such cancellations do not generate any unphysical tadpole diagrams and scaling by $1/u_{0}^{2}$ would be incorrect. Yet, expanding out the complicated NRQCD Hamiltonian in (2) in terms of links $U_{\mu}(x)$ is excessively expensive for numerical calculations. Consequently, link-pair cancellations are only taken into account separately for each derivative, $(\Delta^{(2n)})^{m}$, or field strengths ($E_{i}(x)$ or $B_{i}(x)$) appearing in the action. This is called partial cancellation MovingNRQCD ; Dowdall:Upsilon . The difference between the complete and partial cancellation prescriptions was empirically shown not to be sizable MovingNRQCD . Formulae for the partially-cancelled derivatives are given in Appendix A.

By using the partially-cancelled mean-field improvement procedure just described, one can find the ${\mathcal{O}(\alpha_{s})}$ tadpole counterterms for the one-loop quantities described in Sec. II.1. For the NRQCD action without the ${\mathcal{O}(p^{6})}$ operators, the computation of the tadpole counterterms was checked in two separate calculations. The first was performed analytically, and the second using a Mathematica script. Both calculations reproduced the results of MovingNRQCD (in the case where the additional parameter used there, $v$, is set to zero) and Dowdall:Upsilon ; Gulez:A0V0 . We extended the numerical code to include the ${\mathcal{O}(p^{6})}$ operators. The one-loop tadpole counterterms are given in Appendix B. After the (unimproved) one-loop quantities have been found in lattice perturbation theory, we can add the appropriate tadpole counterterms to determine the improved values. This will be discussed further in Sec. III.2.

As can be seen in Appendix B, the mean-field counterterms obtained from using the ${\mathcal{O}(p^{6})}$ NRQCD action contain higher-orders of $1/m_{b}$ relative to the counterterms obtained from using the ${\mathcal{O}(p^{4})}$ NRQCD action. This is a consequence of partially-cancelling the derivative operators $(\Delta^{(2n)})^{m}$, whose counterterms are given in Appendix A. The impact of this becomes pronounced as $am_{b}$ is reduced as will become evident in Section III.

In this study, we choose to account for the unphysical tadpole contributions using two different prescriptions. The first prescription proceeds via the partially-cancelled mean-field improvement procedure just described. The one-loop tadpole counterterms given in Appendix B, which depend on $1/am_{b}$, remove the unphysical tadpoles. As seen in Sec. III.2, the improved values give smaller absolute renormalisations compared to the unimproved case, and exhibit a longer plateau over a larger range in $am_{b}$ indicating stable behaviour in the effective field theory. However, the tadpole counterterms from using the ${\mathcal{O}(p^{6})}$ action diverge faster as $am_{b}\to 0$ due to the higher-order terms in $1/am_{b}$, and therefore the tadpole-improved one-loop results obtained from the ${\mathcal{O}(p^{6})}$ NRQCD action also diverge faster (see Sec. III.2). This could be slightly inconvenient for ensembles with increasingly small lattice spacings, such as the super-fine ensembles currently in use Fermilab:Vub , which have a lattice spacing of $a\approx 0.06\mathrm{fm}$.

Because of this, we explore an alternative improvement procedure based on the fattening of gauge-links Orginos:Fat3 ; HIPPY . The Fat$7$-smeared link Orginos:Fat7 introduces staples of up to seven-link paths to completely remove the tree-level couplings to gluons with high transverse-momentum modes equal to $\pm\pi$. As the Fat$7$ link is computationally expensive, alternative fat-links have been designed based on three- or five-link staples, called Fat$3$ and Fat$5$ respectively. These latter links reduce the couplings to gluons with high transverse-momentum and suppress unphysical tadpole diagrams Orginos:Fat7 . This will be discussed further in section III.2. Therefore, we also explore, for the first time, how a Fat$3$-smeared NRQCD action correct to ${\mathcal{O}(p^{4})}$ affects the renormalisation of kinetic couplings. Here the fattened links are projected back onto $U(3)$ HIPPY (not $SU(3)$).

The last piece of information needed to use the tadpole-improved one-loop coefficients in a non-perturbative computation is the scale, $q^{*}$, at which to evaluate the strong coupling constant.

II.3 Determining the scale of $\alpha_{s}$

The Brodsky-Lepage-Mackenzie procedure BLM ; Lepage:Pert determines an optimal $q^{*}$ for $\alpha_{V}$, the coupling defined using the heavy quark potential PhysRevD.28.228 , by examining the momentum flowing through a gluon in the one-loop Feynman diagram. In this prescription, one studies the one-loop integral of a fully dressed gluon within a particular diagram, then uses the running of $\alpha_{V}(q)$ to find a mean-value $q^{*}$ which reproduces the integral. To do this, one expands the running of $\alpha_{V}(q)$ as a polynomial in $\log(q^{2}/q^{*2})$ and assumes that the leading order log-moments are the dominant contributions. However, in certain areas of parameter space, the leading order log-moments can be anomalously small and give unphysically large or small erroneous $q^{*}$. This was noticed in Morningstar:KinCouplings after which Hornbostel:Scale determined $q^{*}$ when the zeroth and first log-moments are anomalously small via

[TABLE]

where $\langle\langle\log(q^{2})\rangle\rangle=\langle f(q)\log(q^{2})\rangle/\langle f(q)\rangle$ indicates the weighted average, $f(q)$ is the integrand of the one-loop Feynman diagram, and $\sigma^{2}=\langle\langle\log^{2}(q^{2})\rangle\rangle-\langle\langle\log(q^{2})\rangle\rangle^{2}$. The appropriate choice of $\pm$ in Eq. (28) is usually clear based on requiring $q^{*}$ to be continuous and physically sensible, although the ambiguity can be removed by calculating higher log-moments Hornbostel:Scale . When $\sigma^{2}>0$, only the first term in Eq. (28) is used. However, when $\sigma^{2}<0$, Eq. (28) takes into account the anomalies to first order.

The unphysical tadpole diagrams contribute to the scale $q^{*}$, using the mean-field improvement prescription described in Sec. II.2 will alter its value. When the tadpole counterterm $c_{tad}\alpha_{V}(q^{*}_{tad})$ is added to the one-loop contribution $c_{a}\alpha_{V}(q^{*}_{a})$, the second-order formula (28) is altered to Hornbostel:Scale

[TABLE]

Again, if $\sigma^{2}>0$, then only the first-order term in Eq. (29) is needed and used, while if $\sigma^{2}<0$ both terms are needed to yield physical results.

Theoretically, we expect $a\Lambda_{QCD}<aq^{*}<\pi$ as the one-loop corrections take into account UV modes neglected by imposing a momentum cutoff. Even though the corrected second order formula given in Eq. (29) was used, unphysical values of $q^{*}$ for certain values of $am_{b}$ in the one-loop quantities were sometimes obtained. In these cases, although rare, it was usually clear that the issue was due to the $0^{th}-2^{nd}$ log-moments being anomalously small. To rectify this issue, we use the simple $n^{th}$-order formula given by Hornbostel:Scale

[TABLE]

Leaving the tadpole pieces out of Eq. (31) gives the higher-order tadpole-unimproved scale. As we do not mean-field improve the Fat$3$-smeared one-loop quantities, the above formulae with the tadpole pieces set to zero are used to find $q^{*}$ in the case of Fat3-smeared links.

III Perturbative Determination of One-loop Quantities

III.1 Perturbative Computational Details

Due to the complexity of the NRQCD action that we utilise, an efficient computational methodology is needed to calculate the Feynman integrals of the one-loop formulae given in (21) and (22). Fortunately, the theory behind the automatic generation of Feynman rules for complex lattice actions exists Luscher:Pert ; HIPPY . Here, we employ the automated lattice perturbation theory routines HiPPy and HPsrc HIPPY ; HIPPY2 . These routines have been thoroughly tested and used in previous perturbative calculations Dowdall:Hl ; Dowdall:Upsilon ; Horgan:Tads ; MovingNRQCD ; Hammant:2013 . Given that we will produce results for a number of different NRQCD actions, these automated packages are ideal.

We automatically generate the Feynman rules for a specific NRQCD action (along with the Symanzik-improved gluonic action Horgan:Tads ; Dowdall:Upsilon ) using the HiPPy package. We can then construct the Feynman diagrams in a generic fashion using the HPsrc package, which will use these Feynman rules to numerically evaluate the diagram, along with its derivatives thanks to automated differentiation techniques vonHippel:Taylur ; HIPPY .

In the matching procedure both the continuum and lattice contributions to the dispersion relation are separately infrared (IR) finite. However, intermediate steps on the lattice may produce IR divergences which cancel when evaluating the one-loop quantities. To regulate the IR divergences, we use twisted boundary conditions (TBCs) on a finite-volume lattice where the momentum integral is replaced by a summation over momentum modes Luscher:Pert . TBCs introduce a lower momentum cutoff by removing the zero mode from the gluon propagator. Specifically, we employ triple-twist boundary conditions with an appropriate squashing factor in the untwisted temporal direction (used to broaden peaks of the integrand, thus removing most of the dependence on $L$ Luscher:Pert ). Computational details of both concepts are described in Horgan:Tads ; HIPPY , and we refer the reader to those articles for further details. All numerical results are IR finite as expected. As the dispersion relation is UV finite, this allows us to directly equate results obtained on the lattice to those obtained in the continuum. Furthermore, we test that the gauge-invariant quantities are independent of the gluon propagator gauge parameter by working in both Feynman gauge and Landau gauge. All perturbative results presented, except for the Landau-gauge mean-field parameter $u_{0}^{(2)}$, will be in Feynman gauge.

The one-loop contributions to the self-energy are shown in Figure 1. Care must be taken when numerically evaluating the rainbow diagram so that the pole of the heavy-quark propagator does not cross the temporal integration contour. Details of our implementation of the contour shift can be found in MovingNRQCD .

In this study we will always take the spatial box length to be $4\leq L\leq 16$ and choose a temporal extent of $T=16L$. This allows the pole structure to be resolved in greater detail and reduces finite-$T$ effects. As the one-loop integration is carried out by direct summation of the twisted momentum modes, numerical results are exact Horgan:Tads . We follow the approach suggested by Lepage:Code in order to fit exact data. Here, our exact results from TBCs can be expressed as a polynomial in $1/L$ Luscher:Pert , yet we are only interested in knowing the constant term (corresponding to the infinite-volume result). We may use priors to model the polynomial dependence and then marginalise MarjFit from the exact data the part of the polynomial that we are not interested in. Using a finite-degree polynomial of order $N$ to model the exact results, we find that $N=20$ is a suitable choice and check that all results are unchanged with its variation. Marginalising the last $N-N_{L}$ terms of this polynomial into the exact data and then performing a Bayesian fit Lepage:Code ; Lepage:Fitting to a polynomial of degree $N_{L}$ successfully determines the desired constant parameter of the polynomial. As is common with marginalised Bayesian fits MarjFit , marginalising all but one or two fit parameters produces stable and precise results and has seen wide success FastFits ; Hughes:Hindered . Even though we produce successful fits when marginalising all but the constant term, we choose $N_{L}=5$ and ensure that there is no sensitivity to this.

In the following section, we will give perturbative results for three NRQCD actions: (i) the ${\mathcal{O}(p^{6})}$ NRQCD action with stability parameter $n=4$ as described in Sec. II; (ii) the ${\mathcal{O}(p^{4})}$ NRQCD action with stability parameter $n=4,6$ and $8$; and (iii) a Fat$3$-smeared ${\mathcal{O}(p^{4})}$ NRQCD action with stability parameter $n=4$ and no mean-field improvement. For a fixed quark and gauge action the one-loop coefficients depend only on the input parameter $am_{b}$. We calculate results for a range of $am_{b}$ values, enabling interpolation to values not explicitly calculated that may be useful for lattice calculations. This also allows us to demonstrate the functional dependence on $am_{b}$ graphically to see where the divergent behaviour begins as $am_{b}$ goes to zero.

III.2 Perturbative Results and Analysis

We calculate $W_{1}$ and $W_{2}$ for both the ${\mathcal{O}(p^{4})}$ and ${\mathcal{O}(p^{6})}$ actions with $n=4$ (including all log moments) with spatial extent $L=4,6,8,10,16$. We then successfully fit this data using the methodology described in Sec. III.1. Figure 2 shows an example of this, with the raw $W_{1}$ data at multiple values of $1/L$ overlaid against the fit curve. In fact, we found Bayesian fitting to a polynomial so successful that we only needed data with $L=4,6,8,10$ to obtain the constant term to sub-percent precision in general. Consequently, we calculate data for $W_{1}$ and $W_{2}$ for an ${\mathcal{O}(p^{4})}$ action with $n=6,8$ and with $L=4,6,8,10$, as well as all $Z_{m}^{(1)}$, $W_{0}$ and Fat3-smeared results (including log moments). The short computational time needed to calculate $u_{0}^{(2)}$ (and its log moments), meant that we were able to do this on lattices of size $L=4,6,8,10,12,14,16$. We present the infinite-volume results for the mean-field unimproved $W_{1}$ and $W_{2}$ in Figure 3. Also shown on each figure is a smooth interpolating curve between the results. This interpolating curve was chosen to be a polynomial in $1/am_{b}$ in order to reproduce the static limit as $m_{b}\to\infty$. It is expected that all one-loop quantities diverge as $am_{b}\to 0$ for our improved NRQCD action, indicating a breakdown of NRQCD, and that is clearly illustrated in our figures.

The difference between $W_{1}$ and $W_{2}$ in Figure 3 and $\tilde{c}^{(1)}_{1}$ and $c^{(1)}_{5}$ in Figure 4 is purely the conversion factors given in Equations (24) and (25). In these plots, a clear observation is that the results are very insensitive to an increase in $n$ over the mass ranges that we are interested in for non-perturbative calculations on the lattice ($1<am_{b}<5$). Therefore, as there is no clear benefit to increase $n$ in the perturbative results, future non-perturbative calculations can choose $n=4$ for all $am_{b}$ in this range. The Fat$3$ smearing works as expected to remove the unphysical tadpoles (as outlined in Sec. II.2), indicated by a reduction in the absolute size of the one-loop corrections. There is a significant improvement, in terms of longer plateau in $am_{b}$ and sharper divergence at smaller $am_{b}$, when using the ${\mathcal{O}(p^{6})}$ NRQCD action over the ${\mathcal{O}(p^{4})}$. This improved behaviour in the couplings leads to the expectation that the second-order couplings are also well-behaved.

In Figure 5, we then include the mean-field tadpole corrections for all results (except those for Fat$3$ smearing data) with the formulae explicitly given in Appendix B. The infinite-volume values are given in Table 1. Table 1 shows why we do not need tadpole-improvement when smeared links are used. The one-loop coefficient in $u_{0}$ is much smaller in the smeared cases reflecting the fact that tadpole effects are much smaller and the mean smeared link is much closer to $1$. Little is then gained by dividing by it.

As can be seen in Figure 5, mean-field improvement noticeably reduces the magnitude of the one-loop coefficients, where it is applied. Due to the higher-order $1/am_{b}^{n}$ terms in the ${\mathcal{O}(p^{6})}$ mean-field counterterms, as described in Sec. II.2, the one-loop couplings with the ${\mathcal{O}(p^{6})}$ action now diverge earlier as $am_{b}\to 0$. This is not a desirable feature. This common behaviour is seen in all mean-field improved data we present. Interestingly, the absolute value of $c_{5}^{(1)}$ is significantly reduced when using a ${\mathcal{O}(p^{6})}$ action. $c_{5}^{(1)}$ is the coupling which removes the rotational-symmetry breaking operator $\Delta^{(4)}$ at one-loop. Therefore, it is indicative that the ${\mathcal{O}(p^{6})}$ action will reduce $SO(3)$ symmetry breaking in non-perturbative calculations also, as will be discussed in Sec. IV.2.

To fully determine the one-loop shift to the kinetic couplings, the scale at which to evaluate $\alpha_{s}$ in the V-scheme needs to be found. We give the mean-field improved $aq^{*}$ and the Fat$3$ smeared $aq^{*}$ in Appendix C. To determine the physical scale, $q^{*}$, we use $a^{-1}=1.3$, $1.6$, $2.2$ and $3.3$ GeV corresponding to very coarse, coarse, fine and superfine MILC ensembles used by the HPQCD collaboration Dowdall:Upsilon . To run the strong coupling in a particular renormalisation scheme, an initial condition needs to be chosen. Here, we use $\alpha_{s}^{\overline{\text{MS}}}(M_{Z},n_{f}=5)$ taken from the Particle Data Group PDG:2016 , where $\alpha_{s}$ is defined in the $\overline{\text{MS}}$-scheme, $M_{Z}$ is the $Z$-boson mass and $n_{f}=5$ is the number of active flavours. To use this with our data, we perturbatively remove the $b$-quarks’ contribution to the running AddQuarkSea , with $m_{b}(m_{b})=4.164(23)$ GeV MarjFit , convert to the $V$-scheme Schroder:Convert ; Lepage:Pert and run to $q^{*}$ QuarkMass:2014 . Finally, we combine $\alpha_{V}(q^{*})$ with the one-loop coefficient to give the full one-loop coefficient. These are plotted in Figure 6.

We show data for the tadpole-improved $Z_{m}^{(1)}$ in Figure 7. Without mean-field improvement, the ${\mathcal{O}(p^{4})}$ and ${\mathcal{O}(p^{6})}$ data overlap very closely and are not shown due to this. The modified behaviour after mean-field improvement is due to the different tadpole-corrections, where again, the ${\mathcal{O}(p^{6})}$ tadpole-corrections cause faster divergences as $am_{b}\to 0$. The mean-field corrections work as expected to reduce the absolute value of $Z_{m}^{(1)}$, e.g., $Z_{m}^{(1)}(am_{b}=1.9)$ is reduced from $1.57$ (without improvement) to $0.45$. The full one-loop correction is shown in the lower plot in the same figure.

$W_{0}$ is observed to have opposite sign to $Z_{m}^{(1)}$, but has similar qualitative features as those just described, e.g., the mean-field unimproved result is typically is a factor of $2-3$ in magnitude larger than the mean-field improved values and there is a clear plateau and a sharp divergence at small $am_{b}$. The tadpole-improved $W_{0}$ is shown in Figure 8.

Because $Z_{m}^{(1)}$ and $W_{0}$ have opposite sign, the one-loop shift in the zero of energy $\delta C$ is found to be very close to zero in all cases as seen in Figure 9.

We note that our ${\mathcal{O}(p^{4})}$ $\tilde{c}^{(1)}_{1}$, $c_{5}^{(1)}$ and $Z_{m}^{(1)}$ differ by small but significant amounts from those in Ref. Dowdall:Upsilon . Ref. Dowdall:Upsilon used Monte Carlo integration combined with numerical derivatives, which they note leads to unstable behaviour when there are large peaks in the IR region. Consequently subtraction functions were used Eike . In our study, we avoid these complications by using TBCs as a gauge-invariant IR regulator, and automatic differentiation for the derivatives, which avoids the numerical instabilities arising from finite-differencing schemes HIPPY ; vonHippel:Taylur .

Finally, the tadpole-improved results for the one-loop coefficients plotted here are given in Appendix C. Subtracting the mean-field corrections (given in Appendix B) from this data gives the results before tadpole improvement.

IV Non-Perturbative Kinetic masses

Here we test how improving the NRQCD action as in Sec. II and III affects the reliability and accuracy of energies of bottomonium mesons obtained from non-perturbative calculations.

The static mass (the energy corresponding to zero spatial momentum) in lattice NRQCD is shifted due to the removal of the mass term from the Hamiltonian Dowdall:Upsilon , where we found the one-loop shift, $C$, in Sec. II.1. Consequently, one can only determine static mass differences fully nonperturbatively. However, one can still obtain kinetic masses Dowdall:Upsilon ; Lewis:Rad1 entirely non-perturbatively via a fully relativistic dispersion relation as

[TABLE]

where $a\Delta E$ is the energy difference between the meson with momentum $a{{\bf\bm{P}}}$ and the meson at rest. The kinetic mass depends on the internal kinematics of the hadron, and hence on the kinetic terms in the NRQCD action. For example, changing the coefficient of the $(\Delta^{(2)})^{2}/8am_{b}^{3}$ term, $c_{1}$, from $1$ to $1+{\mathcal{O}(\alpha_{s})}$ will modify the amount of the internal kinetic energy that is incorporated into the meson’s kinetic mass, effectively correcting for an ${\mathcal{O}(\alpha_{s})}$ mismatch between the static and kinetic masses from this operator’s contribution to the binding energy Dowdall:Upsilon . The change would be expected to be ${\mathcal{O}(\alpha_{s}B)}$ where $B$ is the binding energy of ${\mathcal{O}(500)}$ MeV. This could in principle be as large as $150-200$ MeV but in practice was found to be much smaller and around $80$ MeV on coarse and fine lattices (because $c_{1}^{(1)}$ is small) Dowdall:Upsilon .

Therefore, the kinetic mass is the ideal candidate on which to test our improvement of the kinetic part of the action. Furthermore, the kinetic mass is typically utilised to tune the $b$-quark mass Dowdall:Upsilon ; Lewis:Rad1 ; Meinel ; Dowdall:bquark and thus if sizable improvement is seen, this would indicate that improving the kinetic action would benefit future calculations, where a highly accurate calculation with a reliable error budget requires knowledge of at least the ${\mathcal{O}(\alpha_{s})}$ corrections to the matching coefficients.

In a rotationally invariant theory, the symmetry group is the semi-direct product of the rotational group $SO(3)$ with three translations. The little group of the symmetry group, used to classify energy eigenstates in terms of invariant quantities (e.g., $J^{P}$ at zero-momentum and helicity $\lambda$ at non-zero momentum), is broken by a finite-volume lattice Moore:2005 . The symmetry of the lattice discretisation, which breaks $SO(3)$ symmetry at small distances, does not need to be the same as the symmetry of the finite volume, which breaks rotational symmetry at larger distances Thomas:Helicity . Here we consider a cubic lattice in a finite cubic box with PBCs, and so both the lattice and the boundary break the full $SO(3)$ rotational symmetry of the continuum to the (double cover) of the octahedral group, $O_{h}^{D}$. The lattice irreducible representations (irreps) for a cubic finite-volume on a cubic lattice depend on the allowed momenta types Moore:2005 ; Thomas:Helicity (as not all lattice-momenta are related by an octahedral symmetry) and we reproduce them in Table 2 for convenience. The energy eigenstates of the lattice Hamiltonian (as obtained from non-perturbative lattice QCD calculations) are classified according to representations of the lattice symmetry group.

We denote the energy computed on the lattice for a $\eta_{b}$ meson with spatial momentum $P$ as $E_{\eta_{b}}(|aP|)$. Then $E_{\eta_{b}}(|a{{\bf\bm{P}}}|)$ computed with the same $a^{2}{{\bf\bm{P}}}^{2}$ but with $a{{\bf\bm{P}}}$ which lie in different lattice little groups (e.g., $(3,0,0)$ which has little group $Dic_{4}$ and $(2,2,1)$ which has little group $C_{4}$ ) do not need to yield the same energy within errors. However, as the infinite-volume continuum limit is taken and full $SO(3)$ symmetry is restored, these energies should converge. Improving the lattice NRQCD action, both by adding in higher-order ${\mathcal{O}(p^{6})}$ terms and one-loop radiative corrections, should reduce $SO(3)$ symmetry breaking and produce the desired infinite-volume continuum energies more accurately at a given value of the lattice spacing. Examining the non-perturbative energies should indicate this to be the case.

Improving the NRQCD action will reduce the breaking of $SO(3)$ symmetry due to a cubic lattice. This is because higher-order rotational-symmetry breaking operators (which vanish as $a\to 0$) will be increasingly taken into account correctly, e.g., the $\sum_{i}\Delta_{i}^{4},\sum_{i,j}\Delta_{i}^{2}\Delta_{j}^{4},\sum_{i}\Delta_{i}^{6}$ operators in Eq. (2). It is perhaps indicative that including ${\mathcal{O}(p^{6})}$ operators reduces rotational symmetry breaking, as we found in Sec. III.2 that the one-loop coupling $c_{5}^{(1)}$, which is constructed in (25) to remove the rotational-symmetry breaking $\sum_{i}p_{i}^{4}$ terms from the dispersion relation to one-loop, gets reduced when improving to the ${\mathcal{O}(p^{6})}$ NRQCD action.

In the following we will describe our non-perturbative computational setup as well as discuss how the data from the kinetic masses illustrates the reduction of $SO(3)$ symmetry breaking when improving the kinetic parts of the NRQCD action.

IV.1 Non-Perturbative Computational Setup

Our computational setup is similar to that in Refs. Dowdall:Upsilon ; Hughes:Hindered and we point the reader to those texts for specific details. However, we give a brief overview. We use gauge ensembles generated by the MILC collaboration MILC:Configs with the tadpole-improved Lüscher-Weisz gauge action Hart:GluonImprovement with $2+1+1$ dynamical flavours of HISQ sea quarks HISQAction . Details of these ensembles are given in Table 3. We use ensembles at three values of the lattice spacing, approximately 0.15 fm, 0.12 fm and 0.09 fm, so that we can test the changing impact of lattice discretisation effects.

Details of the covariant derivative and chromo-magnetic/electric field implementation in our NRQCD action can be found in Dowdall:Upsilon . Each of these must be tadpole-improved using the same improvement procedure as in the perturbative calculation of the matching coefficients in Sec. III.2. We present kinetic masses using the mean-field improvement procedure where, as in the perturbative results, we take $u_{0}$ as the mean trace of the gluon field in Landau gauge, calculated in Dowdall:Upsilon ; Dowdall:BMeson . The $u_{0}$ values used for each ensemble are given in Table 4. We also give in Table 4 the values that we use for the bare $b$ quark mass $am_{b}$ on each ensemble.

The lattice two-point correlator most naturally encodes information on meson energies. We use bilinear $b\overline{b}$ interpolating operators, listed in Table 5 with $\Gamma=i\gamma_{5},\gamma^{j}$, which overlap onto definite $J^{PC}=0^{-+},1^{--}$ energy eigenstates at rest, respectively, in the infinite-volume continuum version of our theory (which is rotationally invariant) Thomas:Helicity . In Hughes:Hindered , as well as Thomas:Helicity , it has been shown that at nonzero momentum, $\mathcal{O}^{\gamma^{5}}(p)$ is a helicity operator which creates a definite $\lambda=0^{-}$ energy eigenstate, but $\mathcal{O}^{\gamma^{i}}(p)$ creates an admixture of $\lambda=0^{+},\pm 1$, where these $\lambda$ get contributions from $J^{P}$ values as listed in the third column of Table 5. The $\pm$ superscript on the $\lambda=0$ case represents the eigenvalue $\tilde{\eta}=P(-1)^{J}$ from the $\hat{\Pi}$ symmetry (a parity transformation followed by a rotation to bring the momentum direction back to the original direction) Thomas:Helicity .

Again, following Dowdall:Upsilon , we simultaneously fit multi-exponential functions to the bottomonium meson correlator at rest and with momentum $a{{\bf\bm{P}}}$. Doing so allows the correlations between the ground states energies to be correctly taken into account when computing the kinetic mass. We take priors of $0.1(1.0)$ on the amplitudes, priors on the ground state energies are estimated from previous results and given a suitably wide width Dowdall:Upsilon , and priors on energy splittings are taken to be $E_{n+1}-E_{n}=500(250)$ MeV. To help invert the covariance matrix a singular value decomposition is used with a tolerance of $10^{-5}$ ChainedFit ; Hughes:Hindered . We present fit results, following Dowdall:Upsilon , for fits including eleven exponentials for Set $1$, nine exponentials for Set $3$, and seven exponentials for Set $5$.

IV.2 Non-Perturbative Results and Analysis

We generate data for $E_{\eta_{b}}(|a{{\bf\bm{P}}}|)$ and $E_{\Upsilon}(|a{{\bf\bm{P}}}|)$ with momenta $a{{\bf\bm{P}}}=(0,0,0)$, $(1,0,0)$, $(1,1,0)$, $(1,1,1)$, $(2,0,0)$, $(2,1,1)$, $(2,2,1)$ and $(3,0,0)$ in multiples of $2\pi/L$. As discussed above, helicity classifies the energy eigenstates of the infinite-volume continuum NRQCD theory at non-zero momentum. Therefore, compared to the zero-momentum case, additional $J^{P}$ states can contribute to the correlator data at non-zero momentum. The authors of Hughes:Hindered found that when fitting to a $3\times 3$ matrix of smeared correlators, the first excited state in the fit at non-zero momentum was the $\chi_{b1}(1P)$, $h_{b}(1P)$ for the operators $\mathcal{O}^{\gamma^{5}}(x)$, $\mathcal{O}^{\gamma^{i}}(x)$ respectively. At zero momentum, the first excited state was the $\eta_{b}(2S)$, $\Upsilon(2S)$ respectively. By using the same smearing types and correlators as those authors, we check that the additional states are present at non-zero momentum when using a $3\times 3$ fit. However even when a fit does not resolve the additional (first excited) state accurately we find that the ground state is uncontaminated and still precise. Further, the finite-volume lattice breaks $SO(3)$ symmetry and allows mixing with higher $J^{P}$ states within each of the lattice irreps given in Table 2. As in Hughes:Hindered , we find no signal for any mixing in the low-lying spectrum. We conclude that our ground state energies are reliably determined.

Each $E_{\Upsilon}(|a{{\bf\bm{P}}}|)$ extracted from our lattice calculation has larger errors than those on $E_{\eta_{b}}(|a{{\bf\bm{P}}}|)$ because of the slightly poorer signal-to-noise ratio. The statistical errors also grow with momentum. Consequently, $\Delta E(|a{{\bf\bm{P}}}|)$ has larger absolute errors as $a^{2}{{\bf\bm{P}}}^{2}$ grows, but the relative error on $\Delta E(|a{{\bf\bm{P}}}|)$ is larger for small $a^{2}{{\bf\bm{P}}}^{2}$. As a result, the kinetic masses with smaller $a^{2}{{\bf\bm{P}}}^{2}$ have larger errors, which then stabilise.

On each ensemble we examine the kinetic mass, given by Eq. (32), in order to see how the kinetic mass changes for a given $am_{b}$ as a function of momentum, as we improve the NRQCD action. Since the energies and kinetic masses should only depend on the magnitude of the spatial momentum, rotational symmetry breaking effects show up most clearly as a difference between the energies corresponding to momenta $(3,0,0)$ and $(2,2,1)$ in units of $2\pi/L$, and also as a kinetic mass that depends on $a^{2}{{\bf\bm{P}}}^{2}$.

One feature of the results is that the kinetic mass for the $\eta_{b}$ is slightly larger than that of the $\Upsilon$, rather than being lower to reflect the ordering of the masses seen in experiment. This was also seen in Dowdall:Upsilon and explained there as the result of not including relativistic corrections to the $\sigma\cdot B/2m$ term in the NRQCD action (the term with coefficient $c_{4}$ in Eq. (3)). Such corrections (spin-dependent terms at $\mathcal{O}(v^{6})$) would allow the effect of the $\sigma\cdot B$ term to be correctly incorporated in the kinetic mass and solve this problem Dowdall:Hyperfine ; Dowdall:ErratumHF . The strategy adopted in Dowdall:Upsilon to mitigate this problem was to use the spin-averaged kinetic mass, which is less sensitive to these effects, to tune the $b$ quark mass. The spin-averaged kinetic mass is given by

[TABLE]

Figure 10 illustrates this feature by showing results for the $\Upsilon$ and $\eta_{b}$ kinetic masses on set 5, for the $p^{6}+\alpha_{s}p^{4}$ action. We also show the spin-averaged kinetic mass. The solid lines show the corresponding experimental values. In a full nonperturbative calculation we would want to tune the $am_{b}$ value for each action separately to match the spin-averaged kinetic mass to experiment. Here however we keep the same $am_{b}$ value for each action on a given ensemble (with only an approximate tuning) so that we can compare how the kinetic mass changes.

Data for the spin-averaged kinetic masses from the tree-level ${\mathcal{O}(p^{4})}$ and ${\mathcal{O}(p^{6})}$ NRQCD actions, both with and without ${\mathcal{O}(\alpha_{s}p^{4})}$ corrections, are presented in Figures 11 and 12 for the ensembles described in Table 3. Errors are statistical only. Since these plots have the same vertical scale we can see a reduction in the size of $ap^{6}$ and $\mathcal{O}(\alpha_{s})$ effects as the lattice spacing is reduced, from the fact that the range of results become more compressed from Sets 1 to 5.

In Figure 13 we plot the differences between the energies of $\eta_{b}$ states with momentum $(2,2,1)$ and $(3,0,0)$, in units of $2\pi/L$, on each ensemble using different actions. We see that the largest $SO(3)$ breaking occurs for the ${\mathcal{O}(p^{4})}$ NRQCD action. This breaking is reduced when including the ${\mathcal{O}(\alpha_{s}p^{4})}$ kinetic couplings, and then reduced further by the ${\mathcal{O}(p^{6})}$ NRQCD action. The least $SO(3)$ breaking occurs for the ${\mathcal{O}(\alpha_{s}p^{4},p^{6})}$ NRQCD action. This improvement is sizable for the very coarse ensemble, Set 1, while for the coarse and fine ensembles the improvement is visible, but small. Further, the breaking on the coarse and fine ensembles goes from being a significant effect to a non-significant ($2\sigma$) one after improvement. Using this improved action allows for a more accurate and reliable determination of the kinetic mass, and hence also of the tuned $b$-quark mass in high precision calculations.

In Figure 14, we show the speed of light squared, c^{2}=\Big{(}(\Delta E+M_{\text{kin}})^{2}-M_{\text{kin}}^{2}\Big{)}/P^{2}, computed on Set 5 with the ${\mathcal{O}(\alpha_{s}p^{4},p^{6})}$ NRQCD action against $P^{2}a^{2}$, where good agreement with the value of 1 is seen.

V Discussion and Conclusions

In this work we have made the next round of improvement to the HPQCD collaboration’s formulation of the NRQCD action to allow increasingly accurate non-perturbative calculations in the future. The key results presented herein include:

• •

Determining the required operators which need to be added to the NRQCD action in order to give a correct heavy-quark dispersion relation to ${\mathcal{O}(p^{6})}$, presented in Sec. II.

• •

Determining the one-loop coefficients of the ${\mathcal{O}(p^{4})}$ kinetic couplings, namely $c^{(1)}_{1},c^{(1)}_{6}$ and $c_{5}^{(1)}$, in automated lattice perturbation theory using twisted boundary conditions as an IR regulator. We also present results for the one-loop (bare-to-pole) heavy-quark mass renormalisation $Z_{m}^{(1)}$ and zero-point energy $W_{0}$ which can be combined to give the one-loop energy shift (from neglecting the quark mass term in the NRQCD action) of a $b$-quark. This one-loop energy shift can be added to the non-perturbatively determined simulation energies to give a numerical value, which after converting to GeV, can be compared to the experimentally determined masses. All perturbative results are shown in Sec. III.2.

• •

Determining the full one-loop radiative correction of these quantities by finding the scale $q^{*}$ of $\alpha_{s}$ defined in the $V$-scheme. In doing this we use the higher order methodology which takes into account the anomalously small leading-order moments in order to obtain physical $q^{*}$ as described in Sec. II.2.

• •

Determining the one-loop quantities for three different NRQCD action formulations, namely a NRQCD action that gives a heavy-quark dispersion relation correct (i) to ${\mathcal{O}(p^{4})}$ and (ii) to ${\mathcal{O}(p^{6})}$. These actions employ a mean-field tadpole improvement procedure. For reasons described in Sec. II.2 we also explore, for the first time, a (iii) Fat$3$ smeared NRQCD action with the quark dispersion relation correct to ${\mathcal{O}(p^{4})}$ which does not require mean-field improvement. The Fat$3$ results are encouraging and show stable behaviour against $am_{b}$, indicating that the use of this or a similar smearing may be the way forward in future, rather than tadpole-improvement.

• •

Varying the stability parameter with $n=4,6$ and $8$ to show that, as shown in Sec. III.2, the ${\mathcal{O}(\alpha_{s}p^{4})}$ kinetic couplings in the NRQCD action are insensitive to this choice. Thus, if future calculations need to compensate a decrease in lattice spacing (which allows higher momentum fluctuations) with an increase in $n$, they can do so reliably.

• •

Testing how the improvement of the NRQCD action, both in terms of additional ${\mathcal{O}(p^{6})}$ operators and one-loop radiative kinetic coefficients, affects the non-perturbatively obtained kinetic masses, c. f. Figs. 11, 12 and 13. The impact of the $p^{6}$ terms and the radiative corrections on the kinetic masses obtained is small, particularly on the finer lattices. We find a significant reduction in $SO(3)$ symmetry breaking when using the improved actions on the very coarse ensemble, Set 1, which decreases as the lattice spacing is reduced. On the fine lattice, Set 5, $SO(3)$ symmetry breaking has been reduced to the point that the energy splitting, shown in Figure 13, is very nearly consistent with zero.

Taken together, NRQCD allows increasingly accurate and precise numerical calculations to be performed by including higher-order operators, in combination with determining the matching coefficients using perturbation theory. We have taken both these steps in this work. Furthermore, NRQCD is numerically cheap compared to its relativistic counterparts, being an initial-value, rather than a boundary-value problem. The outlook for NRQCD in the high-precision era is promising and this work helps ensure that this NRQCD formalism will continue to be an active contributor.

ACKNOWLEDGMENTS

We would like to thank Christopher Thomas for the many insightful discussions on finite-volume lattice effects. We are also grateful to the MILC collaboration for the use of their gauge configurations. This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U. S. Department of Energy, Office of Science, Office of High Energy Physics. This work was funded in part by UK Science and Technology Facilities Council (STFC) grants ST/L000385/1, ST/L000466/1 and ST/P000681/1. This work used the DiRAC Data Analytic system at the University of Cambridge, operated by the University of Cambridge High Performance Computing Service on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant (ST/K001590/1), STFC capital grants ST/H008861/1 and ST/H00887X/1, and STFC DiRAC Operations grant ST/K00333X/1. DiRAC is part of the National E-Infrastructure.

Appendix A Derivative Conventions

In this section we define our convention for the discretised derivative operators for use in the perturbative and non-perturbative calculations. The forward, backward and partially-cancelled second-order $\Delta^{(2)}_{j}$ operator are given by (note that all gauge-links are implicitly mean-field improved so that $U_{i}(x)$ is replaced by $U_{i}(x)/u_{0}$)

[TABLE]

Then our partially corrected operators are

[TABLE]

The additional terms in Eqs. (40), (41) and (42) are needed for the partial cancellation as described in Sec. II.2. As can be seen, when the operators are transformed to momentum space, the additional terms in these partially-cancelled operators allow mixing down of higher-order coefficients to lower-order tadpole counterterms. For the smeared operators, no mean-field improvement is performed (i.e. $u_{0}$ is set to 1) so the additional terms vanish and the links are replaced by their smeared counterparts.

Appendix B Tadpole Counterterms from Mean-field Improvement

In this appendix we will give explicit formulae for the tadpole counterterms used to remove the unphysical tadpole contributions, as described in Sec. II.2, when using a mean-field improvement procedure. These formulae are utilised to produce the one-loop mean-field improved quantities discussed in Sec. III.2. Features of these formulae have been discussed in Sec. II.2.

For a ${\mathcal{O}(p^{4})}$ NRQCD action (e.g., using Eq. (2) with $\delta H_{p^{6}}=0$) with partial cancellation, the tadpole counterterms are

[TABLE]

[TABLE]

For a ${\mathcal{O}(p^{6})}$ NRQCD action (e.g., using Eq. (2)) with partial cancellation, the tadpole counterterms are

[TABLE]

[TABLE]

Appendix C Numerical Results

In this section, we give numerical values for the tadpole-improved one-loop coefficients in Table 6. The tadpole counterterms given in Appendix B can be used with the mean-field improved data to produce the raw results. We also give $aq^{*}$ in Table 7 for each quantity which is used to determine the the value of $\alpha_{V}$. Lastly, a subset of the full one-loop radiative corrections relevant for heavy-quark non-perturbative calculations are given in Tables 8 and 9. Other values can be read off the Figures in Sec. III.2.