Determination of the QCD coupling from the static energy and the free energy

TL;DR

This paper provides two precise determinations of the strong coupling constant lpha_s at the Z boson mass scale using static energy and free energy lattice data, with detailed error analysis and updated results.

Contribution

It offers new lattice-based measurements of lpha_s with improved accuracy and shorter distance data, updating previous determinations and including comprehensive error estimates.

Findings

lpha_s(M_Z)=0.11660^{+0.00110}_{-0.00056} from static energy

lpha_s(M_Z)=0.11638^{+0.00095}_{-0.00087} from free energy data

Inclusion of new lattice data at smaller spacings and distances

Abstract

We present two determinations of the strong coupling \(\alpha_s\). The first one is from the static energy at three-loop accuracy, and may be considered an update of earlier determinations by some of us. The new analysis includes new lattice data at smaller lattice spacings, and reaches distances as short as \(0.0237\,{\rm fm}\). We present a comprehensive and detailed estimate of the error sources that contribute to the uncertainty of the final result, $\alpha_s(M_Z)=0.11660^{+0.00110}_{-0.00056}$. The second determination is based on lattice data for the singlet free energy at finite temperature up to distances as small as \(0.0081\,{\rm fm}\), from which we obtain $\alpha_s(M_Z)=0.11638^{+0.0009 5}_{-0.00087}$.

Tables

Table 1. Table 1: Parameters for the fine T = 0 𝑇 0 T=0 ensembles. In the seventh column we indicate the number of correlator measurements performed (MEAS).

$m_l=m_s/20$:
$\beta$	a (fm)	$N_{\sigma },N_{\tau }$	$a{m}_s$	$m_{\pi }L$	#TUs	#MEAS	Ref.
7.373	0.060	${64}^3\times 48$	0.0250	2.3	4623	1000	Bazavov et al. (2014b)
7.596	0.049	${64}^4$	0.0202	2.6	4757	1000	Bazavov et al. (2014b)
7.825	0.040	${64}^4$	0.0164	2.0	4768	1000	Bazavov et al. (2014b)
$m_l=m_s/5$:
$\beta$	a [fm]	$N_{\sigma },N_{\tau }$	$a{m}_s$	$m_{\pi }L$	#TUs	#MEAS	Ref.
8.000	0.035	${64}^4$	0.01299	3.6	4616	1000	Bazavov et al. (2018a)
8.200	0.029	${64}^4$	0.01071	3.1	4616	1000	Bazavov et al. (2018a)
8.400	0.025	${64}^4$	0.00887	2.6	4616	1000	Bazavov et al. (2018a)

Table 2. Table 6: Parameters for the coarse T = 0 𝑇 0 T=0 ensembles. In the seventh column we indicate the number of correlator measurements performed.

$m_l=m_s/20$:
$\beta$	a (fm)	$N_{\sigma },N_{\tau }$	$a{m}_s$	$m_{\pi }L$	#TUs	#MEAS	Ref.
6.664	0.117	${32}^4$	0.0514	3.0	4000	400	Bazavov et al. (2012b)
6.740	0.109	${48}^4$	0.0476	4.2	4000	800	Bazavov et al. (2012b)
6.800	0.103	${32}^4$	0.0448	2.7	4000	400	Bazavov et al. (2012b)
6.880	0.095	${48}^4$	0.0412	3.7	8200	820	Bazavov et al. (2012b)
6.950	0.089	${32}^4$	0.0386	2.3	19400	1940	Bazavov et al. (2012b)
7.030	0.082	${48}^4$	0.0356	3.2	7800	780	Bazavov et al. (2012b)
7.150	0.070	${64}^3\times 48$	0.0320	2.7	4047	667	Bazavov et al. (2012b)
7.280	0.065	${64}^3\times 48$	0.0284	2.5	3978	650	Bazavov et al. (2012b)
$m_l=m_s/5$:
$\beta$	a [fm]	$N_{\sigma },N_{\tau }$	$a{m}_s$	$m_{\pi }L$	#TUs	#MEAS	Ref.
7.030	0.082	${48}^4$	0.0712	6.4	2990	598	Bazavov et al. (2018a)
7.825	0.040	${64}^4$	0.01542	4.0	2980	298	Bazavov et al. (2018a)

Table 3. Table 7: Parameters of N τ = 12 subscript 𝑁 𝜏 12 N_{\tau}=12 ensembles with aspect ratio 4 4 4 and m l = m s / 20 subscript 𝑚 𝑙 subscript 𝑚 𝑠 20 m_{l}=m_{s}/20 . Adjacent correlators are separated by 10 10 10 TUs.

$m_l=m_s/20$:
$\beta$	$a{m}_s$	T (MeV)	#TUs	Ref.
7.373	0.0250	273	85120	Bazavov et al. (2014b)
7.596	0.0202	334	98010	Bazavov et al. (2014b)
7.650	0.0202	350	3230	Bazavov et al. (2013)
7.825	0.0164	408	134600	Bazavov et al. (2014b)
8.570	0.008376	770	6320	Bazavov et al. (2018b)
8.710	0.007394	866	6490	Bazavov et al. (2018b)
8.850	0.006528	974	6340	Bazavov et al. (2018b)
9.060	0.004834	1162	7430	Bazavov et al. (2018b)
9.230	0.004148	1340	7280	Bazavov et al. (2018b)
9.360	0.003691	1495	7910	Bazavov et al. (2018b)
9.490	0.003285	1667	9780	Bazavov et al. (2018b)
9.670	0.002798	1938	7650	Bazavov et al. (2018b)
$m_l=m_s/5$:
$\beta$	$a{m}_s$	T (MeV)	#TUs	Ref.
8.000	0.01299	474	71670	Bazavov et al. (2018b)
8.200	0.01071	562	71390	Bazavov et al. (2018b)
8.400	0.00887	667	71170	Bazavov et al. (2018b)

Table 4. Table 8: Parameters of the new N τ = 16 subscript 𝑁 𝜏 16 N_{\tau}=16 ensembles with aspect ratio 4 4 4 and expectation values of bare Polyakov loops. Adjacent correlators are separated by 10 10 10 TUs.

$m_l=m_s/20$:
$\beta$	$a{m}_s$	T (MeV)	#TUs	$L^{\mathrm{bare}}$, Ref.
7.825	0.00164	306	67960	0.002962(06)
8.570	0.008376	577	10400	Bazavov et al. (2018b)
8.710	0.007394	650	10190	Bazavov et al. (2018b)
8.850	0.006528	731	4480	Bazavov et al. (2018b)
9.060	0.004834	872	41870	0.015993(20)
9.230	0.004148	1005	3610	Bazavov et al. (2018b)
9.360	0.003691	1121	3530	Bazavov et al. (2018b)
9.490	0.003285	1250	6790	Bazavov et al. (2018b)
9.670	0.002798	1454	42060	0.025530(29)
$m_l=m_s/5$:
$\beta$	$a{m}_s$	T (MeV)	#TUs	$L^{\mathrm{bare}}$
8.000	0.001299	356	11460	Bazavov et al. (2018b)
8.200	0.001071	422	10660	Bazavov et al. (2018b)
8.400	0.000887	500	64370	0.007819(10)

Table 5. Table 9: For many data sets we have to restrict the fit windows in the ground state extraction due to instabilities at large τ / a 𝜏 𝑎 \tau/a . Otherwise, the fit window reaches up to τ / a = 12 𝜏 𝑎 12 \tau/a=12 .

$\beta =7.373$:
$\text{mathbold}r/a$	$\max \left(\tau /a\right)$
(1,0,0)	10
(2,1,0)	10
(2,2,2)	10
(3,3,1)	10
(4,0,0)	10
(4,1,0)	10
(4,2,0)	10

Table 6. Table 10: Systematic uncertainty estimates for the discretization artifacts from the static energy at T = 0 𝑇 0 T=0 or the singlet free energy with N τ = 16 subscript 𝑁 𝜏 16 N_{\tau}=16 or 12 12 12 at r / a = 3 𝑟 𝑎 3 r/a=3 . The estimates have been calculated individually for each data set. The fluctuations within a factor 2 or 3 between the uncertainty estimates for the different data sets are due to statistical fluctuations of the data.

$\beta$	$N_{\tau }$	$r_1{\delta }_{\mathrm{abs}}^{\mathrm{syst}}$	${\delta }_{\mathrm{rel}}^{\mathrm{syst}}\left(\frac{r}{a}=3\right)$
7.596	64	0.0041	0.00103
7.825	64	0.0041	0.00093
8.000	64	0.0038	0.00076
8.200	64	0.0037	0.00066
8.400	64	0.0034	0.00054
7.825	16	0.0161	0.00360
8.000	16	0.0179	0.00362
8.200	16	0.0161	0.00290
8.400	16	0.0108	0.00172
8.570	16	0.0055	0.00079
8.710	16	0.0057	0.00074
8.850	16	0.0071	0.00086
9.060	16	0.0078	0.00082
9.230	16	0.0084	0.00079
9.360	16	0.0098	0.00085
9.490	16	0.0111	0.00089
9.670	16	0.0054	0.00038
7.596	12	0.0052	0.00128
7.650	12	0.0065	0.00156
7.825	12	0.0052	0.00116
8.000	12	0.0038	0.00078
8.200	12	0.0033	0.00060
8.400	12	0.0028	0.00045
8.570	12	0.0030	0.00044
8.710	12	0.0041	0.00054
8.850	12	0.0079	0.00095
9.060	12	0.0097	0.00102
9.230	12	0.0114	0.00108
9.360	12	0.0066	0.00057
9.490	12	0.0058	0.00046
9.670	12	0.0049	0.00035

Table 7. Table 11: Tadpole factors u 0 subscript 𝑢 0 u_{0} for the lattices used in the improvement procedure. Tadpole factors have been obtained using zero temperature or finite temperature lattices. Temperature effects are always at most 0.00001 0.00001 0.00001 .

$\beta$	$u_0$	$\beta$	$u_0$	$\beta$	$u_0$
6.608	0.86823	7.030	0.88173	8.200	0.90604
6.664	0.87026	7.100	0.88363	8.400	0.90909
6.700	0.87152	7.150	0.88493	8.570	0.91152
6.740	0.87288	7.200	0.88621	8.710	0.91341
6.770	0.87388	7.280	0.88817	8.850	0.91522
6.800	0.87485	7.373	0.89035	9.060	0.91778
6.840	0.87613	7.596	0.89517	9.230	0.91973
6.880	0.87736	7.650	0.89627	9.360	0.92116
6.910	0.87827	7.825	0.89962	9.490	0.92253
6.950	0.87945	8.000	0.90274	9.670	0.92435

Equations

$$r^2{\frac{dE(r)}{dr}|}_{r=r_1}=1,$$

$$F(r)$$

