Full O(a) improvement in EQCD

TL;DR

This paper achieves full O(a) improvement in EQCD by numerically determining the scalar mass renormalization, enhancing the precision of lattice-continuum matching for studying nonperturbative hot QCD effects.

Contribution

It provides the first complete O(a) matching coefficient for EQCD, enabling more accurate lattice simulations of hot QCD phenomena.

Findings

Determined the unknown scalar mass renormalization coefficient.

Updated the EQCD phase diagram.

Facilitates precision studies of transverse momentum diffusion.

Abstract

EQCD is a 3D bosonic theory containing SU(3) and an adjoint scalar, which efficiently describes the infrared, nonperturbative sector of hot QCD and which is highly amenable to lattice study. We improve the matching between lattice and continuum EQCD by determining the final unknown coefficient in the O(a) matching, an additive scalar mass renormalization. We do this numerically by using the symmetry-breaking phase transition point of the theory as a line of constant physics. This prepares the ground for a precision study of the transverse momentum diffusion coefficient C(qperp) within this theory. As a byproduct, we provide an updated version of the EQCD phase diagram.

Tables

Table 1. Table 1: 3D EQCD parameters for four typical scenarios.

$T$	$n_f$	$x$	$y$
$250\mathrm{MeV}$	3	$0.08896$	$0.452423$
$500\mathrm{MeV}$	3	$0.0677528$	$0.586204$
$1\mathrm{GeV}$	4	$0.0463597$	$0.823449$
$100\mathrm{GeV}$	5	$0.0178626$	$1.64668$

Table 2. Table 2: Results of our five EQCD simulation sets.

$x$	$y_{\mathrm{crit},\mathrm{cont}}$	$\delta y_{3loop}/g_{3d}^{2}a$
$0.0463596$	$0.9293\left(13\right)$	$-0.467\left(19\right)$
$0.0677528$	$0.67627\left(85\right)$	$-0.298\left(10\right)$
$0.08896$	$0.54092\left(76\right)$	$-0.1750\left(74\right)$
$0.13$	$0.4043\left(18\right)$	$-0.037\left(18\right)$
$0.2$	$0.2961\left(15\right)$	$0.004\left(15\right)$

Table 3. Table 3: Covariance matrix of the grand fit.

$\mathrm{cov}(C_i,C_j)$	$C_0$	$C_1$	$C_2$	$C_3$
$C_0$	$0.001700$	$-0.02997$	$0.1101$	$-3.563\cdot {10}^{-8}$
$C_1$	$-0.02997$	$0.5451$	$-2.046$	$1.129\cdot {10}^{-6}$
$C_2$	$0.1101$	$-2.046$	$7.899$	$-1.079\cdot {10}^{-5}$
$C_3$	$-3.563\cdot {10}^{-8}$	$1.129\cdot {10}^{-6}$	$-1.079\cdot {10}^{-5}$	$3.025\cdot {10}^{-5}$

Table 4. Table 4: Simulation parameters for the pure scalar simulation.

$\lambda a$	$N_x{N}_y{N}_z$	total statistics
$1/2$	${16}^3$	$402900$
$1/3$	${24}^3$	$1490120$
$1/4$	${32}^3$	$2487440$
$1/6$	${48}^3$	$2931750$
$1/8$	${64}^3$	$2513280$

Table 5. Table 5: Parameters and results of EQCD simulations.

$g_{3d}^{2}a$	$x_{\mathrm{cont}}$	$N_x{N}_y{N}_z$	$y_{\mathrm{crit},\mathrm{cont}}$	statistics $y_{\mathrm{crit}}$	​$\mathrm{Tr}{\Phi }_{\mathrm{brok}}^2/g_{3d}^2$	$\mathrm{Tr}{\Phi }_{\mathrm{symm}}^2/g_{3d}^2$
$1/6$	$0.0463597$	${72}^2\times 192$	$0.824773\left(45\right)$	$300690$	$14.0477\left(40\right)$	$-0.22234\left(28\right)$
$1/8$	$0.0463597$	${96}^2\times 256$	$0.855935\left(85\right)$	$373030$	$15.4237\left(89\right)$	$-0.23064\left(43\right)$
$1/12$	$0.0463597$	${144}^2\times 384$	$0.88387\left(28\right)$	$37170$	$17.716\left(12\right)$	$-0.23390\left(77\right)$
$1/16$	$0.0463597$	${192}^2\times 512$	$0.89504\left(62\right)$	$1300$	$17.714\left(14\right)$	$-0.23748\left(74\right)$
$1/24$	$0.0463597$	${192}^2\times 512$	$0.9072\left(21\right)$	$1430$	$17.585\left(57\right)$	$-0.2355\left(17\right)$ 444For the sake of feasibility, we did not scale the volume here, accordingly. However, ref. Hietanen:2008tv states that no finite volume effects occur for $\beta =\frac{6}{g_{3d}^{2}a}$ being smaller than the smallest extend of the lattice.
$1/4$	$0.0677528$	${48}^2\times 128$	$0.59100\left(19\right)$	$13110$	$6.7502\left(59\right)$	$-0.11392\left(62\right)$
$1/6$	$0.0677528$	${72}^2\times 192$	$0.62149\left(12\right)$	$150000$	$7.6155\left(72\right)$	$-0.12043\left(55\right)$
$1/8$	$0.0677528$	${96}^2\times 256$	$0.63661\left(11\right)$	$200000$	$7.960\left(14\right)$	$-0.12101\left(78\right)$
$1/12$	$0.0677528$	${144}^2\times 384$	$0.64963\left(29\right)$	$38990$	$8.287\left(15\right)$	$-0.12400\left(80\right)$
$1/16$	$0.0677528$	${192}^2\times 512$	$0.6546\left(11\right)$	$1910$	$8.4103\left(61\right)$	$-0.1228\left(12\right)$
$1/3$	$0.08896$	${36}^2\times 96$	$0.47232\left(24\right)$	$10000$	$4.0728\left(22\right)$	$-0.03277\left(47\right)$
$1/4$	$0.08896$	${48}^2\times 128$	$0.49119\left(14\right)$	$100000$	$4.4343\left(53\right)$	$-0.03489\left(59\right)$
$1/6$	$0.08896$	${72}^2\times 192$	$0.50940\left(15\right)$	$150000$	$4.7728\left(35\right)$	$-0.0344\left(10\right)$
$1/8$	$0.08896$	${96}^2\times 256$	$0.51736\left(20\right)$	$97500$	$4.9249\left(58\right)$	$-0.0336\left(11\right)$
$1/12$	$0.08896$	${144}^2\times 384$	$0.52565\left(45\right)$	$44530$	$5.052\left(10\right)$	$-0.0321\left(13\right)$
$1/3$	$0.13$	${36}^2\times 96$	$0.37184\left(94\right)$	$10000$	$2.3337\left(64\right)$	$0.1112\left(15\right)$
$1/4$	$0.13$	${48}^2\times 128$	$0.38341\left(83\right)$	$100000$	$2.3899\left(64\right)$	$0.1126\left(19\right)$
$1/6$	$0.13$	${72}^2\times 192$	$0.39323\left(21\right)$	$150000$	$2.4703\left(75\right)$	$0.1190\left(38\right)$
$1/8$	$0.13$	${96}^2\times 256$	$0.39702\left(31\right)$	$109230$	$2.5371\left(82\right)$	$0.1246\left(39\right)$
$1/12$	$0.13$	${144}^2\times 384$	$0.3982\left(11\right)$	$120600$	$2.618\left(13\right)$	$0.1332\left(65\right)$
$1/3$	$0.2$	${36}^2\times 96$	$0.27743\left(41\right)$	$12340$	$1.3115\left(74\right)$	$0.578\left(32\right)$
$1/4$	$0.2$	${48}^2\times 128$	$0.28595\left(34\right)$	$21030$	$1.283\left(13\right)$	$0.630\left(34\right)$
$1/6$	$0.2$	${72}^2\times 192$	$0.29214\left(37\right)$	$66160$	$1.2850\left(81\right)$	$0.650\left(24\right)$
$1/8$	$0.2$	${96}^2\times 256$	$0.29377\left(48\right)$	$122080$	$1.323\left(11\right)$	$0.713\left(37\right)$
$1/12$	$0.2$	${144}^2\times 384$	$0.29393\left(97\right)$	$3330$	$1.390\left(14\right)$	$0.707\left(26\right)$