$$\displaystyle+\frac{\alpha_{s}(1/r)}{4\pi}\Big{(}\tilde{a}_{1}-2\beta_{0}\Big{)}$$

$$\displaystyle+\frac{\alpha_{s}^{2}(1/r)}{(4\pi)^{2}}\Big{(}\tilde{a}_{2}-4\tilde{a}_{1}\beta_{0}-2\beta_{1}\Big{)}$$

$$\displaystyle+\frac{\alpha_{s}^{3}(1/r)}{(4\pi)^{3}}\Big{(}\tilde{a}_{3}-6\tilde{a}_{2}\beta_{0}-4\tilde{a}_{1}\beta_{1}-2\beta_{2}$$

$$\displaystyle+a_{3}^{L}\ln\frac{C_{A}\alpha_{s}(1/r)}{2}\Big{)}+\mathcal{O}(\alpha_{s}^{4},\alpha_{s}^{4}\ln^{2}\alpha_{s})\Bigg{]},$$

$$F(r,\nu =1/r)$$

$$\displaystyle+\frac{\alpha_{s}(1/r)}{4\pi}\Big{(}\tilde{a}_{1}-2\beta_{0}\Big{)}$$

$$\displaystyle+\frac{\alpha_{s}^{2}(1/r)}{(4\pi)^{2}}\Big{(}\tilde{a}_{2}-4\tilde{a}_{1}\beta_{0}-2\beta_{1}\Big{)}$$

$$-\frac{{\alpha }_s^2(1/r)}{(4\pi {)}^2}\frac{a_3^L}{2{\beta }_0}\ln \frac{{\alpha }_s({\mu }_{us})}{{\alpha }_s(1/r)}$$

$$+\frac{{\alpha }_s^2(1/r){\alpha }_s({\mu }_{us})}{(4\pi {)}^3}{a}_3^L\ln \frac{C_A{\alpha }_s(1/r)}{2r{\mu }_{us}}$$

$$\displaystyle+\frac{\alpha_{s}^{3}(1/r)}{(4\pi)^{3}}\Big{(}\tilde{a}_{3}-6\tilde{a}_{2}\beta_{0}-4\tilde{a}_{1}\beta_{1}-2\beta_{2}\Big{)}$$

$$\displaystyle+\mathcal{O}(\alpha_{s}^{4})\Bigg{]}.$$

$${\alpha }_s(M_Z,N_f=5)=0.11660_{-0.00056}^{+0.00110},$$

$$\delta {\alpha }_s(M_Z,N_f=5)=(41{)}^{\mathrm{stat}}(21{)}^{\mathrm{lat}}(10{)}^{r_1}{(}_{-13}^{+95}{)}^{\mathrm{soft}}(28{)}^{\mathrm{us}},$$

$${\Lambda }_{\overline{\mathrm{MS}}}^{N_f=3}=314.0_{-8.0}^{+15.5}\mathrm{MeV},$$

$$\delta {\Lambda }_{\overline{\mathrm{MS}}}^{N_f=3}=(5.8{)}^{\mathrm{stat}}(3.0{)}^{\mathrm{lat}}(1.7{)}^{r_1}{(}_{-1.8}^{+13.4}{)}^{\mathrm{soft}}(4.0{)}^{\mathrm{us}}\mathrm{MeV}.$$

$${\delta }^{\sqrt{2}}{\alpha }_s(M_Z,N_f=5)=(41{)}^{\mathrm{stat}}(21{)}^{\mathrm{lat}}(10{)}^{r_1}{(}_{-13}^{+37}{)}^{\mathrm{pert}},$$

$${\delta }^{\sqrt{2}}{\Lambda }_{\overline{\mathrm{MS}}}^{N_f=3}=(5.8{)}^{\mathrm{stat}}(3.0{)}^{\mathrm{lat}}(1.7{)}^{r_1}{(}_{-1.8}^{+5.2}{)}^{\mathrm{pert}}\mathrm{MeV}.$$

$$F_S(r,T)=-T\ln \left(\frac{1}{N_c}⟨\mathrm{Tr}\left[W(r)W^{†}(0)\right]⟩\right).$$

$$F_S(r,T)=V_s(r,{\mu }_{us})+\delta F_S(r,T,{\mu }_{us}),$$

$${\alpha }_s(M_Z,N_f=5)=0.11638_{-0.00087}^{+0.00095},$$

$$\delta {\alpha }_s(M_Z,N_f=5)=$$

$$(80{)}^{\mathrm{stat}}(21{)}^{\mathrm{lat}}(17{)}^{T>0}(10{)}^{r_1}{(}_{-06}^{+40}{)}^{\mathrm{soft}}(15{)}^{\mathrm{us}},$$

$${\Lambda }_{\overline{\mathrm{MS}}}^{N_f=3}=310.9_{-12.3}^{+13.5}\mathrm{MeV},$$

$$\delta {\Lambda }_{\overline{\mathrm{MS}}}^{N_f=3}=$$

$$(11.3{)}^{\mathrm{stat}}(3.0{)}^{\mathrm{lat}}(2.4{)}^{T>0}(1.7{)}^{r_1}{(}_{-0.8}^{+5.6}{)}^{\mathrm{soft}}(2.1{)}^{\mathrm{us}}\mathrm{MeV}.$$

$$E_{\mathrm{tree}}^{\mathrm{lat}}(r)=-\frac{4}{3}{g}_0^2\int \frac{d^{3}k}{(2\pi {)}^3}{e}^{i\mathbf{kr}}{D}_{00}(k),$$

$$E_{\mathrm{tree}}^{\mathrm{lat}}(r)=-\frac{4}{3}{g}_0^2\frac{1}{4\pi r_I}$$

$$E_{\mathrm{imp}}(r_I)=E^{\mathrm{lat}}(r),$$

$$E_{\mathrm{imp}}(r)=E^{\mathrm{lat}}(r)+{\kappa }_I\left(\frac{\text{mathbold}r}{a},\beta \right)\frac{r_I-r}{r{r}_I}.$$

$$K\left(\frac{\text{mathbold}r}{a},\beta \right)=\frac{E^{\mathrm{lat}}\left(\text{mathbold}r,a(\beta )\right)}{E^{\mathrm{cont}}(r_I)}-1,$$

$$E_{\mathrm{imp}}^{\mathrm{lat}}(r,\beta )=\frac{E^{\mathrm{lat}}\left(\frac{\text{mathbold}r}{a},a(\beta )\right)}{1+K\left(\frac{\text{mathbold}r}{a},\beta \right)},$$

$$\displaystyle\Braket{K\left(\frac{\mathbold r}{a},\alpha_{s}^{\rm lat}\right)}=\sum\limits_{i=1}^{N}\kappa_{i}\left(\frac{\mathbold r}{a}\right)\left(\alpha_{s}^{\rm lat}\right)^{i}.$$

$$C_F=\frac{N_c^2-1}{2N_c},$$

$$\frac{d{\alpha }_s(\nu )}{d\ln \nu }$$

$$=-2{\alpha }_s\left[{\beta }_0\frac{{\alpha }_s}{4\pi }+{\beta }_1{\left(\frac{{\alpha }_s}{4\pi }\right)}^2+\cdots \right],$$

$${\beta }_0$$

$${\beta }_1$$

$${\beta }_2$$

Full text

††preprint: TUM-EFT 111/18; INT-PUB-19-028;

I Introduction

A precise determination of the strong coupling $\alpha_{s}$ is of key importance both for the theory of strong interactions, which is determined by such fundamental parameter, and for investigations of new physics beyond the Standard Model. Cross section calculations at the Large Hadron Collider, for example suffer from the uncertainties related to our limited knowledge of $\alpha_{s}$. The last decades have witnessed an impressive effort in the $\alpha_{s}$ extraction from an ample range of physical observables with a broad sweep of different methods. Notwithstanding all these efforts, the Particle Data Book (PDG) world average value for $\alpha_{s}$ in 2018, $\alpha_{s}(M_{Z})=0.1181\pm 0.0011$, has an overall uncertainty that has almost doubled with respect to the PDG average in 2013, $\alpha_{s}(M_{Z})=0.1185\pm 0.0006$ Tanabashi et al. (2018). This is due to the fact that the uncertainty of several determinations that enter the average is dominated by errors of theoretical origin, which are often difficult to precisely assess. It is important therefore to use a variety of different ways to extract $\alpha_{s}$ and to validate each extraction at the best of the state of the art.

In this paper we aim at an improved determination of $\alpha_{s}$ both by updating earlier extractions from the static energy and by proposing a new method that involves the singlet free energy.

The static energy is an observable up to an additive constant, and it is a function of the distance $r$ between the static quark and the static antiquark. In the limit of massless dynamical quarks, it depends only on $\alpha_{s}$. It can be calculated on the lattice for any distance $r=\sqrt{x^{2}+y^{2}+z^{2}}$, where $x$, $y$ and $z$ are integer multiples of the lattice spacing. For short distance $r$ it can be calculated in perturbation theory in QCD as a function of $\alpha_{s}$ using nonrelativistic QCD effective field theories (EFT). Noteworthy, perturbation theory is accurate for this quantity at three loops, and the tree-level result is already sensitive to $\alpha_{s}$. The comparison between the perturbative expression and lattice data at distances small enough to be accessible to perturbation theory is a good way to obtain a precise determination of $\alpha_{s}$. For this endeavor it is crucial to have lattice data covering an interval of sufficiently small distances in order to be sensitive to the minute details of the curvature of the static energy, namely, with sufficiently fine lattice spacing. At the present day, these lattices still pose a major challenge. In this paper we will use lattices Bazavov et al. (2018a) with an extraordinarily fine lattice spacing $a=0.0249\,{\rm fm}$ to achieve a systematically improved extraction of $\alpha_{s}$.

Additionally we exploit a new idea. One reason for which it is challenging to reach such fine lattice spacings in lattice QCD simulations with dynamical quarks is that one has to simultaneously maintain the control over finite volume effects arising from the propagation of the lightest hadronic modes, namely, the Goldstone bosons, at the pion scale. A lattice simulation at high enough temperature avoids this infrared problem, and thus enables reaching much finer lattice spacings using smaller volumes. We use finite temperature lattices with unprecedentedly fine lattice spacing $a=0.00848\,{\rm fm}$ Bazavov et al. (2018b). The singlet static free energy is again a function of the static quark-antiquark distance and has been calculated on the lattice Bazavov et al. (2018b) and perturbatively Berwein et al. (2017). The comparison between the two offers a novel and independent method to extract a precise determination of $\alpha_{s}$.