Table 6. Table 6: Continuum-extrapolated Tr ​ Φ 2 g 3 ​ d 2 Tr superscript Φ 2 superscript subscript 𝑔 3 d 2 \frac{\mathrm{Tr}\,\Phi^{2}}{g_{3\mathrm{d}}^{2}} in the symmetric phase at criticality, 𝒪 ​ ( a ) 𝒪 𝑎 \mathcal{O}(a) operator improvement of Tr ​ Φ 2 g 3 ​ d 2 Tr superscript Φ 2 superscript subscript 𝑔 3 d 2 \frac{\mathrm{Tr}\,\Phi^{2}}{g_{3\mathrm{d}}^{2}} and continuum- extrapolated difference of Tr ​ Φ 2 g 3 ​ d 2 Tr superscript Φ 2 superscript subscript 𝑔 3 d 2 \frac{\mathrm{Tr}\,\Phi^{2}}{g_{3\mathrm{d}}^{2}} in the two phases at criticality.

$x$	$\mathrm{Tr}{\Phi }_{\mathrm{symm},\mathrm{cont}}^2/g_{3d}^2$	$\delta \mathrm{Tr}{\Phi }_{3loop}^2/g_{3d}^{4}a$	$\Delta \mathrm{Tr}{\Phi }_{\mathrm{cont}}^2/g_{3d}^2$
$0.0463597$	$-0.2350\left(24\right)$	$-0.081\left(44\right)$	$18.232\left(27\right)$
$0.0677528$	$-0.1231\left(22\right)$	$-0.019\left(30\right)$	$8.692\left(13\right)$
$0.08896$	$-0.0272\left(27\right)$	$-0.072\left(25\right)$	$5.177\left(11\right)$
$0.13$	$0.151\left(11\right)$	$-0.256\left(70\right)$	$2.502\left(24\right)$
$0.2$	$0.763\left(76\right)$	$-0.66\left(85\right)$	$0.689\left(33\right)$

Table 7. Table 7: Covariance matrix of the grand fit of the missing 3 loop δ ​ Φ 2 g 3 ​ d 2 𝛿 superscript Φ 2 superscript subscript 𝑔 3 d 2 \frac{\delta\Phi^{2}}{g_{3\mathrm{d}}^{2}} -contribution.

$\mathrm{cov}(C_i,C_j)$	$C_{\Phi a}$	$C_{\Phi b}$	$C_{\Phi c}$
$C_{\Phi a}$	$21694$	$-267617$	$5.3859\cdot {10}^{-6}$
$C_{\Phi b}$	$-267617$	$3.4914\cdot {10}^6$	$-1.432\cdot {10}^{-4}$
$C_{\Phi c}$	$5.3859\cdot {10}^{-6}$	$-1.432\cdot {10}^{-4}$	$9.0\cdot {10}^{-4}$

Equations

$$C({\mathbf{b}}_{\perp })=\underset{L\to \infty }{\lim }-\frac{1}{L}\ln W[(0,0,0);(0,0,{\mathbf{b}}_{\perp });(L,L,{\mathbf{b}}_{\perp });(L,L,0)],$$

$$S_{\mathrm{EQCD,c}}=\int\mathrm{d}^{3}x\,\left(\frac{1}{2g_{3\mathrm{d}}^{2}}\mathrm{Tr}\,F^{ij}F^{ij}+\mathrm{Tr}\,D^{i}\Phi D^{i}\Phi+m_{\mathrm{D}}^{2}\mathrm{Tr}\,\Phi^{2}+\lambda\big{(}\mathrm{Tr}\,\Phi^{2}\big{)}^{2}\right)\,.$$

$$\int {\mathrm{d}}^{3}x\mathrm{Tr}{D}^i\Phi D^i\Phi \to 2Z_{\Phi }\sum _{\mathbf{x},i}\mathrm{Tr}\left({\Phi }_{\mathrm{L}}^2(\mathbf{x})-{\Phi }_{\mathrm{L}}(\mathbf{x})U_i(\mathbf{x}){\Phi }_{\mathrm{L}}(\mathbf{x}+a\hat{i})U_i^{†}(\mathbf{x})\right)$$

$$S_{\mathrm{EQCD},\mathrm{L}}=$$

$$Z_2=$$

$$\delta y=$$

$$\displaystyle-\frac{1}{16\pi^{2}}\left[10x\bigg{(}\frac{3}{2}\Sigma^{2}+3\Sigma\xi-6\delta\bigg{)}+\bigg{(}\ln\bigg{(}\frac{6}{g_{3\mathrm{d}}^{2}a}\bigg{)}+\zeta-3\Sigma\xi\bigg{)}\big{(}60x-20x^{2}\big{)}\right.$$

$$\displaystyle\left.\hskip 40.00006pt+9\bigg{(}\frac{7}{8}\Sigma^{2}-\frac{\Sigma\pi}{6}+\frac{31\Sigma\xi}{6}+2\kappa_{1}-\kappa_{4}+4\rho-4\delta\bigg{)}\right]\,,$$

$$\delta y_{3\mathrm{loop}}=g_{3\mathrm{d}}^{2}a\Big{(}C_{0}+C_{1}x+C_{2}x^{2}+C_{3}x^{3}\Big{)}.$$

$$y(z)=y_{\mathrm{crit},\mathrm{est}}+\Delta y\cos (2\pi z/L_z),$$