The rest of the paper is organized as follows. In Section II we give an account of the gauge ensembles and lattice correlators that we use and discuss the relevant systematic effects in the lattice data. In Section III after briefly recalling the procedure used in Ref. Bazavov et al. (2014a) to extract $\alpha_{s}$ from the static energy, we compare between the lattice data and the weak-coupling result, discuss our estimates of systematic uncertainties, and finally obtain an updated extraction of $\alpha_{s}$. In Section IV we discuss the singlet free energy at finite temperature, and outline the relation with the static potential and static energy at zero temperature. We discuss the relevant scale hierarchies, compare with the lattice, and proceed to extract $\alpha_{s}$ at even shorter distances. Section V contains some discussion of our results, as well as a comparison with previous related works. Finally, in Sec. VI we present a short summary of the main results and conclude. In Appendix A we discuss the gauge ensembles used in this study. In Appendix B we discuss discretization artifacts at short distances in detail. In Appendix C we list the coefficients appearing in the perturbative results used in this paper.

II Lattice setup

We require the lattice result of the static energy at the smallest available distances in order to compare with the weak-coupling calculations. For this reason we employ the 2+1 flavor ensembles using the highly improved staggered quark (HISQ) action Follana et al. (2007) that have been generated for studies of the QCD Equation of State Bazavov et al. (2014b, 2018a). The spatial volume is given by $V_{3}=(aN_{\sigma})^{3}$, and the physical length of the Euclidean time direction is $aN_{\tau}$. We summarize the zero temperature ensembles111 The lattice observables (gluon action density, and sea quark condensates) contributing to the Equation of State have divergent, additive contributions, which depend only on the lattice spacing. Renormalization group invariant observables are obtained by subtracting the $T=0$ observables from the $T>0$ observables at each lattice spacing.

with $N_{\tau}=64$ in Tab. 1, which correspond to two different pion masses in the continuum limit. Namely, these ensembles have different light sea quark masses, $m_{l}=m_{s}/20$ or $m_{l}=m_{s}/5$. The strange quark mass $m_{s}$ is at its physical value, namely, such that the mass of the hypothetical $\eta_{s\bar{s}}$ meson is reproduced as $m_{\eta_{s\bar{s}}}=686\,{\rm MeV}$. The lattice spacing $a$ has been fixed by the $r_{1}$ scale, which is defined by the equation

[TABLE]

see Refs. Bazavov et al. (2014b, 2018a) for the details. Here, and in the following $E(r)$ denotes the QCD static energy calculated on the lattice. We use the value $r_{1}=0.3106(17)\,{\rm fm}$ determined from the pion decay constant. Since the quark mass dependence of $r_{1}/a$ is small, the $m_{l}=m_{s}/20$ and $m_{l}=m_{s}/5$ results can be combined to produce a parametrization of $r_{1}/a$ as a function of the bare lattice gauge coupling $\beta$ Bazavov et al. (2018a). The tunneling between different topological sectors is suppressed for these gauge ensembles. Nevertheless, no statistically significant difference has been observed for bulk observables other than the light quark condensate Bazavov et al. (2018a) in different topological sectors. For the distances we consider in this study the static energy Weber et al. (2018) does not vary significantly between the different topological sectors.

In addition, we use the finite temperature results of the singlet free energy with $N_{\sigma}/N_{\tau}=4$ and $N_{\tau}=12$, or $16$. These ensembles correspond to the thermal QCD medium at temperatures $T=1/(aN_{\tau})$. The finite temperature ensembles have been generated using lattice parameters that would correspond to the same two pion masses (at zero temperature) in the continuum limit as the zero temperature ensembles Bazavov et al. (2018b). For these we use the same parametrization of $r_{1}/a$. We give an account of these ensembles in Tables 7 and 8 in Appendix A.

The static energy can be calculated using Wilson line correlators in Coulomb gauge or Wilson loops. We use the former because it is more convenient for practical reasons. For additional crosschecks we also studied Wilson loops for $\beta=7.825$. When extracting the static energy from Wilson loops the spatial lines have to be smeared in order to obtain a reasonable signal-to-noise ratio and suppress excited states. We used one, two and five levels of HYP smearings Hasenfratz and Knechtli (2001) in our calculations. In the studies of the Equation of State Bazavov et al. (2014b, 2018a), the static energy was extracted from two-exponential fits to the Wilson line correlators using a fixed range of Euclidean time $\tau/a$ that was chosen independently from $r/a$. This previous calculation is not adequate for our analysis. Here we improve the extraction of the static energy in the following way.

• •

We explicitly include the correlation matrix of the static energy in the fits to determine $\alpha_{s}$. For this reason we have to recalculate the static energy and preserve its statistical correlations.

• •

We explicitly distinguish at short distances between data with distinct path geometries instead of averaging over them. This increases the number of data for each $\beta$ compared to the earlier analysis Bazavov et al. (2014a).

• •

At small $\tau$ the correlators of Wilson lines are contaminated by excited states, while at large $\tau$ the statistical fluctuations are large, and the effective masses sometimes do not show a clear plateau, see Fig. 1 for some of the worst-case scenarios. Therefore a careful choice of the lower end, $\tau_{\rm min}$, and upper end, $\tau_{\rm max}$, of the fit interval is needed to obtain reliable results of the static energy. The corresponding fit interval $[\tau_{\rm min}:\tau_{\rm max}]$ will depend on the distance $r$. We estimated the systematic errors due to the choice of the fit interval. The procedure of extracting the static energy is demonstrated in Fig. 1, where the solid lines show the static energy and the dashed lines the corresponding total uncertainty. We summarize the fit windows used for obtaining the ground state in Tab. 9 in Appendix A.

At short distances, namely, for $r\lesssim 0.14\,{\rm fm}$, we observe the non-monotonic behavior of the effective mass at small $\tau/a$, typically, within $\tau\lesssim 0.14\,{\rm fm}$. This non-monotonic behavior is not restricted to the Wilson line correlator in Coulomb gauge, but also appears in Wilson loops with spatially-smeared spatial Wilson lines, see Fig. 2. For this reason, these non-monotonicities are due to the gauge ensembles with improved gauge action instead of being due to the details of the interpolating operator. The effective masses obtained from correlators of Wilson lines in Coulomb gauge and from Wilson loops with HYP smeared spatial lines converge to the same plateau at large $\tau$.

The current and previous results for the static energy are consistent within errors. For the given reasons the new result is more precise at short distances and has a better estimate of the errors. It supersedes the previous calculation Bazavov et al. (2014b, 2018a).

The last major source of systematic uncertainties of the lattice data is due to discretization artifacts at short distances, which are caused by the breaking of rotational symmetry from O(3) in the continuum to the cubic group $W_{3}$ on the lattice. Even for an improved gauge action these artifacts are significantly larger than any of the other uncertainties of the lattice result. The artifacts express themselves in the form of a variation of the lattice static energy around the continuum static energy at any given distance $r$. The sign and the size of this variation depends on the geometry of the shortest paths connecting the two sites at distance $r$ on the lattice. This variation is typically large for on-axis distances, i.e., $\mathbold r/a=(n,0,0)$, while being smaller for off-axis distances, and shrinks to the point of numerical irrelevance as the distances become larger than $r/a\gtrsim 5$. These lattice artifacts could be understood qualitatively in leading-order perturbation theory, i.e. at tree level. The analysis presented in Appendix B shows that the large artifacts at $r/a=1$ are about $8\%$ for the Wilson action and $4\%$ for the Lüscher–Weisz action that we use, see Fig. 12. The tree-level discretization effects decrease rapidly at larger distances, and become very small for $r/a\geq\sqrt{8}$. One can use the tree-level results to reduce lattice artifacts also in the interacting theory. This procedure is called the tree-level improvement and discussed in Appendix B.

The statistical errors of the static energy obtained on the lattice increase when increasing the separation $r$. For $r/a\geq\sqrt{8}$ the statistical errors in the static energy are typically larger than the discretiation errors after tree-level improvement. Thus for these distances discretization errors can be neglected. For $\sqrt{5}\leq r/a\leq\sqrt{8}$ the statistical and discretization errors are comparable. In this region one can deal with the discretization errors by estimating them and combining with the statistical errors. For very short distances $r/a<\sqrt{5}$ the discretization effects are much larger that the statistical errors and should be estimated carefully and corrected for, before the lattice results can be compared to the weak-coupling calculations. The corresponding procedure is described in Appendix B.

III Extracting $\alpha_{s}$ from the static energy

The static energy of a quark-antiquark pair is directly accessible in lattice calculations. Our lattice calculation of the static energy for $N_{f}=2+1$ light flavors has been discussed in Sec. II. In the weak-coupling regime the static energy is known up to three loops, i.e., up to order $\alpha_{s}^{4}$ in the small $\alpha_{s}$ expansion. The resummation of $\ln\alpha_{s}$ terms is also known at subleading order. Hence, since realistic lattice calculations for $N_{f}=2+1$ are available at short enough distances where the small $\alpha_{s}$ expansion is reliable, $\alpha_{s}$ can be extracted by comparing the lattice data for the static energy with its perturbative expression for the case of three (approximately) massless quarks.

This program was initially carried out in Ref. Bazavov et al. (2012a) and considerably improved in Ref. Bazavov et al. (2014a), to a large extent due to the incorporation of shorter distance data on finer lattices. Here, we further improve the $\alpha_{s}$ extraction by including three even finer lattices. The abundance of data at short distances allows us to perform the analysis omitting the data with $r/a<\sqrt{5}$ that are susceptible to significant discretization artifacts. It also permits us to quantify the impact of using such data with or without appropriate corrections. Rather than making a separate analysis for each lattice spacing and averaging the final results, here we make a global fit to all lattice data leaving as a free parameter a normalization constant for each lattice spacing. This approach is well-motivated from the earlier observation that no remaining lattice spacing dependence can be resolved from the three lattice spacings in Ref. Bazavov et al. (2014a) and that the scale uncertainties for all considered lattice spacings are similar. We also take into account the statistical correlations between the fluctuations of data at different distances on the lattice. Eventually the effect of these correlations turns out to be much smaller than the statistical uncertainty, even though it reduces the uncertainty of the extracted value of $\alpha_{s}$.

We focus on data for $r\lesssim 0.14\,{\rm fm}$ or $r\lesssim 0.45\,r_{1}$, see Fig. 3. The shortest distance that we can access, namely, due to a single lattice spacing222 Actually, this is 4% smaller than the lattice spacing itself due to the tree-level improvement.

on our finest lattice (at zero temperature) is $0.0237\,{\rm fm}$. In Ref. Bazavov et al. (2014a) only data with $0.057\,{\rm fm}\lesssim r\lesssim 0.16\,{\rm fm}$ was used for the extraction. For that range it was shown that the small $\alpha_{s}$ expansion was reliable. Since we use only shorter distances in this paper, the perturbative expansion remains reliable. In Ref. Bazavov et al. (2014a), separate analyses with data for each gauge ensemble were performed in the range $r/a\geq\sqrt{2}$. In this work, thanks to the finer lattices and the combined analysis, we can afford to vary the minimal distance $r/a$ and quantify the corresponding uncertainty. Data on coarser lattices in the same $r$ range as well as data on finer lattices may be excluded due to restriction of the $r/a$ range in the analysis. This is shown with more detail in Fig. 3. It was also shown in Ref. Bazavov et al. (2014a) that there is no sensitivity to nonperturbative contributions giving rise to power-like corrections. This is important, since sea quark mass effects, finite volume effects, or effects due to an incorrect expectation value of the topological charge would make themselves noticeable in such contributions333 These contributions are also related to the second renormalon of the singlet potential, which needed to be accounted for explicitly in Refs. Takaura et al. (2019a, b) because the range of data used in those references extends to much larger distances. , if they were not negligibly small.

To obtain the perturbative result for the static energy we follow Ref. Bazavov et al. (2014a). We use the perturbative result for the force, $F(r)=dE/dr$, and integrate it to obtain the energy. The constant contribution obtained this way is irrelevant and subsumed into the normalization constants that we have to fix anyway when comparing the weak-coupling results to the lattice data at each $\beta$. A choice of the renormalization scale $\nu$ (also called the soft scale, see below) proportional to $1/r$ avoids large logarithms of the form $\ln(\nu r)$ Bazavov et al. (2014a). We call the choice $\nu=1/r$ our standard choice for the renormalization scale. For this choice we have Bazavov et al. (2014a)

[TABLE]

where the first three terms correspond to tree-level, one-loop, and two-loop order, respectively, the fourth term corresponds to three-loop order. The strong coupling, $\alpha_{s}$, is understood, here and in the following, as renormalized in the $\overline{\rm MS}$ scheme. At any time the strong coupling can be traded for the corresponding QCD scale. At three flavors in the $\overline{\rm MS}$ scheme the QCD scale is $\Lambda_{\rm\overline{MS}}^{N_{f}=3}$. The coefficients in Eq. (III) are defined in Ref. Garcia i Tormo (2013) and reproduced in Appendix C. At three loops there is a contribution proportional to $\ln\alpha_{s}$ as well as a non-logarithmic term. The origin of these terms can be understood using the effective field theory approach, namely, potential non-relativistic QCD (pNRQCD) Brambilla et al. (2000, 1999). The term, pointed out in Ref. Appelquist et al. (1977), proportional to $\ln\alpha_{s}$ comes from the energy scale $\sim\alpha_{s}/r$, which is called the ultra-soft scale to distinguish it from the energy scale $1/r$, which is referred to as the soft scale. Equation (III) is accurate at next-to-next-to-next-to-leading order (N3LO).

For sufficiently small values of the coupling, i.e., at small enough distances the ultra-soft logarithm $\ln(C_{A}\alpha_{s}(1/r)/2)$ eventually becomes large and can be resummed. This can be done using the renormalization group equations in pNRQCD Pineda and Soto (2000); Brambilla et al. (2005, 2009). The expression for the force after performing the resummation of the ultra-soft logarithms reads Bazavov et al. (2014a)

[TABLE]

We choose as a standard value for the ultra-soft scale $\mu_{us}=C_{A}\alpha_{s}(1/r)/(2r)$ where $C_{A}=N_{c}$. The coefficients in Eq. (III) are defined in Ref. Garcia i Tormo (2013) and reproduced in Appendix C.

At leading order in $\ln\alpha_{s}$ we can set $\mu_{us}$ to $1/r$ and Eq. (III) reduces to Eq. (III). The first four terms of Eq. (III) include all terms of the form $\alpha_{s}^{3+n}\ln^{n}\alpha_{s}$, i.e. it is accurate at next-to-next-to-leading logarithmic order (N2LL). Keeping all terms in Eq. (III), we obtain an expression that is accurate at three-loop order with leading ultra-soft resummation.

In Fig. 4, we show the two-loop contribution to the force, i.e., the third term of Eq. (III), along with different $\alpha_{s}^{4}$ contributions as function of the distance $r$. Here and in the rest of the paper we use, if not differently specified, the standard scales $\nu=1/r$ and $\mu_{us}=C_{A}\alpha_{s}(1/r)/(2r)$. One can see that the $\alpha_{s}^{4}\ln\alpha_{s}$ contribution is never larger than the non-logarithmic $\alpha_{s}^{4}$ contribution in the entire distance range used in our study, and both contributions are smaller than the two-loop contribution. If we resum the ultra-soft logarithms (i.e. by replacing the last line of Eq. (III) by the fourth and fifth lines of Eq. (III)), the corresponding term becomes slightly larger in absolute value, see the dot-dashed line in Fig. 4, but it is still significantly smaller than the two-loop contribution. Therefore, as already noticed in Bazavov et al. (2014a), the leading ultrasoft logarithms can be counted as being of the same order as the other non-logarithmic three-loop terms, while subleading logarithms may be included in the uncertainty due to the missing four-loop contributions. In Fig. 4, we see also that there is a partial cancellation between the soft and ultra-soft contributions at order $\alpha_{s}^{4}$. This cancellation, however, may be accidental and may not happen at higher orders. Our final result will come from comparing lattice data with perturbation theory at three-loop order with leading ultra-soft resummation, and we will use the difference between the N3LO expression of Eq. (III) and the three-loop expression with leading ultra-soft resummation from Eq. (III) to estimate the size of unknown higher-loop contributions. Other ways to estimate the higher-loop contributions will be discussed in the following and will enter eventually in our final error assessment.

We proceed as follows. On the lattice data side, we put all six data sets (see Tab. 1) together by adding an arbitrary subtraction constant to each set that is determined by the fit. On the continuum side, we generate a grid of $\alpha_{s}(M_{Z},N_{f}=5)$ input values ranging from $0.1140$ to $0.1200$ in steps of $0.0001$. We use the perturbative result for the force444 We use perturbative running and decoupling to convert each grid value of $\alpha_{s}(M_{Z},N_{f}=5)$ to the corresponding $\Lambda_{\rm\overline{MS}}^{N_{f}=3}$ value, which we use in the analytical expressions of the force, Eqs. (III) and (III). The details of the running and decoupling are the same as in Ref. Bazavov et al. (2014a). Namely, we use four-loop running, with the charm quark mass equal to $1.6\,{\rm GeV}$ and the bottom quark mass equal to $4.7\,{\rm GeV}$. The effects of higher-order terms in the running are negligible with the current accuracies. , given in Eq. (III), restricting ourselves to the three-loop with leading ultra-soft resummation expression, for the reasons discussed above. For each of the $\alpha_{s}$ input values we calculate the static energy $E(r)$ by numerically integrating Eq. (III) and fit the lattice data to it, thus determining the subtraction constants. Hence, the subtraction constants for different sets are obtained in independent fits, since the different ensembles are statistically independent. Of course the subtraction constants depend on the input value of $\alpha_{s}$. For the determination of the subtraction constants we consider fits with or without the full correlation matrix, which does not lead to statistically distinguishable results. The smallest $\chi^{2}/{\rm d.o.f.}$ determines the preferred value of $\alpha_{s}(M_{Z})$ by the lattice data. We calculate the $\chi^{2}/{\rm d.o.f.}$ with or without the full correlation matrix. Typically, the minimum of the $\chi^{2}/{\rm d.o.f.}$ varies by at most $\delta^{\rm corr}=\pm 0.0001$ upon inclusion of the full correlation matrix, which is always smaller than the statistical error. We consider this difference as part of the statistical error that is already accounted for. The subtraction constants in the correlated fit are also consistent with the subtraction constants in the uncorrelated fits within fractions of the statistical error. We repeat this procedure for several values of $\max(r)$, vary the smallest distances $\min(r/a)$ in units of the lattice spacing considered in the analysis, and vary the treatment of discretization artifacts at the shortest distances. Eventually we repeat the entire procedure for expressions with different values of the soft scale, and for the two-loop and three-loop (without resumming the ultra-soft contribution) expressions.

Discretization artifacts are significant at distances $r/a\lesssim 3$, with increasing importance at smaller and smaller distances as detailed in Sec. II. We use different strategies to handle the residual discretization artifacts, see Appendix B for details:

1. (i)

ignoring them, i.e., we engage in no further corrections beyond the tree-level improvement. In this case we restrict the fit to $r/a\geq\sqrt{8}$ and omit $r/a=\sqrt{12}$; 2. (ii)

correcting for them, by following a strategy similar to Ref. Bazavov et al. (2014a), i.e., nonperturbative improvement (see Appendix B); 3. (iii)

accounting for them, by enlarging the statistical errors to fully cover the discretization artifacts, i.e., tree-level improvement and systematic error estimates. This approach is suitable for $r/a\geq\sqrt{5}$.

First, we analyze the data for $r/a\geq\sqrt{8}$ without considering the discretization artifacts beyond the tree-level improvement. In these fits we omit the data with $r/a=\sqrt{12}$ due to their exceptionally large discretization artifacts. We generally obtain $\chi^{2}/{\rm d.o.f.}$ between 0.3 and 0.6 for the best fits, which rises steeply outside a narrow window around the central value of $\alpha_{s}$. We estimate the systematic error due to discretization artifacts from the difference to fits with $r/a\geq\sqrt{8}$ with tree-level improvement or with nonperturbative improvement as being typically about $\delta^{\rm lat}=\pm 0.00021$. Tab. 2 clearly shows that any errors due to discretization artifacts are safely contained in the statistical uncertainty. We point out that the analysis for $\max(r)\leq 0.12\,{\rm fm}$ and $r/a\geq\sqrt{8}$ does not depend on gauge ensembles with the smaller sea quark mass, i.e., uses exclusively ensembles with $m_{l}=m_{s}/5$.

Second, for $1<r/a<\sqrt{8}$ the discretization artifacts are too large to be ignored. The data can still be used after correcting for the discretization artifacts through the nonperturbative improvement. The $\chi^{2}/{\rm d.o.f.}$ is about the same, unless data at large distances are discarded. In this case, the $\chi^{2}/{\rm d.o.f.}$ becomes even smaller, due to the conservative error estimates for the corrections. We do not observe any larger deviation from the results for $r/a\geq\sqrt{8}$, and, hence, estimate that the uncertainty due to the corrections is generally about the same as for $r/a\geq\sqrt{8}$, namely, $\delta^{\rm lat}=\pm 0.00021$. The discretization errors are safely contained in the statistical uncertainty. The results for fits in four representative fit intervals555 We considered many more, but restrict the discussion to a representative set of fits. using nonperturbatively corrected discretization errors are summarized in Tab. 3. In each fit interval all results are consistent. We find slightly lower central values for fits with only 10 data or less, which, however, are never statistically significant. We show the tree-level or nonperturbatively improved lattice data together with the three-loop with leading ultra-soft resummation result with standard scales in Fig. 5. Given the large variation over the $r$ range fine details are difficult to resolve. It can be seen that the fit to the data with large $r/a$ misses the data with $r/a=1$, i.e., the first data for each $\beta$, unless the discretization artifacts are taken care of with an approach beyond the tree-level improvement.