$$y_{\mathrm{L},\mathrm{new}}=y_{\mathrm{L},\mathrm{old}}+c_{\mathrm{B}}\cdot \left(\frac{\frac{1}{V}{\sum }_{\mathbf{x}}\mathrm{Tr}{\Phi }^2-\mathrm{Tr}{\Phi }_{\mathrm{symm}}^2}{\mathrm{Tr}{\Phi }_{\mathrm{brok}}^2-\mathrm{Tr}{\Phi }_{\mathrm{symm}}^2}-0.5\right),$$

$$y_{\mathrm{crit}}\left(g_{3\mathrm{d}}^{2}a\right)=\sum _{j=0}^{N_{\mathrm{poly}}}{y}_j(g_{3\mathrm{d}}^{2}a{)}^j.$$

$$|y_i|\le \frac{y_0}{2^i},$$

$$C_3=\frac{\delta {\tilde{y}}_{3loop}}{\lambda a}=0.0151(55).$$

$$\frac{\delta y_{3\mathrm{loop}}}{g_{3\mathrm{d}}^{2}a}(x)=0.0151(55)x^3-31.8(28)x^2+10.80(74)x-0.886(41),$$

$$x{y}_{\mathrm{crit}}=\frac{3}{8{\pi }^2}$$

$$\mathcal{B}\equiv \frac{⟨{\left(\mathrm{Tr}{\bar{\Phi }}^2\right)}^2⟩}{⟨\mathrm{Tr}{\bar{\Phi }}^2{⟩}^2},$$

$$\frac{\mathrm{Tr}{\Phi }_{\mathrm{cont}}^2}{g_{3\mathrm{d}}^2}=Z_m\left(Z_{\Phi }\mathrm{Tr}{\Phi }_{\mathrm{L}}^2-\frac{\delta {\Phi }^2}{g_{3\mathrm{d}}^2}\right),$$

$$Z_m^{-1}=1+\frac{g_{3\mathrm{d}}^{2}a}{4\pi }\left(3N\xi +\frac{N\Sigma }{6}-(N^2+1)\xi x\right).$$

$$\frac{\delta {\Phi }^2}{g_{3\mathrm{d}}^2}$$

$$C_{\Phi c}=2(N^2+1)C_4,C_4=0.5630(4),$$

$$C_{\Phi a}$$

$$C_{\Phi b}$$

Full text

1 Introduction

At low energy scales, and therefore at low temperatures, the coupling of QCD becomes large and the theory’s behavior becomes nonperturbative. Therefore we should not be surprised if perturbation theory does not work for thermodynamical or dynamical properties as one approaches the QCD crossover temperature, $T\sim 150\>\mathrm{MeV}$ Bazavov:2014pvz ; Borsanyi:2013bia , from above. However, it came as a surprise just how badly perturbation theory works at scales up to many times the transition temperature. For instance, thermodynamical properties such as the pressure have an expansion in $g$ the strong coupling which is known up to $g^{6}\ln(g)$ Shuryak:1977ut ; Kapusta:1979fh ; Toimela:1983407 ; Arnold:1994eb ; Zhai:1995ac ; Kajantie:2002wa , and while the leading-order behavior is within $30\%$ of the lattice result above $360\>\mathrm{MeV}$ Bazavov:2014pvz ; Borsanyi:2013bia , the series of corrections does not converge even for $T=100\>\mathrm{GeV}$, a scale where perturbation theory should work well Laine:2003ay . This problem was first understood broadly by Linde Linde:1980ts , and was diagnosed more completely starting in the mid-1990s with the work of Braaten and Nieto Braaten:1995cm , who showed that the perturbative expansion could be understood as a two-step process. Treating the problem in Euclidean space, the time direction is periodic with period $\beta=1/T$, corresponding in frequency space to a tower of discrete frequencies $\omega_{n}=2\pi nT$, the Matsubara frequencies (fermions have $\omega_{n}=(2n+1)\pi T$). One can integrate out all but the $n=0$ modes of the spatial gauge field $A_{i}$ and its temporal component $A_{0}$: dimensional reduction Appelquist:1981vg ; Nadkarni:1982kb . This results in a 3-dimensional theory of an SU(3) gauge field and an adjoint scalar ($A_{0}$ behaves as a scalar and we will call it $\Phi$ henceforth), which has been christened “3D Electrostatic QCD” or EQCD for short.

Explicit loopwise-order calculations of the matching between full thermal QCD with quarks and EQCD Kajantie:1997tt ; Laine:2005ai indicate that integrating out the nonzero Matsubara frequencies is a well-behaved perturbative procedure down to temperatures of order $300\>\mathrm{MeV}$. It is the behavior of EQCD itself which is not under perturbative control. But one can solve EQCD nonperturbatively on the lattice, and this appears to generate much closer approximations to the full 4D thermodynamics than perturbation theory alone Laine:2005ai ; Hart:2000ha ; Hietanen:2008xb .

EQCD is a 3D theory of bosons only, which is relatively easy to treat on the lattice; for most thermodynamical quantities, high-quality results were already available 20 years ago. But in 2008 Caron-Huot showed CaronHuot:2008ni that EQCD is also the right effective theory for computing a dynamical transport property, $C(q_{\perp})$ the differential rate at which a highly relativistic colored particle exchanges transverse momentum with the medium Bass:2008rv ; CasalderreySolana:2007qw . This has important applications in hard particle suppression and jet modification in the hot QCD medium. For instance, the frequently discussed “transport coefficient” $\hat{q}$ is defined as the $q_{\perp}^{2}$ moment of $C(q_{\perp})$, $\hat{q}=\int\frac{d^{2}q_{\perp}}{(2\pi)^{2}}q_{\perp}^{2}\>C(q_{\perp})$. To calculate $C(q_{\perp})$, one notes that its transverse-space Fourier dual $C({\boldsymbol{b}}_{\perp})$ is defined in terms of a Wilson loop with two spatial and two lightlike edges CasalderreySolana:2007qw :

[TABLE]

where $W(a;b;c;d)$ is the (fundamental representation) Wilson loop connecting the four spacetime points $(a;b;c;d)$, and we have written $a=(a_{t},a_{z},\boldsymbol{a}_{\perp})$. What Caron-Huot showed is that $W$ can be replaced with a similar Wilson loop in EQCD, with the $L$-length edges replaced by a modified Wilson line which also incorporates the scalar field $\Phi$, up to corrections which should be highly suppressed in any regime where dimensional reduction works. It would be extremely interesting to pursue a high-quality determination of $C(q_{\perp})$ within EQCD, to extract directly important nonperturbative dynamical information about the behavior of hot QCD, of direct relevance for experiment.