We show in Fig. 6 the differences between the normalized lattice data and the three-loop with leading ultra-soft resummation result with standard scales, either using nonperturbatively improved data, or data with only tree-level improvement. In the considered range the nonperturbative improvement is irrelevant.

We show the differences between the normalized, nonperturbatively improved lattice data and the three-loop with leading ultra-soft resummation result with standard scales in Fig. 7 for fits with $r/a\geq 1$. The difference is generally smaller than $0.007/r_{1}$ (about 4 MeV) and fluctuates around zero with no apparent trend. At distances smaller than $0.15\,r_{1}$, the difference tends to be marginally below zero for $\max(r)=0.055\,{\rm fm}$ and slightly above zero otherwise. The central value of $\alpha_{s}$ is $0.00027$ lower for fits restricted to the shortest distance interval ($\max(r)=0.055\,{\rm fm}$) than for $\max(r)=0.131\,{\rm fm}$. This difference may be indicative of the imperfect removal of the discretization artifacts, but is still compatible with the estimate of the discretization error using $\min(r/a)=\sqrt{8}$, and is more than covered by statistical uncertainties.

We perform another set of tests to investigate the systematic uncertainties of the lattice result. Namely, we test for a bias due to combining multiple lattice spacings into a single analysis that assumes that the data are within uncertainties consistent with the continuum limit. We separately analyze the data at the finest lattice spacing, $\beta=8.4$. Here we also explore whether using data with $r/a\geq\sqrt{5}$ or $\sqrt{8}$ with errors that are inflated by estimates of the discretization artifacts has a significant impact on our results. We use the same uncertainty estimate for the influence of the discretization artifacts that we obtained previously, since it is larger than any estimates that we would obtain by comparing the results in this set of tests. Table 4 shows that the separate analysis for the finest lattice spacing is consistent within statistical errors or within the even smaller estimate of the uncertainties due to the treatment of discretization artifacts. While the preferred central value in this analysis appears to be marginally higher than in the analysis that combines multiple spacings, this is not relevant for our final result.

In order to estimate the perturbative error of the small $\alpha_{s}$ expansion, i.e., the error made by neglecting contributions beyond three loops, we proceed similarly to Ref. Bazavov et al. (2014a) using three methods to quantify the uncertainty of the perturbative calculation. Namely, we add and subtract a term with the typical size of the soft four-loop contribution, we analyze the dependence on the soft scale $\nu$, and we use the difference between the three-loop order with leading ultra-soft resummation and the three-loop result. We estimate the perturbative error in each method as the difference between $\alpha_{s}(M_{Z})$ calculated in that way and the $\alpha_{s}(M_{Z})$ value obtained from the three-loop with leading ultra-soft resummation expression at standard scales. In Ref. Bazavov et al. (2014a), the error due to the soft scale variation was estimated by taking $\nu=\sqrt{2}/r$ and $\nu=1/\sqrt{2}r$. In this study, we enlarge the range of scale variation by considering $\nu=2/r$ and $\nu=1/(2r)$. In doing so, we realize that the $\nu$ dependence of the extracted $\alpha_{s}$ is not monotonous: it has a minimum close to $\nu=1/\sqrt{2}r$, see Fig. 8. For the lowest soft scale $\nu=1/(2r)$ the use of perturbation theory becomes inadequate much earlier than for other choices. Namely, for $\max(r)\gtrsim 0.1\,{\rm fm}$, the fits do not describe the data in a satisfying manner, and $\chi^{2}/{\rm d.o.f.}$ increases much more rapidly for $\nu=1/(2r)$ than for the other choices of the soft scale. This is expected on the grounds of the observation that the 2014 analysis showed that the scale $\nu=1/(\sqrt{2}r)$ was adequate up to about $r=r_{1}/2$. If we exclude those larger distances the upper perturbative error is always given in terms of the soft scale variation to $\nu=2/r$. We observe that for $\max(r)\gtrsim 0.09\,{\rm fm}$ the soft scale variation to $\nu=1/(\sqrt{2}r)$ contributes to a larger estimate of the lower perturbative error than the soft higher-order terms, in line with the 2014 results. However, for smaller distances $\max(r)\lesssim 0.09\,{\rm fm}$ the estimate of the lower perturbative error due to inclusion of soft higher-order terms exceeds the estimate of the lower perturbative error due to soft scale variation. We take the most extreme among the error estimates from the soft contribution as the upper and lower perturbative errors due to the soft contribution. We also estimate the perturbative error by evaluating the difference between the three-loop with leading ultra-soft resummation result and the fixed-order three-loop result. The fixed-order three-loop result always yields a lower central value than the three-loop with leading ultra-soft resummation result. In order to be conservative, we apply the difference as a symmetric error. This has an even larger effect on the lower perturbative error than the soft scale variation or the soft higher-order terms. Moreover, the soft contribution to the lower perturbative error decreases for smaller values of $\max(r)$ more rapidly than the ultra-soft contribution. We stress that the lower perturbative error is always smaller than the statistical error for $r\lesssim 0.1\,{\rm fm}$. We add this symmetric perturbative error due to (inclusion or not of) ultra-soft resummation to the perturbative error due to the soft contribution in quadrature. The reasoning for adding these errors in quadrature is that it is unknown whether there is a similar partial compensation between soft and ultra-soft contributions at higher orders as it happens at three loops, and as such, we have to treat the corresponding uncertainties as independent. Lastly, we compare the three-loop with leading ultra-soft resummation and the two-loop results (first three lines of Eq. (III)). As can be read off from Tabs. 2 and 3, the difference to the two-loop result never exceeds $+0.00025$ and decreases for smaller values of $\max(r)$, i.e., it is smaller than the statistical errors and smaller than the other effects due variation of soft scale, soft higher order terms, or variation of the ultra-soft resummation. The $\chi^{2}/{\rm d.o.f.}$ does not change significantly between using the two-loop or three-loop with leading ultra-soft resummation results. Hence, we confirm the criterion for having the lattice data in the perturbative regime. We observe that the smaller $\max(r)$ is, the smaller the variation of the central value of $\alpha_{s}$ between fits with different forms of the weak-coupling results becomes. Tab. 3 shows clearly that, for a given $\min(r/a)$, the perturbative errors are dramatically reduced at smaller distances, as expected, while the statistical error increases as less data are used to constrain the fits.

Let us summarize the considerations of the preceding paragraphs. We have to use $\max(r)\lesssim 0.1\,{\rm fm}$ to perform the full scale variation and keep the perturbative uncertainties fully under control. We should ideally use significantly more than 10 data point to limit the impact of the imperfectly treated discretization artifacts. Given the considerations of the preceding paragraphs, we take the result for $1\leq r/a\leq 5$ and $\max(r)=0.073\,{\rm fm}$, namely, $\alpha_{s}(M_{Z})=0.11660$ as our final result, which corresponds to $r_{1}\Lambda_{\rm\overline{MS}}^{N_{f}=3}=0.4943$. The uncertainty of the scale $r_{1}$ is $\pm 0.0017\,{\rm fm}$, which yields an error of $\delta^{\rm scale}=\pm 1.7\,{\rm MeV}$ for $\Lambda_{\rm\overline{MS}}^{N_{f}=3}$, and $\delta^{\rm scale}=\pm 0.00010$ for $\alpha_{s}(M_{Z},N_{f}=5)$. Therefore, the final result and full error budget of our zero temperature lattice calculation are given as

[TABLE]

or in terms of $\Lambda_{\rm\overline{MS}}^{N_{f}=3}$ as

[TABLE]

We have added the statistical error and the lattice discretization error of the static energy, the total error of the $r_{1}$ scale, and the perturbative error in quadrature.

In order to compare the current analysis to the previous analysis Bazavov et al. (2014a), we use the smaller window $[1/(\sqrt{2}r),\sqrt{2}/r]$ for the variation of the soft scale $\nu$, and do not account for the uncertainty arising from the difference between resumming or not the leading ultra-soft logarithms to obtain

[TABLE]

In this case, the perturbative error is not dominant anymore. Thus, the presented analysis has approximately halved the uncertainties of Ref. Bazavov et al. (2014a). Nevertheless our final errors are only 10% (upper error) and 30% (lower error) smaller than the ones in Bazavov et al. (2014a), since we have accounted for the other possible sources of uncertainty listed above. The central values of Eq. (4) and of the final result in Ref. Bazavov et al. (2014a) coincide.

IV Extracting $\alpha_{s}$ from the singlet free energy

In this section, we consider the extraction of the strong coupling from the singlet free energy at non-zero temperature, as it is expected that at small distances medium effects are small. We define the singlet free energy in terms of the correlation function of two thermal Wilson lines in Coulomb gauge

[TABLE]

At distances much smaller than the inverse temperature $rT\ll 1$, we can write using pNRQCD Berwein et al. (2017)

[TABLE]

where $\mu_{us}$ is the ultra-soft scale. The form of the thermal correction depends on the scale hierarchy. One could consider the case $1/r\gg\alpha_{s}/r\gg T\gg m_{D}\sim gT$ or the case $1/r\gg T\gg m_{D}\sim gT\gg\alpha_{s}/r$. In the former case $\mu_{us}\sim\alpha_{s}/r$ and $\delta F_{S}(r,T,\mu_{us})=\delta E_{US}(\mu_{us})+\Delta F_{S}(r,T)$ with $E_{US}(\mu_{us})$ being the ultra-soft contribution to the static energy in the vacuum. In the latter case $\mu_{us}\sim T$ and $\delta F_{s}(r,T,\mu_{us})$ has been calculated to order $g^{5}$, i.e., see Eqs. (16) – (19) in Ref. Bazavov et al. (2018b). In the latter case the cancellation of the ultra-soft factorization scale dependence cannot be verified because of the unknown $g^{6}$ contribution to $\delta F_{s}(r,T,\mu_{us})$. Since $V_{s}(r,\mu_{us})$ has a term $\sim\alpha_{s}^{3}\ln(\mu_{us}r)/r$, however, the difference between the $T=0$ static energy and singlet free energy, $E(r)-F_{S}(r,T)$ should have a term $\sim\alpha_{s}^{3}\ln(rT)/r$. This complicates the extraction of the strong coupling from the singlet free energy. The matching between NRQCD and pNRQCD also induces a term $\sim g^{6}T$ in $F_{S}(r,T)$ for both scale hierarchies Berwein et al. (2017), which also needs to be considered.

The singlet free energy has been studied on the lattice in Ref. Bazavov et al. (2018b) using a wide temperature range and several lattice spacings, i.e., several temporal extents $N_{\tau}$. The shortest distance that we can access, due to a single lattice spacing on our finest lattice at $T>0$, is $0.00814\,{\rm fm}$. For our analysis the relevant data correspond to $N_{\tau}=10,~{}12$ and $16$, since $rT$ has to be small. From the analysis of Ref. Bazavov et al. (2018b) we know that thermal effects are small for $rT\lesssim 0.3$. To understand the temperature dependence of the singlet free energy in more detail we show the lattice results for the difference666 The discretization artifacts between both quantities cancel exactly at tree level. Nonetheless, we still use the tree-level improved distance in Fig. 9, in order to permit the visual distinction between data that are inequivalent in the full QCD result. Moreover, we matched the perturbative result to the lattice data at $r/a=1$ using a constant that mimics the effect of the unknown $g^{6}T$ term.

between the singlet free energy at $T>0$ and the static energy at $T=0$ for $\beta=8.4$ (corresponding to the finest zero temperature lattice) in Fig. 9. For other $\beta$ values the results are similar.

From the figure we see that for $r/a\leq\sqrt{6}$ the difference approaches a constant proportional to the temperature according to the above expectations and no temperature effects beyond this constant can be seen in this range within the errors of the lattice results. Therefore, we conclude that in this regime the scale hierarchy $1/r\gg\alpha_{s}/r\gg T$ is appropriate. We treat the finite temperature data in this range as if at zero temperature and fit them with the three-loop with leading ultra-soft resummation result for the static energy with standard scales. Alternatively one could use the $N_{\tau}=12$ data in the range where thermal effects appear to be small together with the two-loop result for the static energy, where the problem of the US scale dependence does not enter at all. This, however, would lead to a larger theoretical uncertainty.

For distances $r/a>\sqrt{6}$ we see some temperature dependence. For the lattice data with $N_{\tau}=12$ this temperature dependence is to some extent captured by the known $g^{5}$ result for $\delta F_{S}(r,T)$, but for $N_{\tau}=10$ and $16$ it is off. On the one hand, for $N_{\tau}=16$, i.e., for a lower temperature, the temperature effects are larger than what is expected based on the $g^{5}$ result. Since the deviation of the temperature effects from a constant or from the $g^{5}$ result is quite similar to the size of the typical discretization artifacts or the statistical errors of the data for $r/a>\sqrt{6}$, it is unclear whether this is truly an effect due to the finite temperature. On the other hand, for $N_{\tau}=10$, where the same temperature corresponds to a coarser lattice, the temperature effects are in the opposite direction and may be caused by temperature-dependent discretization errors Bazavov et al. (2018b).

The known $g^{5}$ result consists of two contributions with opposite sign but similar magnitude for the temperatures under consideration, and, therefore, the overall temperature effect is small due to cancellations. As these two contributions are due to either non-static or static Matsubara modes at orders $g^{4}$ or $g^{5}$, respectively, it is clear that this is an accidental cancellation in the temperature window of our lattice simulations. Furthermore, there is no reason to expect that discretization artifacts in both contributions are similar enough to achieve a cancellation between them. In fact, it has been shown that the $\mathcal{O}(a^{2})$ discretization errors appear to be more pronounced in the difference $E-F_{S}$ than in $E$ or $F_{S}$ individually Bazavov et al. (2018b).

For these reasons, we conclude that for $r/a>\sqrt{6}$ we do not have sufficient understanding of thermal effects to use $F_{S}$ for extraction of $\alpha_{s}$ with fully controlled uncertainties. Nonetheless, the cancellations that appear to be at work both in the $g^{5}$ result and in the lattice data suggest that a fit of the lattice data with the zero temperature result may still be possible as a cross-check. Hence, we again treat the finite temperature data in this range as if at zero temperature and fit them with the three-loop with leading ultra-soft resummation result with standard scales.

In order to determine $\alpha_{s}$ from the singlet free energy at $T>0$ we proceed as follows. We analyze the $T=0$ static energy result ($N_{\tau}=64$) in the same $r/a$ intervals for which we expect that temperature effects in the singlet free energy are under control or are small due to accidental cancellations. We use the same weak-coupling result, namely, the three-loop with leading ultra-soft resummation result with standard scales to obtain $\alpha_{s}$ from the singlet free energy. We estimate the uncertainty due to discretization artifacts to be the same as in the zero temperature analysis, i.e., $\delta^{\rm lat}=\pm 0.00021$. Similarly, we estimate the uncertainty due to temperature effects from the difference of the central values of $\alpha_{s}$ for any combination of $N_{\tau}$ values in each fit interval.

First, we perform fits in the range where temperature effects are expected to be constant and the scale hierarchy $1/r\gg\alpha_{s}/r\gg T\gg m_{D}\sim gT$ applies. We report the result of representative fits with $r/a\leq 2$ in Tab. 5, which corresponds to $rT\lesssim 0.17$ or $0.13$ for $N_{\tau}=12$ or $16$, respectively. The central values of $\alpha_{s}$ in fits to the singlet free energy with $N_{\tau}=12$ or $16$ bracket the result for the static energy at zero temperature, with deviations being smaller than any other uncertainties. Hence, we conclude that we do not resolve nonconstant $T>0$ effects as expected. As for the $T=0$ analysis, we see a marginal decrease of the central value of $\alpha_{s}$ upon restriction to shorter distances. This decrease is about the same magnitude as in the $T=0$ analysis, but is delayed to significantly shorter distances, which supports our interpretation that it is due to imperfections of the nonperturbative improvement procedure. Within the smallest statistical uncertainties or within the even smaller estimates of the lattice discretization uncertainties, all of these results are consistent with the zero temperature analysis of the previous section, while the individual estimates of perturbative errors are systematically smaller, compare with Tab. 3.

Next, we perform fits in the range where temperature effects beyond a constant appear to be small due to accidental cancellations. We assume the same scale hierarchy $1/r\gg\alpha_{s}/r\gg T\gg m_{D}\sim gT$, noting that it is at present not verifiable if the ultra-soft factorization scale dependence cancels in this case. Namely, we perform the fits with $r/a\leq 3$ and report the results in Tab. 5 as well777 Using an even larger fit range up to $r/a\leq\sqrt{12}$ leads to the same conclusions, i.e., the nonconstant thermal effects are still numerically irrelevant. . This corresponds to $rT\lesssim 0.25$ or $0.19$ for $N_{\tau}=12$ or $16$ respectively. $\chi^{2}/{\rm d.o.f.}$ is practically unchanged for $N_{\tau}=12$, but increases slightly for $N_{\tau}=16$. Remarkably, the agreement with the previous $T=0$ analysis, using up to $r/a\leq 5$, does not change significantly as $\max(r)$ becomes larger. This suggests that the differences between the results are not caused by effects of the finite temperature, but rather by the more severe influence of the imperfections of the nonperturbative improvement procedure in the $T=0$ analysis with $r/a\leq 2$ or $r/a\leq 3$. Hence, the estimate of the $\delta^{T>0}$ errors for $r>0.05\,{\rm fm}$ cannot be separated from the $\delta^{\rm lat}$ errors due to lattice discretization uncertainties.

Let us summarize the considerations of the preceding paragraphs. For $T>0$ we have to use $r/a\lesssim 2$ to guarantee the cancellation of ultra-soft factorization scale dependence, although we do not see any indication that it does not work for somewhat larger distances, i.e., $r/a\leq 3$. The weak-coupling picture suggests that cancellation between different medium effects is responsible for the small thermal modification of the singlet free energy. While the medium effects according to the weak-coupling picture are not unambiguously resolved in the data, the data seem to be affected by an accidental cancellation in the temperature window under consideration. Perturbative uncertainties are dramatically reduced at smaller distances $\max(r)$. In Fig. 10 we show our final $T=0$ result with the reported $\chi^{2}/{\rm d.o.f.}$ for two different fits to the singlet free energy data. The data show no nontrivial thermal effects in the ranges considered and are consistent with the central value of the $T=0$ result.

Given the considerations of the preceding paragraphs, we take the $N_{\tau}=12$ result for $1\leq r/a\leq 2$ and $\max(r)=0.030\,{\rm fm}$, namely, $\alpha_{s}(M_{Z})=0.11638$ as our final result, which corresponds to $r_{1}\Lambda_{\rm\overline{MS}}^{N_{f}=3}=0.4900$. Furthermore, it has a largely independent error budget. We note that the finite temperature result for these distances does not depend on gauge ensembles corresponding to the larger sea quark mass, i.e., uses exclusively ensembles with $m_{l}=m_{s}/20$. The scale uncertainty is the same as in the $T=0$ analysis, $\delta^{\rm scale}=\pm 1.7\,{\rm MeV}$ for $\Lambda_{\rm\overline{MS}}^{N_{f}=3}$, and $\delta^{\rm scale}=\pm 0.00010$ for $\alpha_{s}(M_{Z},N_{f}=5)$. Therefore, the final result and full error budget of our finite temperature lattice calculation are given as

[TABLE]

or in terms of $\Lambda_{\rm\overline{MS}}^{N_{f}=3}$ as

[TABLE]

We have added the statistical error, the lattice discretization error of the static energy, the finite temperature error due to using the singlet free energy, the total error of the $r_{1}$ scale, and the perturbative error in quadrature.

V Discussion

In the sections III and IV, we have presented two extractions of the strong coupling constant from two different observables, the static energy on zero temperature lattices, and the singlet free energy on finite temperature lattices.

Both results, i.e., $\alpha_{s}(M_{Z},N_{f}=5)=0.11660^{+0.00110}_{-0.00056}$, and $\alpha_{s}(M_{Z},N_{f}=5)=0.11638^{+0.00095}_{-0.00087}$, are in excellent agreement. This is despite the facts that they are based on several statistically independent sets of gauge ensembles corresponding to the QCD vacuum or the quark-gluon plasma and have been obtained using two different sea quark masses. The systematic uncertainties due to the lattice discretization and due to the lattice scale are not independent, but do not play a significant role in either of the two error budgets. Moreover, the two results are obtained in different $r$ ranges, $[0.0237\,{\rm fm},0.0734\,{\rm fm}]$, or $[0.0081\,{\rm fm},0.0301\,{\rm fm}]$, which overlap only in a narrow window. As such, we may even consider the perturbative uncertainties of the two analyses as practically independent. We stress that every single fit in our analysis is contained already within the purely statistical errors of either the final $T=0$ or the final $T>0$ results as long as the soft scale is at $\nu=1/r$.