Panero, Rummukainen, and Schäfer have made a first exploration of this quantity in EQCD Panero:2013pla . However it appears that the numerics are much more challenging than for standard thermodynamical quantities, especially for small ${\boldsymbol{b}}_{\perp}$ values. Therefore it is essential to minimize lattice spacing errors, since the numerical cost of a study increases as approximately $a^{-5}$, with $a$ the lattice spacing. In a 3D theory it should be possible to perform a lattice-continuum matching which is free of $\mathcal{O}(a)$ errors. Without doing so, the leading errors are of order $\mathcal{O}(a/b_{\perp})$; with $a$-errors removed, the leading errors are $\mathcal{O}(a^{2}/b_{\perp}^{2})$, a very significant improvement for small $b_{\perp}$. At the level of the Lagrangian parameters of EQCD, such an improvement was performed a long time ago Moore:1997np – except that one parameter remains undetermined at the $\mathcal{O}(a)$ level. When evaluating a particular composite operator such as the modified Wilson line, one must also determine the $\mathcal{O}(a)$ corrections to the operator; but this was done for the modified Wilson line operator a few years ago DOnofrio:2014mld . So we would have everything we need to proceed with a lattice study, free from any $\mathcal{O}(a)$ errors, if we only knew the lattice-continuum matching of that one remaining Lagrangian parameter. Specifically, while most Lagrangian terms receive $\mathcal{O}(a)$ corrections at one loop, the $\Phi$-field mass-squared parameter receives (known) $\mathcal{O}(1/a)$ corrections at one loop and $\mathcal{O}(\ln(a))$, $\mathcal{O}(1)$ corrections at two loops Kajantie:1997tt . The unknown $\mathcal{O}(a)$ errors arise at the 3-loop level DOnofrio:2014mld .

This paper will undertake the necessary technical development of determining this $\mathcal{O}(a)$ correction in the lattice-continuum matching of 3D EQCD theory. The direct diagrammatic evaluation appears too difficult to pursue; so we will use an alternative method to extract the corrections. Since the $m^{2}_{3\mathrm{d}}$-renormalization is the only missing $\mathcal{O}(a)$-contribution, it can be determined by fitting to a line of constant physics. EQCD features a phase transition, which we will utilize to obtain such a line.111Note that the phase transition in EQCD is unphysical in the sense that it is not related to any phase transition in 4D thermal QCD. Indeed, 4D QCD with physical quark masses has a crossover at zero quark chemical potential Bazavov:2014pvz ; Borsanyi:2013bia ; the pure-glue theory at high temperature has a $Z_{3}$-breaking phase structure which is related to the EQCD phase transition in a 1-loop analysis Kajantie:1997tt , but because EQCD lacks a true $Z_{3}$ symmetry, this is essentially a coincidence (see however Ref. Vuorinen:2006nz ).

In the remainder of the paper, we present our investigation of the matching problem. Section 2 sets the theoretical stage. Section 3 presents our approach to determining the 3-loop mass renormalization indirectly from lines of constant physics. We present our results in Section 4 and leave conclusions and outlook to Section 5. A few odds and ends appear in two appendices.

For the impatient reader, here is a very short summary. The theory EQCD has two parameters, the mass-squared $y$ and scalar self-coupling $x$ (the gauge coupling just sets a scale). For a given value of the self-coupling $x$, there is a critical $y$-value where a phase transition occurs. We find this point at several lattice spacings, and extrapolate to the continuum behavior; the slope of the fit at $a=0$ is precisely the $\mathcal{O}(a)$ mass correction which must be compensated. Perturbative arguments show that the resulting slope should depend on $x$ as a third-order polynomial. Determining this at several $x$-values allows us to fit all polynomial coefficients, which are presented with their covariance matrix in Eq. (13) and Table 3.

2 Theoretical setup

The theory EQCD is a 3-dimensional SU(3) gauge theory with gauge field $A_{i}^{a}T^{a}$ ($i=1,2,3$ and $a=1\ldots 8$) with a real adjoint scalar field $\Phi$ which can be understood as the dimensional reduction of the 4D Euclidean $A_{0}$ field component. The continuum action is

[TABLE]

The parameter $m_{D}^{2}$ has logarithmic scale dependence which we resolve in the same way as in Kajantie:1997tt . We will use the coupling $g_{3\mathrm{d}}^{2}$, which has dimensions of energy, to set the scale, and we work in terms of the dimensionless ratios $x\equiv\lambda/g_{3\mathrm{d}}^{2}$ and $y\equiv m_{D}^{2}(\bar{\mu}=g_{3\mathrm{d}}^{2})/g_{3\mathrm{d}}^{4}$, originally introduced by Farakos:1996705 ; Kajantie:1997tt .

We will not present the full details of our lattice implementation or update algorithms, since they are almost identical to Kajantie:1995kf . We use the standard Wilson gauge action and nearest-neighbor scalar gradient or “hopping” term. The only crucial difference to the presented $SU(2)$ + fundamental scalar-case concerns the treatment of the hopping term in the gauge field update. It arises from the scalar kinetic term, which translates into

[TABLE]

in the lattice formulation, where $\Phi_{\mathrm{L}}$ is the rescaled, dimensionless lattice version of the adjoint scalar field, $Z_{\Phi}$ is a field renormalization factor, and $U_{i}({\boldsymbol{x}})$ is the standard gauge link at lattice site ${\boldsymbol{x}}$ in direction $i$. In contrast to the fundamental scalar case treated in Kajantie:1995kf , the present hopping term is non-linear in the link. Therefore, it has to be incorporated into the link update via a Metropolis step. Our scalar update, on the other hand, is a mixture of heatbath updates with the $x\mathrm{Tr}\,\Phi^{4}$ term included by Metropolis accept/reject, and the overrelaxation update introduced in Ref. Kajantie:1995kf . We update sites in checkerboard order. Our code was modified from the OpenQCD-1.6 package openQCD .

Now we return to the parameters of the continuum and lattice actions. For this choice of parameters, 1-loop relations between the lattice gauge and scalar couplings and their continuum values, and two-loop relations for the scalar mass, are known; we use the expressions from DOnofrio:2014mld .222The paper is written for general gauge groups, where there are two independent scalar self-couplings. These are equivalent in $SU(3)$, so we take $x_{2}=0$ in their notation. Note that in the lattice action in DOnofrio:2014mld , $x_{1}$ and $x_{2}$ actually have to be interchanged for consistency with the rest of that paper. Also, since the normalization of the lattice scalar field is arbitrary, we have chosen $Z_{\Phi}=1$, that is, we normalize our hopping term to have unit norm. The matching between the lattice and continuum is such that we know the lattice $x$ and $g_{3\mathrm{d}}^{2}$ parameters up to $\mathcal{O}(a^{2}g_{3\mathrm{d}}^{4})$ corrections. Effects from higher-dimension operators (present in the Wilson action and nearest-neighbor hopping term) are also of $\mathcal{O}(a^{2})$. We also know the multiplicative rescaling between $y$ and $y_{\mathrm{latt}}$ to the same precision, and we know the $\mathcal{O}(1/a)$ and $\mathcal{O}(1,\ln(a))$ additive contributions to $y$. Only the (3-loop) $\mathcal{O}(a)$ additive contribution to $y$ is unknown. Therefore any $\mathcal{O}(a)$ difference in a physical result between lattice treatments at different lattice spacings must arise due to this additive contribution.