There has been an observation in pure Yang–Mills theory that the approach to the continuum limit of $r_{0}\Lambda_{\rm\overline{MS}}^{N_{f}=0}$ extracted in the ${Q\bar{Q}}$-scheme888 The gauge coupling in the ${Q\bar{Q}}$-scheme is defined as $\alpha_{{Q\bar{Q}}}(1/r)=r^{2}F(r)/C_{F}$. appears to change if data with $\alpha_{{Q\bar{Q}}}^{3}\lesssim 0.01$ are considered Husung et al. (2018). Namely, shorter distances corresponding to smaller $\alpha_{{Q\bar{Q}}}$ seemed to be related to smaller $r_{0}\Lambda_{\rm\overline{MS}}^{N_{f}=0}$ in the continuum limit, with about 10% lower continuum value in the pure Yang–Mills theory analysis. We show such a comparison for our data in Fig. 11. The fits with $\max(r)\lesssim 0.05\,{\rm fm}$ are safely in the $\alpha_{{Q\bar{Q}}}^{3}\lesssim 0.01$ range, and are consistent with either a stable value or a modest drop. To account for possible artifacts in our analysis due to the use of small $r/a$, we have included the systematic error estimate $\delta^{\rm lat}=\pm 0.00021$.

Our result from the static energy supersedes the result of the previous analysis Bazavov et al. (2014a), $\alpha_{s}(M_{Z},N_{f}=5)=0.1166^{+0.0012}_{-0.0008}$, obtained using three gauge ensembles with light quark masses $m_{l}=m_{s}/20$. Although it reproduces the previous central value, it includes more conservative estimates of the known sources of systematic error. Nevertheless it achieves to reduce the uncertainties through the use of finer lattice spacings.

In Ref. Takaura et al. (2019b) the strong coupling was calculated from the static energy using two different analyses of spatially-smeared Wilson loop data on three gauge ensembles with a pion mass of about $300\,{\rm MeV}$, which are both consistent with our results. The presented two-step analysis (obtaining the continuum limit first and then extracting the strong coupling) yielded $\alpha_{s}(M_{Z},N_{f}=5)=0.1166^{+0.0021}_{-0.0020}$, while the second analysis using a single global fit yielded a higher value $\alpha_{s}(M_{Z},N_{f}=5)=0.1179^{+0.0015}_{-0.0014}$.

Our result is slightly smaller than the PDG average Tanabashi et al. (2018) $\alpha_{s}(M_{Z},N_{f}=5)=0.1181(11)$ as well as the recent FLAG average Aoki et al. (2019) $\alpha_{s}(M_{Z},N_{f}=5)=0.11823(81)$. It is also smaller than the result of the ALPHA collaboration Bruno et al. (2017) $\alpha_{s}(M_{Z},N_{f}=5)=0.11852(84)$, which uses step scaling and the Schrödinger functional method. These three values are, however, correlated. Our result is smaller than the results of the HPQCD collaboration using pseudoscalar quarkonium correlators or small Wilson loops McNeile et al. (2010); Chakraborty et al. (2015) with three or four sea quark flavors, namely, $\alpha_{s}(M_{Z},N_{f}=5)=0.1183(7)$, or $\alpha_{s}(M_{Z},N_{f}=5)=0.11822(74)$, respectively. Our result agrees with the very recent lattice determination in terms of pseudoscalar quarkonium correlators Petreczky and Weber (2019), $\alpha_{s}(M_{Z},N_{f}=5)=0.1159(12)$, or the preceding analysis of charmonium correlators Maezawa and Petreczky (2016), $\alpha_{s}(M_{Z},N_{f}=5)=0.11622(83)$, or the analyses of higher quarkonium moments Nakayama et al. (2016) $\alpha_{s}(M_{Z},N_{f}=5)=0.1176(26)$. Moreover, our result is consistent within errors with the recent lattice determinations from the hadronic vacuum polarization Hudspith et al. (2018) $\alpha_{s}(M_{Z},N_{f}=5)=0.1181(27)^{+0.0008}_{-0.0022}$, or from the gauge-fixed gluon propagator in Landau gauge Zafeiropoulos et al. (2019), $\alpha_{s}(M_{Z},N_{f}=5)=0.1172(11)$.

Our result is also consistent with phenomenological estimates based on the bottomonium spectrum in NRQCD Mateu and Ortega (2018), $\alpha_{s}(M_{Z},N_{f}=5)=0.1178(51)$, or quarkonia spectra in pNQRCD Peset et al. (2018), $\alpha_{s}(M_{Z},N_{f}=5)=0.1195(53)$, which have errors comparable to an earlier extraction from radiative bottomonium decays Brambilla et al. (2007), $\alpha_{s}(M_{Z},N_{f}=5)=0.119^{+0.006}_{-0.005}$. Finally, a recent study of the static energy using a subset of our data in a much larger $r$ window, with an assumption similar to the one of Ref. Takaura et al. (2019b) and approximations to the weak-coupling result that prevent the quantification of the uncertainties, finds $\alpha_{s}(M_{Z},N_{f}=5)=0.1168$ Mateu et al. (2019), which agrees with our result.

VI Conclusions

In this paper, we have calculated the static energy at zero temperature and the singlet free energy at finite temperature in 2+1 flavor QCD using the HISQ action, and determined the strong coupling. We improved and extended the previous 2+1 flavor HISQ analysis published in Ref. Bazavov et al. (2014a). We included three finer lattice spacings, overhauled the analysis strategy, treated systematic uncertainties more conservatively, and restricted the analysis to much shorter distances that are deeper in the perturbative regime. Our main results are given in Eqs. (4) to (7). Moreover, for the first time we used the singlet free energy on finite temperature lattices to reach much smaller distances than ever used for the analysis of the static energy, obtaining an independent result with a complementary error budget given in Eqs. (12) to (15). Both results agree without even having to invoke the systematic errors of the weak-coupling calculation.

Both of our results for the central value of $\alpha_{s}$ are lower than a few lattice QCD determinations as well as the PDG average. However, the result from the static energy is compatible within errors. Furthermore, our results agree with some other lattice determinations of $\alpha_{s}$ based on moments of quarkonium correlators Petreczky and Weber (2019). Finally, we stress that the extraction of $\alpha_{s}$ from the static energy is the only lattice extraction based, in the continuum part, on an observable known up to three loops and accurate at order $\alpha_{s}^{4}$. This is due to the fact that the static energy is sensitive to $\alpha_{s}$ already at leading order.

Acknowledgments

This work was supported by U.S. Department of Energy under Contract No. DE-SC0012704. J.S. has been supported by the FPA2016-76005-C2-1-P and FPA2016-81114-P projects (Spain), and the 2017-SGR-929 grant (Catalonia). The simulations have been carried out on the computing facilities of the Computational Center for Particle and Astrophysics (C2PAP) in the project “Static quark correlators in lattice QCD at non-zero temperature” (pr83pu) as well as on the SuperMUC at the Leibniz-Rechenzentrum (LRZ) in the project “Static quark-antiquark energy at zero and finite temperature” (pr48le). We thank the Institute for Nuclear Theory at the University of Washington for its kind hospitality and stimulating research environment. This research was supported in part by the INT’s U.S. Department of Energy grant No. DE-FG02-00ER41132. This research was supported by the Munich Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of excellence “Origins” (www.origins-cluster.de). The lattice QCD calculations have been performed using the publicly available MILC code. The data analysis was performed using the R-base and nlme packages R Core Team (2018); Pinheiro et al. (2018).

Appendix A Gauge ensembles and correlators

In this appendix, we discuss additional ensembles that have not been summarized in Tab. 1, since they are used only for auxiliary calculations, and fits to the primary lattice correlators. Further zero temperature ensembles that have only been used for developing a procedure to treat the discretization artifacts are summarized in Tab. 6. They also have been generated by the HotQCD collaboration Bazavov et al. (2012b, 2014b) using the RHMC algorithm and use 2+1 flavors with the strange quark mass at its physical value and the mass of the two degenerate light quarks at 5% of the strange quark mass.

We have also used the gauge configurations generated for the study of quark number susceptibilities at high temperatures Bazavov et al. (2013) and for the study of color screening Bazavov et al. (2018b) by the TUMQCD collaboration. The RHMC algorithm was used for 2+1 quark flavors with the strange quark mass at its physical value and the mass of the two degenerate light quarks at 5% or 20% of the strange quark mass, i.e., corresponding to the zero temperature ensembles. We summarize the nonzero temperature gauge ensembles used in the determination of the strong coupling in Tabs. 7 and 8. We have increased the statistics of some $N_{\tau}=16$ gauge ensembles specifically for this study and updated the Polyakov loop in Tab. 8.

In Tab. 9 we summarize the necessary cuts to the fit ranges for the extraction of the static energy from the lattice correlators.

Appendix B Discretization artifacts

In this appendix, we briefly discuss the discretization artifacts of the static energy calculated on the lattice. We follow the treatment outlined in Ref. Bazavov et al. (2018b) with a few refinements.

It has been known for a long time that the static energy is affected by discretization artifacts at distances that are only slightly larger than the lattice spacing, i.e., $r/a\gtrsim 1$. The ratio of the static energy on the lattice and in the continuum exhibits a distinct pattern of a few percent to permill level variation around 1, which dies out quickly for larger distances, see Fig. 12. This pattern is quite similar in the tree-level calculation and in the QCD result (with or without sea quarks). Hence, one may partly account for these discretization artifacts by dividing the lattice QCD static energy by the ratio of the lattice and continuum static energies from a tree-level calculation, i.e., assuming one-gluon exchange and no running coupling. This procedure is called tree-level improvement. At tree level the static energy is written as

[TABLE]

with $D_{00}(k)$ being the temporal component of the lattice gluon propagator on the lattice. The results on the lattice static energy at tree level calculated for Wilson and Lüscher–Weisz action are shown in Fig. 12. Writing

[TABLE]

defines the so-called improved distance $r_{I}$. In practice, the tree-level improvement is included in the comparison between continuum and lattice results either as

[TABLE]

or – by introducing a fit parameter $\kappa_{I}(\mathbold r/a,\beta)$ for each distinctive path geometry $\mathbold r/a$ between the quark and antiquark at each lattice spacing – as

[TABLE]

Here we use the first approach.

However, this is not sufficient to remove the discretization artifacts entirely. Ideally, one would replace the tree-level improvement by a one- or two-loop improvement, i.e., calculating the correction factors using lattice perturbation theory at one- or two-loop level, and including the tadpole factors. We do not pursue this approach here as it requires extensive calculation beyond the scope of this study. As one can see from Fig. 12 discretization errors are at few percent level for $r/a\leq 2$. For $r/a=\sqrt{5}$ and $\sqrt{6}$ these errors are below a percent and are even smaller for $r/a\geq\sqrt{8}$. Based on previous studies we expect the same pattern of discretization effects in the interacting theory. We show in Fig. 13 that a similar pattern is observed for the HISQ action after the tree-level improvement. For distances $r/a\geq\sqrt{8}$ the statistical errors of the static energy are larger than the discretization errors.

Hence, there are three options to cope with the situation: The first option is to completely avoid the data with significant residual discretization artifacts. This corresponds to the first analysis strategy discussed in Sec. III.

The second option is to estimate the size of the residual discretization artifacts in the data at small $r/a$ and combine a systematic error estimate of similar magnitude with the statistical errors to reduce the statistical weight of these data.

The first quark-antiquark distance that we can access with distinctive path geometries (ignoring permutations) is $r/a=3$, specifically, either through $\mathbold{r}/a=(3,0,0)$, or $(2,2,1)$. We average these two results as $E^{\rm lat}(r=3a)$, and use a Cornell fit to interpolate between the data for $a=r_{I}(\sqrt{8}a),\ 3a,$ and $r_{I}(\sqrt{10}a)$ to the distances $r_{I}((3,0,0)\,a)=2.9788\,a$ and $r_{I}((2,2,1)\,a)=3.0123\,a$. The difference $\Delta^{\rm lat}(\mathbold r/a)=E^{\rm lat}(\mathbold r/a)-E^{\rm int}(\mathbold r/a)$ between the data and the interpolation at $\mathbold{r}/a=(3,0,0)$, or $(2,2,1)$ already includes the tree-level improvement, i.e., $\Delta^{\rm lat}((3,0,0))-\Delta^{\rm lat}((2,2,1))$ is an estimate for the residual discretization artifacts after the tree-level improvement. On the other hand, $E^{\rm lat}((3,0,0))-E^{\rm lat}((2,2,1))$ is an estimate for the discretization artifacts without the tree-level improvement. The latter is usually about 0.1%, the former is usually one order of magnitude smaller. Since the latter certainly overestimates the residual discretization artifacts after tree-level improvement and the former may be accidentally small due to statistical fluctuations, we proceed using the geometric average of the two. We smooth out fluctuations between these geometric averages for different lattice spacings by demanding that the systematic error estimate has to be at least as large as the mean of the systematic error for the two neighboring lattice spacings. We summarize the systematic error estimates for the different lattices in Tab. 10. The absolute size of the estimate of the residual discretization error is typically less than 0.5% of $1/r_{1}$, namely, less than $3\,{\rm MeV}$. This option is suitable for $r/a\geq\sqrt{5}$ at reasonable accuracy and should be applied for all distances with at most two steps along any of the lattice axes.

Finally, we may estimate the residual discretization artifacts in the data for each individual $r/a$, and correct for these errors before comparing lattice data with the weak coupling results. This is the procedure that has been used in Refs. Bazavov et al. (2014a, 2018b). In this procedure, we define the correction factors

[TABLE]

and the nonperturbatively improved static energy

[TABLE]

which will be compared to the continuum static energy $E^{\rm cont}(r_{I})$. In the following, we spell out how we estimate the residual discretization artifacts, and how we correct for these. We use two different estimates for the continuum $E^{\rm cont}(r_{I})$.

1. a)

For the lattices under consideration and $r/a\geq\sqrt{5}$ the relative discretization artifacts are less than 0.1% after the tree-level improvement. We interpolate the static energy on the finest lattices at such distances to obtain a continuum estimate999We extrapolate to 10% shorter distances, see Ref. Bazavov et al. (2018b). that can be used to take ratios of the QCD static energy on coarser lattices and in the continuum.

Our updated interpolation procedure differs in three ways from the procedure outlined in Ref. Bazavov et al. (2018b). Before the interpolation we add a systematic error estimate of 1.5%, 0.5%, 0.5%, 0.5%, 0.04%, 0.02%, and 0.3% to the data on fine lattices at $r/a=1,\ \sqrt{2},\ \sqrt{3},\ 2,\sqrt{5}$, $\sqrt{6}$, and $\sqrt{12}$, respectively. Furthermore, we modify the data before the interpolation by -1.0%, +0.5%, -0.5%, and -0.3% for $r/a=1,\ \sqrt{3},\ 2$, and $\sqrt{12}$, respectively. These estimates of systematic errors are motivated by our previous analysis of the discretization artifacts at short distances Bazavov et al. (2018b). We split the interpolation with a modified Cornell potential into two intervals, a shorter distance (usually 7 points) and a larger distance interval, whose details are decided dynamically depending on the $\chi^{2}/{\rm d.o.f.}$ of the fit, and switch between the two interpolations in an overlap region. In each interval we require fewer higher order terms than for a single fit of all data. We use a spline interpolation to estimate the model dependence. In the next step, we calculate a weighted average $\braket{K(\mathbold r/a,\beta)}$ of the ratios $K(\mathbold r/a,\beta,\beta^{\rm ref})$ obtained with continuum estimates based on different reference lattices (indicated by $\beta^{\rm ref}$) for each $\mathbold r/a$ and $\beta$ as previously.

Using the static energy at zero temperature restricted to $r/a\leq 8$ we vary the reference lattice down to $\beta^{\rm ref}\geq 7.28$ and obtain correction factors for all but the finest lattice ($\beta=8.4$) down to $\beta\geq 6.664$. Using the singlet free energy with $N_{\tau}=12$ restricted to $r/a\leq\sqrt{12}$, namely, $rT\leq 1/\sqrt{12}\approx 0.289$, we obtain correction factors for all but the finest lattice ($\beta=9.67$) down to $\beta\geq 7.373$. 2. b)

We use the three loop result with the leading ultra-soft resummation as a fit form for a continuum estimate up to $r\leq 0.13\,{\rm fm}$. We restrict to this $r$- range to ascertain that the fit is performed in the perturbative regime, which limits this analysis to fine lattice spacings, i.e., to $\beta\geq 7.596$ in practical terms. This fit form is physically more motivated but needs an input value for $\alpha_{s}$. Before the interpolation we add a systematic error estimate of 1.5%, 0.5%, 0.5%, 0.5%, 0.04%, 0.02%, and 0.3% to the data at $r/a=1,\ \sqrt{2},\ \sqrt{3},\ 2,\sqrt{5}$, $\sqrt{6}$, and $\sqrt{12}$, respectively. In order to minimize the bias towards a final determination of $\alpha_{s}$, we vary the range of the input value $\alpha_{s}$ between $0.1155$ and $0.1175$, and determine for each input the normalization constant through a weighted average. Out of the set of inputs we accept either all of the fits with $\chi^{2}/{\rm d.o.f.}\leq 1$ or the three fits with smallest $\chi^{2}/{\rm d.o.f.}$. We estimate the uncertainty of the correction factor by adding the statistical error and the variation of the continuum estimate with the input value of $\alpha_{s}$ in quadrature.

Using the static energy at zero temperature restricted to $r/a\leq 8$ we obtain correction factors for all lattices down to $\beta\geq 7.596$. Using the singlet free energy with $N_{\tau}=12$ restricted to $r/a\leq\sqrt{12}$, namely, $rT\leq 1/\sqrt{12}\approx 0.289$, we obtain correction factors for all lattices down to $\beta\geq 7.825$.

All four estimates of the correction factors are complementary, since they have different systematic uncertainties. In Fig. 14 we show the corresponding correction factors from Cornell fits using lattice data at $r/a\geq\sqrt{5}$ on the reference lattice, or from fits with the weak-coupling result using $r/a\geq\sqrt{1}$. All four results are in very good agreement.

We extrapolate the $\braket{K(\mathbold r/a,\beta)}$ towards the continuum limit using a boosted coupling $\alpha_{s}^{\rm lat}=10/({4\pi\beta u_{0}^{4}})$. We list the corresponding tadpole factors $u_{0}$ in Tab. 11, which have been obtained from the plaquette at $T=0$, or with $N_{\tau}=12$ or $16$.

For the extrapolation we use

[TABLE]

We try extrapolations with up to $N=3$ terms, but never require the third order term, i.e., we cannot fix its coefficient with small uncertainties. While data are quite insensitive to the second order term for $r/a\geq 2$, it cannot be avoided for $r/a<2$. The coefficients of the first and second order term are anti-correlated. Only in the case of the Cornell fit to lattice data at zero temperature we cannot resolve the second order term successfully. We also tested a zeroth order term, which was never necessary.

For each estimate (Cornell vs weak-coupling fits, zero vs finite temperature) we take the median of the correction factors from fits including data starting at $\min(r/a)=\sqrt{1}$ up to $\min(r/a)=\sqrt{5}$. We estimate the overall error from typical uncertainties, or the spread between results for different $\min(r/a)$, whichever is larger. We show the median of the extrapolations for each fit type as function of $\alpha_{s}^{\rm lat}$ as colored bands in Fig. 14. While the agreement between the median of Cornell fits and the median of weak-coupling fits is not satisfactory for fits to $T=0$ data, the different fit results to $T>0$ data agree very well. Moreover, for the correction factor $\braket{K(r/a<2,\beta)}$, where the data are sensitive to the second order term in $\alpha_{s}^{\rm lat}$, the extrapolation of the results obtained with $T>0$ data also extends to the $T=0$ results at larger values of the coupling, i.e, on coarse lattices. Therefore, in order to obtain our final results, we take the mean of the Cornell fit and the weak-coupling fit to the $N_{\tau}=12$ data, and estimate the uncertainty by the largest among the errors of the two results, or half of the spread between the results. Since this uncertainty still underestimates the spread between $\braket{K(r/a=1,\beta)}$ for different fit types and among different $\beta$, we add a systematic error estimate of $0.003$ in quadrature to the errors of this procedure. We summarize the coefficients corresponding to both fit types in Tab. 12 and show the averaged result as a gray band in Fig. 14. While the corrections are significant and constrained quite well up to $\mathbold r/a=(2,1,1)$, the corrections beyond tree-level are consistent with zero for larger distances.

Appendix C Coefficients of the force at three loops

In this appendix, we summarize the constants appearing in the weak-coupling result of the force at three-loop order in Eqs. (III) and (III). The color factors read

[TABLE]

where $N_{c}$ is the number of colors. The QCD beta function is

[TABLE]

where the first three coefficients read Herzog et al. (2017); Luthe et al. (2017)

[TABLE]

$N_{f}$ is the number of active massless quarks. Higher-order coefficients $\beta_{3}$, and $\beta_{4}$ are also known Herzog et al. (2017); Luthe et al. (2017), but are not reproduced for brevity’s sake. We reproduce the coefficients $\tilde{a}_{i}$ as given in Ref. Garcia i Tormo (2013)

[TABLE]

where the coefficients $a_{1}$ and $a_{2}$ read Fischler (1977); Billoire (1980); Schröder (1999)

[TABLE]

The coefficient $a_{3}$ reads Anzai et al. (2010); Smirnov et al. (2010)

[TABLE]

The three numerical values given above are nowadays known analytically Lee et al. (2016); Lee and Smirnov (2016). Lastly, we reproduce the coefficient $a_{3}^{L}$ as given in Ref. Brambilla et al. (1999)

[TABLE]