The phase structure of EQCD was extensively examined in the 90’s, for example in Kajantie:1998yc . The theory has a $\mathbb{Z}_{3}$ symmetry which is broken if $\mathrm{Tr}\Phi^{3}$ takes a nonvanishing value. There is a line of phase transitions separating a large-$y$ region, where $\mathbb{Z}_{3}$ symmetry is preserved, from a small-$y$ region where $\mathbb{Z}_{3}$ symmetry is spontaneously broken. Unlike the transition in $SU(2)$ fundamental Kajantie:1995kf or adjoint Kajantie:1997tt theories, this transition line extends over all $x$ values, since the phases are distinguished by a global discrete symmetry breaking. At small $x$ values the transition is first order; there is a tricritical point, and for large $x$ values it is second order Kajantie:1998yc . Values of $x$ corresponding to dimensional reduction from physical temperatures and quark numbers all land in a region where the transition is first order; they also lie below the critical value $y_{\mathrm{crit}}$, so physical QCD corresponds to metastable points in the EQCD phase diagram. (We emphasize again that the phase transition in EQCD is not related to any thermal phase transitions which may or may not occur in 4D QCD.)

Our methodology will consist of determining, for a given $x$ value, the value $y_{\mathrm{crit}}$ where the phase transition occurs. Doing so at a series of lattice spacings provides a lattice determination of the lattice spacing $a$ dependence of $y_{\mathrm{crit}}$. Since the only $\mathcal{O}(a)$ error remaining in our lattice implementation of the theory is an additive shift to $y$, the slope of $y_{\mathrm{crit}}(a)$ when we extrapolate the lattice spacing $a\to 0$ determines the unknown linear-in-$a$ correction to $y$ at each given $x$ value.

Formally, we know that the $\mathcal{O}(a)$ lattice-continuum additive $\delta y$ contribution arises from 3-loop scalar self-energy diagrams in lattice perturbation theory Moore:1997np . Even without computing these graphs, we can see that they involve 0, 1, 2, and 3 factors of the scalar self-coupling. Therefore, making $N=3$ in the expression from DOnofrio:2014mld explicit, we give the lattice mass-squared in terms of the continuum $y$ value,

[TABLE]

where $\zeta=0.08849$, $\delta=1.942130$ and $4\rho-2\kappa_{1}+\kappa_{4}=-1.968325$. The undetermined $\mathcal{O}(a)$ additive contribution must be parametrically of form

[TABLE]

With results at enough $x$ values, we can perform a polynomial fit to extract these coefficients, and use it to determine the correction at any $x$ value.

Eventually we want to apply EQCD to study QCD. Dimensional reduction at a specific temperature (hence gauge coupling) and number of light fermions leads to a specific $x$ and $y$ value. The 2-loop reduction formulae between high-temperature 3+1 dimensional full QCD and EQCD were worked out by Kajantie et al Kajantie:1997tt ; Laine:2005ai and we use a nonperturbative value of $\Lambda_{\overline{\mathrm{MS}}}$ from Bruno:2017gxd . These lead to the specific $x$ and $y$ values, which we will later investigate for $C(q_{\perp})$ behavior, shown in Table 1. To minimize errors in a future investigation, we will examine the mass renormalization at the $x$-values indicated, except the smallest value where our method will prove ineffective. We will also study larger values of $x$ which do not correspond to any physical QCD regime.

The determination of $y_{\mathrm{crit}}$ faces the usual challenges of supercritical slowing down, associated with determining a first order phase transition point numerically. In the next section we present a methodology for evading this problem.

3 Our method

The standard way of determining $y_{\mathrm{crit}}$ would be by applying multicanonical reweighting in order to enforce tunneling between the two phases Kajantie:1998yc . This is rather inefficient, so we develop another approach to efficiently determine the transition temperature of a first-order phase transition on the lattice.

The main idea is to prepare a lattice configuration where the two phases coexist and are permanently compared to each other at the phase boundaries. If we miraculously guessed the exact value of $y_{\mathrm{crit}}$, the symmetric phase volume would change only via Brownian motion. If our value for $y$ were close to but not exactly $y_{\mathrm{crit}}$, the phase boundaries would feel a small net pressure, and would tend to allow the preferred phase to expand at the expense of the other. This leaves us with two questions:

• •

How do we prepare such configurations?

• •

How can we tune the mass to its critical value and balance the Brownian motion of the phase boundaries?

The true order parameter of EQCD is $\mathrm{Tr}\,\Phi^{3}$, which indicates whether the $\mathbb{Z}_{3}$-symmetry of $\Phi$ is present or broken. However, the phase transition can also be spotted in $\mathrm{Tr}\,\Phi^{2}$ (see Fig. 1), which has smaller fluctuations and leads to a more stable phase discriminator; so we use it in the following. Our approach begins by bounding $y_{\mathrm{crit}}$ by performing a simulation in a modest-sized cubic box, starting from a quite positive $y$ value and decreasing it after each update sweep. At some value, $\mathrm{Tr}\,\Phi^{2}$ abruptly jumps. Then one steadily increases $y$ until the value abruptly falls. This determines upper and lower spinodal $y$-values; $y_{\mathrm{crit}}$ must lie between, typically close to the upper value.

Next we estimate $\mathrm{Tr}\,\Phi^{2}_{\mathrm{symm}}$ and $\mathrm{Tr}\,\Phi^{2}_{\mathrm{brok}}$, the values of $\mathrm{Tr}\,\Phi^{2}_{\mathrm{L}}$ in each of these phases at a mass close to the transition temperature, which we do in separate simulations which are initialized with either vanishing or large constant $\Phi$ values. The method will be rather insensitive to the exact values of these quantities, so it is not important if the determinations are from somewhat incorrect $y$ values.

Next we set up our mass tuning algorithm. We work in a rectangular periodic lattice with one long ($L_{z}$) direction and two equal shorter ($L_{x}=L_{y}$) directions. Initially we make the $y$ (Lagrangian) value $z$-coordinate dependent,

[TABLE]

with $\Delta y$ chosen initially to be large enough that $y_{\mathrm{crit,est}}\pm\Delta y$ are above/below the spinodal values. After a series of update sweeps, the field will find the symmetric phase where $y$ is large and the broken phase where $y$ is small, generating our configuration with both phases and two phase boundaries. Then the magnitude of $\Delta y$ is gradually lowered over a series of update sweeps; if one phase starts to win out over the other, the estimated critical value is adjusted.

With our starting two-phase configuration and estimated $y_{\mathrm{crit}}$ in hand, we proceed to the more accurate determination of $y_{\mathrm{crit}}$. We continue to evolve with a space-uniform $y$ value, but we adjust it after each lattice site-update according to

[TABLE]

where $y_{L}\equiv Z_{3}(y+\delta y_{2\mathrm{loop}})$ and $c_{\mathrm{B}}$ is a small coefficient that controls the strength of the adjustment. The quantity in brackets here is an estimate for the fraction of the volume which lies in the broken phase, based on the known (approximate) values of $\mathrm{Tr}\,\Phi^{2}$ in each phase. Therefore, the adjustment term shifts $y$ upwards (making the symmetric phase more preferred) when more volume is in the broken phase, and downwards (making the broken phase more preferred) if more of the volume is symmetric.

The coefficient $c_{\mathrm{B}}$ is small, $\mathcal{O}\left(\frac{1}{N_{\mathrm{x}}N_{\mathrm{y}}N_{\mathrm{z}}}\right)$, such that the evolution of $y$ is as mild as possible, consistent with enough restorative force to prevent either phase from “winning.” Specifically, whenever $y$ deviates from $y_{\mathrm{crit}}$, there is a net force on the interface, equal to the surface area times $\Delta F$ the free energy difference between phases. At $y=y_{\mathrm{crit}}$, $\Delta F=0$ and there is no net force on the interface. Away from $y=y_{\mathrm{crit}}$, we can expand $\Delta F$ in a Taylor series in $y-y_{\mathrm{crit}}$. At leading order in small $y-y_{\mathrm{crit}}$, the free energy difference will be linear in $y-y_{\mathrm{crit}}$, and the central value of $y$ which maintains coexistence will equal $y_{\mathrm{crit}}$. At quadratic order, $d^{2}F/dy^{2}\neq 0$ means that the restorative force is slightly biased, and we will obtain an incorrect value for $y_{\mathrm{crit}}$. We test for such a distortion by performing a second evolution where $c_{\mathrm{B}}$ is twice as large, to confirm that the central value of $y$ is the same within errors (which it is in all cases we considered).

4 Results

We use the procedure described in the previous section to determine the critical value $y_{\mathrm{crit}}(x,a)$ for several values of the scalar self-coupling $x$, each at several lattice spacings. The exact list of lattices considered is given in Table 5. Because our procedure leads to relatively long autocorrelation in the estimated $y_{\mathrm{crit}}$ value, the errors must be determined via the jackknife method using relatively wide jackknife bins; we vary the bin widths until the error estimates stabilize. We then subtract the known 1- and 2-loop contributions and apply the known multiplicative rescalings Moore:1997np from the results and convert $y_{\mathrm{L,\,crit}}\rightarrow y_{2\mathrm{loop},\,\mathrm{crit}}$. For each $x$ value, we must extrapolate this quantity to zero lattice spacing; the intercept is the continuum critical $y$ value and the slope at the intercept is the desired $\mathcal{O}(a)$ additive correction to the scalar mass.

Because our $y_{\mathrm{crit}}(x,a)$ results are quite precise but the $a$ values are not extremely small, we anticipate that $y_{\mathrm{crit}}(x,a)$ contains corrections beyond linear order in $a$. In principle, we could straightforwardly fit a polynomial of order $N_{\mathrm{poly}}$ in $g_{3\mathrm{d}}^{2}a$ as

[TABLE]

However, as often occurs, ever-higher order coefficients are ever less certain, and including too many coefficients tends to overfit the data and artificially inflates the final fitting errors. In order to extract these coefficients as efficiently as possible from the data, we would like to build in our knowledge about the convergence of the perturbative series to the fit. A useful tool to implement this is constrained curve fitting Davies:1994mp ; Lepage:2001ym . Motivated by a rough estimate of the radius of convergence $(g_{3\mathrm{d}}^{2}a)_{\mathrm{conv}}\approx 0.5$, we make the a priori-guess

[TABLE]

having obtained $y_{0}$ from a standard, unconstrained fit. We then use this estimate to choose the size of a zero-centered chisquare prior on each fitting parameter. The procedure has almost no impact on the determined values of $y_{0}$ and $y_{1}$, where the data is far more constraining than the prior. In practice, a quadratic polynomial is sufficient to give a good fit with a reasonable $\chi^{2}$. The results of these fits are given in Table 2 and the fits themselves are displayed in Fig. 2. We also confirmed by varying the volume that any finite-volume effects are smaller than our statistical error bars.

We caution the reader that, while $y_{0}$ and $y_{1}$ can be interpreted as the continuum critical point and the 3-loop $\mathcal{O}(a)$ additive mass renormalization coefficient, we cannot interpret $y_{2}$ as a 4-loop mass renormalization or use it to further improve the lattice-continuum matching. That is because there are many uncontrolled $\mathcal{O}(a^{2})$ corrections which influence $y_{2}$. For instance, there are unknown 2-loop $\mathcal{O}(a^{2})$ lattice-continuum corrections to $x$, which influence the critical value via $(dy_{\mathrm{crit}}/dx)\delta x$. Similarly, tree-level $\mathcal{O}(a^{2})$ high-dimension operators and $\mathcal{O}(a^{2})$ corrections to $g^{2}$ (which we could interpret as uncertainties in the scale setting) also lead to $\mathcal{O}(a^{2})$ effects in the $y_{\mathrm{crit}}$ value. Because $(dy_{\mathrm{crit}}/dx)\sim x^{-2}$, the $\mathcal{O}(a^{2})$ and other higher-order effects will become severe as we go towards small $x$ values. Therefore small $x$ requires the use of very fine lattices. Furthermore, when $x$ is small, there becomes a hierarchy of mass scales in the problem; $m_{A,\mathrm{brok}}\gg m_{\Phi}\gg m_{A,\mathrm{symm}}$. Both effects make the accurate extraction of $y_{1}$ at small $x$ very numerically demanding. Therefore we did not treat the smallest $x$ value shown in Table 1. Instead, we add two larger values of $x$, $x=0.13$ and $x=0.20$, which are still within the domain where the transition is first order, but which give us a broader $x$ range over which to fit $y_{1}$ as a function of $x$.

Next we use our results for $y_{1}(x)$ to fit its overall $x$ dependence. The parametric form of the $\mathcal{O}(g_{3\mathrm{d}}^{2}a)$-correction DOnofrio:2014mld was given in (7). The $x^{3}$ coefficient corresponds to 3-loop diagrams containing only scalar lines. It is therefore equal to the $\mathcal{O}(a)$ mass renormalization term in the theory in the $g_{3\mathrm{d}}^{2}\to 0$ limit, which is a scalar theory. We explain our (different) procedure to treat this scalar theory in App. A; our analysis leads to the result

[TABLE]

Here $\tilde{y}\equiv m^{2}(\mu=\lambda)/\lambda^{2}$ is the scalar mass, made dimensionless using the scale $\lambda$ rather than the scale $g_{3\mathrm{d}}^{2}$; it equals $y/x^{2}$ up to the effect of the different renormalization scale. We incorporate this result as a prior in fitting a cubic polynomial to the results of Table 2. The resulting fit,

[TABLE]

is displayed in Fig. 3. We report the full error covariance matrix in Table 3. This fit constitutes our main result.

As a corollary, we provide an updated version of the EQCD phase diagram. The version from Kajantie:1998yc does not include continuum-extrapolated critical masses. The intercept of our EQCD fits delivers these critical masses. Additionally, the $x\rightarrow 0$ limit

[TABLE]

is known perturbatively Kajantie:1998yc . We present our data, and this limiting value, in Fig. 4. In addition, to guide the eye333In fact we expect nonanalytical behavior as $x\to 0$, due for instance to the two-loop $\Phi^{2}\ln(\Phi^{2}/\mu)$ terms in the effective potential Arnold:1992rz which give rise to $xy_{\mathrm{crit}}-3/8\pi^{2}\sim x\ln(x)$ corrections to Eq. (14)., we include a cubic fit of $xy_{\mathrm{crit}}$ as a function of $x$, displayed by a dashed line. There is quantitative agreement with the phase diagram in Kajantie:1998yc at small $x$, but at large $x$ we find that the prominent bending down of the $xy_{\mathrm{crit}}$ curve found by Kajantie:1998yc arose because they failed to take a continuum limit.

Kajantie et al Kajantie:1998yc found that the tricritical point occurs at $x=0.25$, beyond which the phase transition becomes of second order. We have not studied $x$ values larger than $x=0.2$, so we cannot make any statement about the location of the tricritical point.

As a byproduct our study also produces values for the discontinuity in $\mathrm{Tr}\,\Phi^{2}$ across the phase transition point, and for the $\mathcal{O}(a)$ additive correction to the $\mathrm{Tr}\,\Phi^{2}$ operator, which was also not previously known. We postpone these secondary results to Appendix B.

5 Conclusion and outlook

We have computed the remaining $\mathcal{O}(a)$ improvement coefficient in the lattice-continuum matching of 3D EQCD (SU(3) gauge theory with an adjoint scalar in 3 space dimensions). We did so by using the first order phase transition point as a fixed physical point. We developed a new methodology to efficiently extract the critical scalar mass where the first order transition occurs. Determining the critical scalar mass $y_{\mathrm{crit}}(x,a)$ as a function of $a$ for fixed $x$ allows an extrapolation to the continuum limit; the linear term in the extrapolation is the desired improvement coefficient. We then performed a grand fit to its known functional form. As a byproduct, a continuum-extrapolated update of the EQCD phase diagram was obtained.

Now that we are in possession of the last missing ingredient, we aim to compute the modified EQCD-Wilson-loop that leads to $C(q_{\perp})$ and extrapolate it to continuum. The continuum extrapolation is drastically facilitated by our completion of the renormalization Panero:2013pla . The jet broadening coefficient $\hat{q}$ can be derived as the second moment of $C(q_{\perp})$. We hope that the resulting nonperturbative information on $C(q_{\perp})$ will be of utility in studying and interpreting jet modification arising from the hot QCD medium in relativistic heavy ion collisions.

Acknowledgments

This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – project number 315477589 – TRR 211. Calculations for this research were conducted on the Lichtenberg high performance computer of the TU Darmstadt. We thank Kari Rummukainen and Aleksi Kurkela for useful conversations, and Daniel Robaina for his patient help with the OpenQCD-1.6 codebase.

Appendix A Algorithm for pure scalar case

The large-$x$ limit of EQCD is the same as the $g_{3\mathrm{d}}^{2}\to 0$ limit (provided we work in terms of $\tilde{y}\equiv y/x^{2}=m^{2}/\lambda^{2}$ and track the lattice spacing in terms of $\lambda a=g_{3\mathrm{d}}^{2}a\,x$). In this limit, we have an 8 (real) component scalar theory with an $\mathcal{O}(8)$ symmetry and a second-order phase transition where this symmetry is spontaneously broken to $\mathcal{O}(7)$. The $x^{3}$ term in the scalar mass renormalization of EQCD arises purely from scalar diagrams which are identical to those in this theory; therefore we can determine this coefficient by studying the $\mathcal{O}(\lambda a)$ corrections to $\tilde{y}$ in this scalar field theory.

A natural approach would be to use, as a line of constant physics, the 2’nd order phase transition point. However this would face the usual problems of critical slowing down and the inaccuracy of establishing the exact transition point. So we choose instead to compare $\tilde{y}$ between different lattice spacings by finding the pseudocritical value where the theory in a specific physical volume encounters a specific pseudocritical criterion. We choose the volume to be $\lambda L=8$ and select as the pseudocritical condition that the 4th order Binder cumulant Binder:1981sa

[TABLE]

where $\bar{\Phi}\equiv\tfrac{1}{N^{3}}\sum_{\boldsymbol{x}}\Phi_{\mathrm{L}}({\boldsymbol{x}})$, takes the value $\mathcal{B}_{\mathrm{pc}}=1.073$. Because the Binder cumulant is dominated by infrared physics and is insensitive to the lattice spacing up to subleading corrections and renormalization effects (which are what we want to study), this should occur at the same physical $y$ value at every lattice spacing up to $\mathcal{O}(a^{2})$ or higher corrections. Our set of simulation parameters can be found in Table 4. The results already appeared in the last panel of Fig. 2.

Appendix B EQCD simulation parameters and phase transition properties

B.1 Parameters and table

In Table 5, we provide parameters and direct results of our EQCD simulations as well as values for $\mathrm{Tr}\,\Phi^{2}$ in both phases at criticality. We converted $\mathrm{Tr}\,\Phi^{2}_{\mathrm{L,crit}}\rightarrow\frac{\mathrm{Tr}\,\Phi^{2}_{\mathrm{crit,cont}}}{g_{3\mathrm{d}}^{2}}$, where we display the latter in the table, via all known contributions up to $\mathcal{O}(g_{3\mathrm{d}}^{2}a)$ in Eq. (18) Moore:1997np . For the sake of readability, we will refer to $\frac{\mathrm{Tr}\,\Phi^{2}_{\mathrm{crit,cont}}}{g_{3\mathrm{d}}^{2}}$ as $\frac{\mathrm{Tr}\,\Phi^{2}}{g_{3\mathrm{d}}^{2}}$ in the following. The raw data was obtained in separate simulations with $V=N_{\mathrm{x}}^{3}$. We ensured that the Monte-Carlo error of the separate simulations dominates the overall error of $\frac{\mathrm{Tr}\,\Phi^{2}}{g_{3\mathrm{d}}^{2}}$, not the uncertainty of $y_{\mathrm{crit}}$. The slightly negative values of $\frac{\mathrm{Tr}\,\Phi^{2}_{\mathrm{symm}}}{g_{3\mathrm{d}}^{2}}$ for small $x$ are expected and arise because the positive mass-squared at the transition point suppresses IR fluctuations, while the renormalization involves subtracting off large positive counterterms (including the massless, free-theory fluctuations).

B.2 Transition strength

With the critical values of $\frac{\mathrm{Tr}\,\Phi^{2}}{g_{3\mathrm{d}}^{2}}$ in both phases from Table 5 in hand, we can infer further interesting features of the first order phase transition. Having several values of $\frac{\mathrm{Tr}\,\Phi^{2}}{g_{3\mathrm{d}}^{2}}$ at the same physical $x$ and different lattice spacings $g_{3\mathrm{d}}^{2}a$, we extrapolate both the difference between phases $\frac{\Delta\mathrm{Tr}\,\Phi^{2}}{g_{3\mathrm{d}}^{2}}$ and the symmetric phase value $\frac{\mathrm{Tr}\,\Phi^{2}_{\mathrm{symm}}}{g_{3\mathrm{d}}^{2}}$ to the continuum. We provide the former in Fig. 5 and the latter in Fig. 6. The continuum limits, and the linear coefficient in the fit for the case of $\frac{\mathrm{Tr}\,\Phi^{2}_{\mathrm{symm}}}{g_{3\mathrm{d}}^{2}}$, are provided in Table 6.

The limiting values of the fits provide us with two interesting pieces of information about the phase transition in this theory. The most interesting is $\frac{\Delta\mathrm{Tr}\,\Phi^{2}}{g_{3\mathrm{d}}^{2}}$, which measures the strength of the phase transition. For small $x$ we can predict this strength perturbatively; the limiting behavior is Kajantie:1998yc $\frac{\Delta\mathrm{Tr}\,\Phi^{2}}{g_{3\mathrm{d}}^{2}}=3/(8\pi^{2}x^{2})$, which we also include in the last frame of Fig. 5. We have provided a cubic fit to guide the eye, but it should not be taken seriously; the strength of the phase transition is a nonperturbative quantity and there is no reason to expect it to take such a simple form. In fact, we know that $\frac{\Delta\mathrm{Tr}\,\Phi^{2}}{g_{3\mathrm{d}}^{2}}\to 0$ as $x\to x_{\mathrm{triple}}$, with a nontrivial critical exponent. Note that in the fit for the $a$ dependence of $\frac{\Delta\mathrm{Tr}\,\Phi^{2}}{g_{3\mathrm{d}}^{2}}$, we have fitted to a polynomial without a linear term; this is because the known $\mathcal{O}(a)$ corrections are sufficient to eliminate such a linear correction in the difference between phases of $\frac{\mathrm{Tr}\,\Phi^{2}}{g_{3\mathrm{d}}^{2}}$.

B.3 Additive operator improvement

On the other hand, the value of either $\frac{\mathrm{Tr}\,\Phi^{2}_{\mathrm{symm}}}{g_{3\mathrm{d}}^{2}}$ or $\frac{\mathrm{Tr}\,\Phi^{2}_{\mathrm{brok}}}{g_{3\mathrm{d}}^{2}}$, by themselves, still contain $\mathcal{O}(a)$ errors, since there is an unknown additive renormalization to the operator $\mathrm{Tr}\,\Phi^{2}$ which arises at 3 loops. Since the correction is additive and both phases were explored at the same $y$ value, these effects cancel in the difference. We took advantage of this cancellation in the last subsection. But now our goal is to use this linear behavior to extract the unknown $\mathcal{O}(a)$ additive corrections to the $\Phi^{2}$ operator. These arise at 3 loops in a perturbative lattice-continuum matching calculation, which is prohibitive; so we will again try to extract them from the data.

Because we are working to 3 loops, we must specify quite carefully how 1-loop multiplicative effects will be implemented, since they can multiply one and two loop additive effects to give 3-loop level contributions, which then differ depending on our exact procedure. Here we will depart slightly from the procedure of Refs Moore:1997np ; Arnold:2001ir ; DOnofrio:2014mld . We write the continuum expectation value as

[TABLE]

where $Z_{m}$ is the 1-loop multiplicative renormalization factor of the $\mathrm{Tr}\,\Phi^{2}$ operator, and $Z_{\Phi}$ accounts for our choice of scalar field normalization on the lattice (see Eq. (3)). Examining the 3-loop diagrams, we find that certain 3-loop effects are absorbed if we define $Z_{m}$, resumming the Dyson series. That is, in Figure 7, we take the operator inserted on the 1-loop diagrams to be the resummed, rather than the bare, operator, which will sum the Dyson series, leading to an expression for $Z_{m}^{-1}$. Slightly rearranging Eq.(32,33) of Ref. Moore:1997np , we find

[TABLE]

Here $\xi=0.152859324966$ and $\Sigma=3.17591153562522$ are standard integrals encountered in the 1-loop lattice-continuum matching. We are also writing the number of colors $N=3$ explicitly, to show the detailed dependence on the number of colors. With this definition, the two-loop and partially 3-loop result for $\delta\Phi^{2}$ is

[TABLE]

The unknown coefficients $C_{\Phi a}$, $C_{\Phi b}$ and $C_{\Phi c}$ capture the remaining $\mathcal{O}(a)$ corrections from 3-loop diagrams which are not iterations of simpler 1 and 2 loop diagrams. Note that Eq. (18) contains several terms proportional to $\ln(g_{3\mathrm{d}}^{2}a)$. These arise from logarithmic divergences in the continuum theory, regulated at our choice of renormalization point $\mu=g_{3\mathrm{d}}^{2}$ but then cut off on the lattice at the scale $1/a$. The term next to $\zeta-\delta$ is the explicit $\mu^{2}$ dependence of $\Phi^{2}/g_{3\mathrm{d}}^{2}$ and is therefore expected; the coefficient $Z_{m}^{-1}$ ensures that it enters Eq. (16) with precisely the right continuum normalization. The log terms proportional to $\xi$ in the last line cancel the $\mu$ dependence of the mass squared ($y$) in the mass-dependent $\mathcal{O}(a)$ shift in the first line.

It remains to determine the three coefficients in the last line. In fact, we only need to fit two of these coefficients; one of them, $C_{\Phi c}$, represents pure-scalar corrections, which can be extracted from Ref. Arnold:2001ir . The reference performs the calculation for an improved hopping term, but repeating the calculation for the nearest-neighbor hopping term we use here555The rest of the $\mathcal{O}(a)$ corrections are not known for improved actions, which is why we do not attempt to use an improved action here. Using improved actions only really helps if one can complete the 2-loop $\mathcal{O}(a^{2})$ matching calculation; this is feasible in a scalar theory, but does not appear practical in a gauge theory as we consider here., we find that

[TABLE]

which differs from the result in the reference, $C_{4,\mathrm{Ref}}=0.2817$, because of the different scalar dispersion between the nearest-neighbor hopping term used here and the improved hopping term used there.666Specifically, in passing from the improved to the unimproved hopping term, Ref. Arnold:2001ir (B38) has $0.30837\to 0.2268854$; (B40) has $.00031757\to.000490546$, the value of $\xi$ changes from $\xi=-.08365\to+.1528593$, and (B41) changes from $.0985(6)\to.1063(4)$.

After performing these subtractions, we can fit the residual linear $a$-dependence of $\mathrm{Tr}\,\Phi^{2}_{\mathrm{L}}$ for each $x$ value we consider, and extract the coefficients $C_{\Phi a},C_{\Phi b}$ from a grand fit in complete analogy with the $m^{2}$ effects we consider in the main text. We find

[TABLE]

again by constrained curve fitting. The grand fit is displayed in the seventh frame of Fig. 6. Unfortunately, it appears that our results fail to constrain these coefficients very much. The full covariance matrix can be found in Tab. 7. From the covariance matrix, we can also see that the error of $C_{\Phi c}$ is by far the smallest, so value and error of $C_{\Phi c}$ were not changed by the constrained fit.