Numerical stochastic perturbation theory applied to the twisted Eguchi-Kawai model

TL;DR

This paper explores the application of numerical stochastic perturbation theory to the four-dimensional twisted Eguchi-Kawai model, achieving high-precision calculations of Wilson loop coefficients and analyzing their behavior as N increases.

Contribution

It demonstrates the use of NSPT with a GHMD algorithm on the TEK model, computing perturbative expansions up to O(g^8) and analyzing large N behavior.

Findings

High-precision coefficients up to O(g^4) match exact values.

Coefficients tend to Gaussian distribution as N increases.

Method can be extended to higher orders.

Abstract

We present the results of an exploratory study of the numerical stochastic perturbation theory (NSPT) applied to the four dimensional twisted Eguchi-Kawai (TEK) model. We employ a Kramers type algorithm based on the Generalized Hybrid Molecular Dynamics (GHMD) algorithm. We have computed the perturbative expansion of square Wilson loops up to $O(g^8)$. The results of the first two coefficients (up to $O(g^4)$) have a high precision and match well with the exact values. The next two coefficients can be determined and even extrapolated to large $N$, where they should coincide with the corresponding coefficients for ordinary Yang-Mills theory on an infinite lattice. Our analysis shows the behaviour of the probability distribution for each coefficient tending to Gaussian for larger $N$. The results allow us to establish the requirements to extend this analysis to much higher order.

Tables

Table 1. Table 1: Simulation parameters and statistics.

$N$	$\hat{L}$	$k$	$N_{\mathrm{MD}}$	Statistics
121	11	3	40	$\mathrm{216\hspace{0.17em}000}$
			32	$\mathrm{175\hspace{0.17em}000}$
			28	$\mathrm{182\hspace{0.17em}000}$
49	7	1	32	$\mathrm{516\hspace{0.17em}000}$
49	7	3	32	$\mathrm{520\hspace{0.17em}000}$
49	7	2	32	$\mathrm{510\hspace{0.17em}000}$
			24	$\mathrm{528\hspace{0.17em}000}$
			20	$\mathrm{504\hspace{0.17em}000}$
16	4	1	32	$\mathrm{1\hspace{0.17em}140\hspace{0.17em}000}$
			20	$\mathrm{1\hspace{0.17em}170\hspace{0.17em}000}$
			16	$\mathrm{1\hspace{0.17em}120\hspace{0.17em}000}$

Table 2. Table 2: Perturbative coefficients for Wilson loops on S 1 , S 2 subscript 𝑆 1 subscript 𝑆 2 S_{1},S_{2} planes. (SU(16), k = 1 𝑘 1 k=1 )

	$S_1$			$S_2$
	Analytic	NSPT	${\chi }^2/\mathrm{DoF}$	Analytic	NSPT	${\chi }^2/\mathrm{DoF}$
$Re{\hat{W}}_1^{11}$	$\mathrm{0.124\hspace{0.17em}511\hspace{0.17em}72}$	$\mathrm{0.124\hspace{0.17em}522\hspace{0.17em}3}$ p m 0.0000093	$1.0$	$=S_1$	$\mathrm{0.124\hspace{0.17em}512}$ p m 0.000012	$0.068$
$Re{\hat{W}}_1^{22}$	$\mathrm{0.328\hspace{0.17em}125\hspace{0.17em}00}$	$\mathrm{0.328\hspace{0.17em}138}$ p m 0.000037	$0.60$	$=S_1$	$\mathrm{0.328\hspace{0.17em}096}$ p m 0.000045	$4.4$
$Re{\hat{W}}_1^{33}$	$\mathrm{0.452\hspace{0.17em}636\hspace{0.17em}72}$	$\mathrm{0.452\hspace{0.17em}662}$ p m 0.000050	$0.35$	$=S_1$	$\mathrm{0.452\hspace{0.17em}601}$ p m 0.000057	$1.8$
$Re{\hat{W}}_2^{11}$	$\mathrm{0.005\hspace{0.17em}055\hspace{0.17em}99}$	$\mathrm{0.005\hspace{0.17em}047\hspace{0.17em}3}$ p m 0.0000096	$0.057$	$\mathrm{0.005\hspace{0.17em}042\hspace{0.17em}42}$	$\mathrm{0.005\hspace{0.17em}055}$ p m 0.000013	$0.52$
$Re{\hat{W}}_2^{22}$	$-\mathrm{0.014\hspace{0.17em}673\hspace{0.17em}13}$	$-\mathrm{0.014\hspace{0.17em}691}$ p m 0.000043	$1.0$	$-\mathrm{0.014\hspace{0.17em}736\hspace{0.17em}99}$	$-\mathrm{0.014\hspace{0.17em}698}$ p m 0.000056	$4.1$
$Re{\hat{W}}_2^{33}$	$-\mathrm{0.045\hspace{0.17em}993\hspace{0.17em}08}$	$-\mathrm{0.045\hspace{0.17em}988}$ p m 0.000062	$2.3$	$-\mathrm{0.045\hspace{0.17em}920\hspace{0.17em}18}$	$-\mathrm{0.045\hspace{0.17em}817}$ p m 0.000074	$2.9$
$Re{\hat{W}}_2^{44}$	$\mathrm{0.158\hspace{0.17em}384\hspace{0.17em}33}$	$\mathrm{0.158\hspace{0.17em}400}$ p m 0.000045	9.4 $\times {10}^{-8}$	$\mathrm{0.157\hspace{0.17em}896\hspace{0.17em}05}$	$\mathrm{0.157\hspace{0.17em}894}$ p m 0.000045	$0.023$
$Re{\hat{W}}_3^{11}$		$\mathrm{0.000\hspace{0.17em}798}$ p m 0.000014	$1.4$		$\mathrm{0.000\hspace{0.17em}775}$ p m 0.000020	$1.8$
$Re{\hat{W}}_3^{22}$		$\mathrm{0.000\hspace{0.17em}114}$ p m 0.000069	$0.38$		$\mathrm{0.000\hspace{0.17em}051}$ p m 0.000089	$2.3$
$Re{\hat{W}}_3^{33}$		$\mathrm{0.000\hspace{0.17em}481}$ p m 0.000096	$0.010$		$\mathrm{0.000\hspace{0.17em}38}$ p m 0.00012	$4.4$
$Re{\hat{W}}_3^{44}$		$-\mathrm{0.023\hspace{0.17em}093}$ p m 0.000053	$0.050$		$-\mathrm{0.022\hspace{0.17em}910}$ p m 0.000056	$1.2$
$Re{\hat{W}}_4^{11}$		$\mathrm{0.000\hspace{0.17em}118}$ p m 0.000025	$2.7$		$\mathrm{0.000\hspace{0.17em}179}$ p m 0.000033	$1.1$
$Re{\hat{W}}_4^{22}$		$-\mathrm{0.000\hspace{0.17em}13}$ p m 0.00012	$2.5$		$\mathrm{0.000\hspace{0.17em}10}$ p m 0.00015	$1.2$
$Re{\hat{W}}_4^{33}$		$-\mathrm{0.000\hspace{0.17em}34}$ p m 0.00016	$0.51$		$-\mathrm{0.000\hspace{0.17em}18}$ p m 0.00020	$3.2$
$Re{\hat{W}}_4^{44}$		$\mathrm{0.000\hspace{0.17em}837}$ p m 0.000083	$0.72$		$\mathrm{0.000\hspace{0.17em}775}$ p m 0.000093	$0.83$
$\left|Im{\hat{W}}_2^{11}\right|$	$\mathrm{0.000\hspace{0.17em}008\hspace{0.17em}51}$	$\mathrm{0.000\hspace{0.17em}009\hspace{0.17em}62}$ p m 0.00000100	$0.37$	$\mathrm{0.000\hspace{0.17em}016\hspace{0.17em}08}$	$\mathrm{0.000\hspace{0.17em}015\hspace{0.17em}2}$ p m 0.0000014	$2.3$
$\left|Im{\hat{W}}_2^{22}\right|$	$\mathrm{0.000\hspace{0.17em}208\hspace{0.17em}33}$	$\mathrm{0.000\hspace{0.17em}208\hspace{0.17em}4}$ p m 0.0000036	$0.90$	$\mathrm{0.000\hspace{0.17em}044\hspace{0.17em}85}$	$\mathrm{0.000\hspace{0.17em}050\hspace{0.17em}7}$ p m 0.0000051	$1.4$
$\left|Im{\hat{W}}_2^{33}\right|$	$\mathrm{0.001\hspace{0.17em}358\hspace{0.17em}65}$	$\mathrm{0.001\hspace{0.17em}366\hspace{0.17em}6}$ p m 0.0000074	$1.1$	$\mathrm{0.001\hspace{0.17em}056\hspace{0.17em}71}$	$\mathrm{0.001\hspace{0.17em}060}$ p m 0.000010	$\mathrm{0.000\hspace{0.17em}16}$
$\left|Im{\hat{W}}_3^{11}\right|$		$\mathrm{0.000\hspace{0.17em}000\hspace{0.17em}4}$ p m 0.0000021	$0.23$		$\mathrm{0.000\hspace{0.17em}004\hspace{0.17em}6}$ p m 0.0000031	$0.13$
$\left|Im{\hat{W}}_3^{22}\right|$		$\mathrm{0.000\hspace{0.17em}002\hspace{0.17em}3}$ p m 0.0000071	$1.7$		$\mathrm{0.000\hspace{0.17em}016}$ p m 0.000010	$0.075$
$\left|Im{\hat{W}}_3^{33}\right|$		$\mathrm{0.000\hspace{0.17em}086}$ p m 0.000015	$0.99$		$\mathrm{0.000\hspace{0.17em}024}$ p m 0.000022	$0.011$
$\left|Im{\hat{W}}_3^{44}\right|$		$\mathrm{0.000\hspace{0.17em}000\hspace{0.17em}7}$ p m 0.0000012	$0.18$		$\mathrm{0.000\hspace{0.17em}001\hspace{0.17em}2}$ p m 0.0000017	$0.035$
$\left|Im{\hat{W}}_4^{11}\right|$		$\mathrm{0.000\hspace{0.17em}000\hspace{0.17em}7}$ p m 0.0000042	$1.5$		$\mathrm{0.000\hspace{0.17em}003\hspace{0.17em}0}$ p m 0.0000059	$0.94$
$\left|Im{\hat{W}}_4^{22}\right|$		$\mathrm{0.000\hspace{0.17em}004}$ p m 0.000013	$0.24$		$\mathrm{0.000\hspace{0.17em}017}$ p m 0.000019	$\mathrm{0.000\hspace{0.17em}84}$
$\left|Im{\hat{W}}_4^{33}\right|$		$\mathrm{0.000\hspace{0.17em}028}$ p m 0.000031	$0.22$		$\mathrm{0.000\hspace{0.17em}018}$ p m 0.000045	$0.65$
$\left|Im{\hat{W}}_4^{44}\right|$		$\mathrm{0.000\hspace{0.17em}005\hspace{0.17em}2}$ p m 0.0000036	$0.20$		$\mathrm{0.000\hspace{0.17em}002\hspace{0.17em}8}$ p m 0.0000051	$2.6$

Table 3. Table 3: Perturbative coefficients for Wilson loops on S 1 , S 2 subscript 𝑆 1 subscript 𝑆 2 S_{1},S_{2} planes. (SU(49), k = 2 𝑘 2 k=2 )

	$S_1$			$S_2$
	Analytic	NSPT	${\chi }^2/\mathrm{DoF}$	Analytic	NSPT	${\chi }^2/\mathrm{DoF}$
$Re{\hat{W}}_1^{11}$	$\mathrm{0.124\hspace{0.17em}947\hspace{0.17em}94}$	$\mathrm{0.124\hspace{0.17em}954\hspace{0.17em}2}$ p m 0.0000053	$0.030$	$=S_1$	$\mathrm{0.124\hspace{0.17em}957\hspace{0.17em}9}$ p m 0.0000065	$0.038$
$Re{\hat{W}}_1^{22}$	$\mathrm{0.341\hspace{0.17em}345\hspace{0.17em}68}$	$\mathrm{0.341\hspace{0.17em}372}$ p m 0.000021	$2.7$	$=S_1$	$\mathrm{0.341\hspace{0.17em}395}$ p m 0.000025	$0.78$
$Re{\hat{W}}_1^{33}$	$\mathrm{0.570\hspace{0.17em}101\hspace{0.17em}77}$	$\mathrm{0.570\hspace{0.17em}171}$ p m 0.000046	$0.34$	$=S_1$	$\mathrm{0.570\hspace{0.17em}128}$ p m 0.000054	$1.6$
$Re{\hat{W}}_1^{44}$	$\mathrm{0.789\hspace{0.17em}492\hspace{0.17em}99}$	$\mathrm{0.789\hspace{0.17em}596}$ p m 0.000071	$0.44$	$=S_1$	$\mathrm{0.789\hspace{0.17em}543}$ p m 0.000082	$1.1$
$Re{\hat{W}}_2^{11}$	$\mathrm{0.005\hspace{0.17em}101\hspace{0.17em}03}$	$\mathrm{0.005\hspace{0.17em}103\hspace{0.17em}7}$ p m 0.0000061	$7.6$	$\mathrm{0.005\hspace{0.17em}100\hspace{0.17em}26}$	$\mathrm{0.005\hspace{0.17em}101\hspace{0.17em}2}$ p m 0.0000077	$0.22$
$Re{\hat{W}}_2^{22}$	$-\mathrm{0.016\hspace{0.17em}678\hspace{0.17em}06}$	$-\mathrm{0.016\hspace{0.17em}682}$ p m 0.000026	$0.013$	$-\mathrm{0.016\hspace{0.17em}680\hspace{0.17em}76}$	$-\mathrm{0.016\hspace{0.17em}677}$ p m 0.000033	$1.0$
$Re{\hat{W}}_2^{33}$	$-\mathrm{0.088\hspace{0.17em}228\hspace{0.17em}03}$	$-\mathrm{0.088\hspace{0.17em}249}$ p m 0.000061	$0.0043$	$-\mathrm{0.088\hspace{0.17em}235\hspace{0.17em}09}$	$-\mathrm{0.088\hspace{0.17em}183}$ p m 0.000074	$0.018$
$Re{\hat{W}}_2^{44}$	$-\mathrm{0.206\hspace{0.17em}182\hspace{0.17em}50}$	$-\mathrm{0.206\hspace{0.17em}18}$ p m 0.00010	$0.12$	$-\mathrm{0.206\hspace{0.17em}240\hspace{0.17em}12}$	$-\mathrm{0.206\hspace{0.17em}13}$ p m 0.00012	$0.034$
$Re{\hat{W}}_3^{11}$		$\mathrm{0.000\hspace{0.17em}774\hspace{0.17em}1}$ p m 0.0000096	$1.1$		$\mathrm{0.000\hspace{0.17em}784}$ p m 0.000013	$1.3$
$Re{\hat{W}}_3^{22}$		$-\mathrm{0.000\hspace{0.17em}029}$ p m 0.000043	$0.021$		$-\mathrm{0.000\hspace{0.17em}107}$ p m 0.000054	$2.6$
$Re{\hat{W}}_3^{33}$		$\mathrm{0.002\hspace{0.17em}795}$ p m 0.000095	$0.068$		$\mathrm{0.002\hspace{0.17em}64}$ p m 0.00012	$0.064$
$Re{\hat{W}}_3^{44}$		$\mathrm{0.019\hspace{0.17em}48}$ p m 0.00016	$0.0025$		$\mathrm{0.019\hspace{0.17em}24}$ p m 0.00019	$0.0023$
$Re{\hat{W}}_4^{11}$		$\mathrm{0.000\hspace{0.17em}177}$ p m 0.000018	$0.48$		$\mathrm{0.000\hspace{0.17em}171}$ p m 0.000024	$0.64$
$Re{\hat{W}}_4^{22}$		$\mathrm{0.000\hspace{0.17em}128}$ p m 0.000077	$0.044$		$\mathrm{0.000\hspace{0.17em}254}$ p m 0.000097	$0.22$
$Re{\hat{W}}_4^{33}$		$\mathrm{0.000\hspace{0.17em}21}$ p m 0.00017	$0.019$		$\mathrm{0.000\hspace{0.17em}55}$ p m 0.00021	$0.43$
$Re{\hat{W}}_4^{44}$		$-\mathrm{0.000\hspace{0.17em}17}$ p m 0.00027	$0.045$		$\mathrm{0.000\hspace{0.17em}38}$ p m 0.00034	$0.28$
$\left|Im{\hat{W}}_2^{11}\right|$	$\mathrm{0.000\hspace{0.17em}000\hspace{0.17em}29}$	$\mathrm{0.000\hspace{0.17em}000\hspace{0.17em}55}$ p m 0.00000059	$0.31$	$\mathrm{0.000\hspace{0.17em}001\hspace{0.17em}02}$	$\mathrm{0.000\hspace{0.17em}000\hspace{0.17em}01}$ p m 0.00000085	$1.7$
$\left|Im{\hat{W}}_2^{22}\right|$	$\mathrm{0.000\hspace{0.17em}014\hspace{0.17em}24}$	$\mathrm{0.000\hspace{0.17em}013\hspace{0.17em}5}$ p m 0.0000027	$0.32$	$\mathrm{0.000\hspace{0.17em}001\hspace{0.17em}99}$	$\mathrm{0.000\hspace{0.17em}004\hspace{0.17em}1}$ p m 0.0000038	$0.35$
$\left|Im{\hat{W}}_2^{33}\right|$	$\mathrm{0.000\hspace{0.17em}018\hspace{0.17em}24}$	$\mathrm{0.000\hspace{0.17em}016\hspace{0.17em}4}$ p m 0.0000063	$1.3$	$\mathrm{0.000\hspace{0.17em}064\hspace{0.17em}97}$	$\mathrm{0.000\hspace{0.17em}054\hspace{0.17em}4}$ p m 0.0000089	$0.097$
$\left|Im{\hat{W}}_2^{44}\right|$	$\mathrm{0.000\hspace{0.17em}383\hspace{0.17em}14}$	$\mathrm{0.000\hspace{0.17em}386}$ p m 0.000011	$1.2$	$\mathrm{0.000\hspace{0.17em}447\hspace{0.17em}56}$	$\mathrm{0.000\hspace{0.17em}438}$ p m 0.000017	$0.14$
$\left|Im{\hat{W}}_3^{11}\right|$		$\mathrm{0.000\hspace{0.17em}000\hspace{0.17em}1}$ p m 0.0000013	$0.0017$		$\mathrm{0.000\hspace{0.17em}002\hspace{0.17em}2}$ p m 0.0000020	$1.5$
$\left|Im{\hat{W}}_3^{22}\right|$		$\mathrm{0.000\hspace{0.17em}005\hspace{0.17em}9}$ p m 0.0000056	$3.9$		$\mathrm{0.000\hspace{0.17em}004\hspace{0.17em}4}$ p m 0.0000076	$0.75$
$\left|Im{\hat{W}}_3^{33}\right|$		$\mathrm{0.000\hspace{0.17em}018}$ p m 0.000014	$1.4$		$\mathrm{0.000\hspace{0.17em}031}$ p m 0.000018	$0.67$
$\left|Im{\hat{W}}_3^{44}\right|$		$\mathrm{0.000\hspace{0.17em}102}$ p m 0.000025	$0.61$		$\mathrm{0.000\hspace{0.17em}088}$ p m 0.000037	$0.053$
$\left|Im{\hat{W}}_4^{11}\right|$		$\mathrm{0.000\hspace{0.17em}000\hspace{0.17em}8}$ p m 0.0000026	$0.044$		$\mathrm{0.000\hspace{0.17em}003\hspace{0.17em}5}$ p m 0.0000039	$0.22$
$\left|Im{\hat{W}}_4^{22}\right|$		$\mathrm{0.000\hspace{0.17em}003}$ p m 0.000011	$2.2$		$\mathrm{0.000\hspace{0.17em}008}$ p m 0.000015	$0.39$
$\left|Im{\hat{W}}_4^{33}\right|$		$\mathrm{0.000\hspace{0.17em}056}$ p m 0.000029	$0.25$		$\mathrm{0.000\hspace{0.17em}068}$ p m 0.000039	$1.8$
$\left|Im{\hat{W}}_4^{44}\right|$		$\mathrm{0.000\hspace{0.17em}108}$ p m 0.000053	$0.86$		$\mathrm{0.000\hspace{0.17em}021}$ p m 0.000079	$0.15$

Table 4. Table 4: Perturbative coefficients for Wilson loops ( 1 / N MD → 0 → 1 subscript 𝑁 MD 0 1/N_{\mathrm{MD}}\to 0 extrapolated). Data in PT Analytic are obtained analytically and data in NPT Fit are obtained by fitting the non-perturbative data with a polynomial of λ 𝜆 \lambda   Perez:2017jyq .

SU(16), $k=1$	$\ell =1$	$\ell =2$	$\ell =3$	$\ell =4$
${\hat{W}}_{\ell }^{11}$	$\mathrm{0.124\hspace{0.17em}518\hspace{0.17em}8}$ p m 0.0000086	$\mathrm{0.005\hspace{0.17em}049\hspace{0.17em}7}$ p m 0.0000084	$\mathrm{0.000\hspace{0.17em}791}$ p m 0.000012	$\mathrm{0.000\hspace{0.17em}138}$ p m 0.000021
${\hat{W}}_{\ell }^{22}$	$\mathrm{0.328\hspace{0.17em}123}$ p m 0.000034	$-\mathrm{0.014\hspace{0.17em}695}$ p m 0.000038	$\mathrm{0.000\hspace{0.17em}094}$ p m 0.000061	$-\mathrm{0.000\hspace{0.17em}05}$ p m 0.00011
${\hat{W}}_{\ell }^{33}$	$\mathrm{0.452\hspace{0.17em}642}$ p m 0.000047	$-\mathrm{0.045\hspace{0.17em}930}$ p m 0.000056	$\mathrm{0.000\hspace{0.17em}448}$ p m 0.000088	$-\mathrm{0.000\hspace{0.17em}29}$ p m 0.00015
${\hat{W}}_{\ell }^{44}$	$0.0$	$\mathrm{0.158\hspace{0.17em}232}$ p m 0.000044	$-\mathrm{0.023\hspace{0.17em}032}$ p m 0.000051	$\mathrm{0.000\hspace{0.17em}816}$ p m 0.000079
	PT Analytic	PT Analytic	NPT Fit
${\hat{W}}_{\ell }^{11}$	$\mathrm{0.124\hspace{0.17em}511\hspace{0.17em}72}$	$\mathrm{0.005\hspace{0.17em}051\hspace{0.17em}46}$	$\mathrm{0.000\hspace{0.17em}826}\left(36\right)$
${\hat{W}}_{\ell }^{22}$	$\mathrm{0.328\hspace{0.17em}125\hspace{0.17em}00}$	$-\mathrm{0.014\hspace{0.17em}694\hspace{0.17em}41}$	$\mathrm{0.000\hspace{0.17em}04}\left(12\right)$
${\hat{W}}_{\ell }^{33}$	$\mathrm{0.452\hspace{0.17em}636\hspace{0.17em}72}$	$-\mathrm{0.045\hspace{0.17em}968\hspace{0.17em}78}$	$\mathrm{0.000\hspace{0.17em}34}\left(16\right)$
${\hat{W}}_{\ell }^{44}$	$0.0$	$\mathrm{0.158\hspace{0.17em}221\hspace{0.17em}57}$	$-\mathrm{0.022\hspace{0.17em}88}\left(6\right)$
SU(49), $k=2$	$\ell =1$	$\ell =2$	$\ell =3$	$\ell =4$
${\hat{W}}_{\ell }^{11}$	$\mathrm{0.124\hspace{0.17em}955\hspace{0.17em}5}$ p m 0.0000047	$\mathrm{0.005\hspace{0.17em}103\hspace{0.17em}0}$ p m 0.0000053	$\mathrm{0.000\hspace{0.17em}777\hspace{0.17em}4}$ p m 0.0000085	$\mathrm{0.000\hspace{0.17em}175}$ p m 0.000016
${\hat{W}}_{\ell }^{22}$	$\mathrm{0.341\hspace{0.17em}380}$ p m 0.000019	$-\mathrm{0.016\hspace{0.17em}680}$ p m 0.000024	$-\mathrm{0.000\hspace{0.17em}055}$ p m 0.000040	$\mathrm{0.000\hspace{0.17em}171}$ p m 0.000071
${\hat{W}}_{\ell }^{33}$	$\mathrm{0.570\hspace{0.17em}156}$ p m 0.000042	$-\mathrm{0.088\hspace{0.17em}227}$ p m 0.000055	$\mathrm{0.002\hspace{0.17em}743}$ p m 0.000084	$\mathrm{0.000\hspace{0.17em}33}$ p m 0.00015
${\hat{W}}_{\ell }^{44}$	$\mathrm{0.789\hspace{0.17em}578}$ p m 0.000066	$-\mathrm{0.206\hspace{0.17em}163}$ p m 0.000097	$\mathrm{0.019\hspace{0.17em}40}$ p m 0.00014	$\mathrm{0.000\hspace{0.17em}01}$ p m 0.00025
	PT Analytic	PT Analytic	NPT Fit
${\hat{W}}_{\ell }^{11}$	$\mathrm{0.124\hspace{0.17em}947\hspace{0.17em}94}$	$\mathrm{0.005\hspace{0.17em}100\hspace{0.17em}77}$	$\mathrm{0.000\hspace{0.17em}883}\left(9\right)$
${\hat{W}}_{\ell }^{22}$	$\mathrm{0.341\hspace{0.17em}345\hspace{0.17em}68}$	$-\mathrm{0.016\hspace{0.17em}678\hspace{0.17em}96}$	$\mathrm{0.000\hspace{0.17em}086}\left(36\right)$
${\hat{W}}_{\ell }^{33}$	$\mathrm{0.570\hspace{0.17em}101\hspace{0.17em}77}$	$-\mathrm{0.088\hspace{0.17em}230\hspace{0.17em}39}$	$\mathrm{0.002\hspace{0.17em}93}\left(5\right)$
${\hat{W}}_{\ell }^{44}$	$\mathrm{0.789\hspace{0.17em}492\hspace{0.17em}99}$	$-\mathrm{0.206\hspace{0.17em}201\hspace{0.17em}71}$	$\mathrm{0.019\hspace{0.17em}59}\left(7\right)$
SU(121), $k=3$	$\ell =1$	$\ell =2$	$\ell =3$	$\ell =4$
${\hat{W}}_{\ell }^{11}$	$\mathrm{0.124\hspace{0.17em}988\hspace{0.17em}7}$ p m 0.0000037	$\mathrm{0.005\hspace{0.17em}103\hspace{0.17em}6}$ p m 0.0000041	$\mathrm{0.000\hspace{0.17em}788\hspace{0.17em}4}$ p m 0.0000069	$\mathrm{0.000\hspace{0.17em}156}$ p m 0.000012
${\hat{W}}_{\ell }^{22}$	$\mathrm{0.342\hspace{0.17em}179}$ p m 0.000015	$-\mathrm{0.016\hspace{0.17em}760}$ p m 0.000019	$-\mathrm{0.000\hspace{0.17em}030}$ p m 0.000032	$\mathrm{0.000\hspace{0.17em}020}$ p m 0.000055
${\hat{W}}_{\ell }^{33}$	$\mathrm{0.575\hspace{0.17em}489}$ p m 0.000034	$-\mathrm{0.090\hspace{0.17em}700}$ p m 0.000043	$\mathrm{0.002\hspace{0.17em}995}$ p m 0.000068	$\mathrm{0.000\hspace{0.17em}03}$ p m 0.00012
${\hat{W}}_{\ell }^{44}$	$\mathrm{0.812\hspace{0.17em}537}$ p m 0.000056	$-\mathrm{0.220\hspace{0.17em}758}$ p m 0.000080	$\mathrm{0.022\hspace{0.17em}06}$ p m 0.00012	$-\mathrm{0.000\hspace{0.17em}39}$ p m 0.00021
${\hat{W}}_{\ell }^{55}$	$\mathrm{1.049\hspace{0.17em}316}$ p m 0.000088	$-\mathrm{0.407\hspace{0.17em}28}$ p m 0.00014	$\mathrm{0.070\hspace{0.17em}62}$ p m 0.00021	$-\mathrm{0.004\hspace{0.17em}59}$ p m 0.00038
${\hat{W}}_{\ell }^{66}$	$\mathrm{1.281\hspace{0.17em}72}$ p m 0.00012	$-\mathrm{0.645\hspace{0.17em}07}$ p m 0.00021	$\mathrm{0.160\hspace{0.17em}12}$ p m 0.00033	$-\mathrm{0.019\hspace{0.17em}33}$ p m 0.00057
	PT Analytic	PT Analytic privcomm
${\hat{W}}_{\ell }^{11}$	$\mathrm{0.124\hspace{0.17em}991\hspace{0.17em}46}$	$\mathrm{0.005\hspace{0.17em}105\hspace{0.17em}92}$
${\hat{W}}_{\ell }^{22}$	$\mathrm{0.342\hspace{0.17em}183\hspace{0.17em}32}$	$-\mathrm{0.016\hspace{0.17em}797\hspace{0.17em}72}$
${\hat{W}}_{\ell }^{33}$	$\mathrm{0.575\hspace{0.17em}499\hspace{0.17em}81}$	$-\mathrm{0.090\hspace{0.17em}734\hspace{0.17em}26}$
${\hat{W}}_{\ell }^{44}$	$\mathrm{0.812\hspace{0.17em}550\hspace{0.17em}00}$	$-\mathrm{0.220\hspace{0.17em}851\hspace{0.17em}02}$
${\hat{W}}_{\ell }^{55}$	$\mathrm{1.049\hspace{0.17em}336\hspace{0.17em}25}$	$-\mathrm{0.407\hspace{0.17em}409\hspace{0.17em}49}$
${\hat{W}}_{\ell }^{66}$	$\mathrm{1.281\hspace{0.17em}751\hspace{0.17em}42}$	$-\mathrm{0.645\hspace{0.17em}311\hspace{0.17em}88}$

Table 5. Table 5: Perturbative coefficients for W ^ ℓ R ​ R subscript superscript ^ 𝑊 𝑅 𝑅 ℓ \hat{W}^{RR}_{\ell} in N → ∞ → 𝑁 N\to\infty . Values in PT Analytic for W ^ ℓ 11 subscript superscript ^ 𝑊 11 ℓ \hat{W}^{11}_{\ell} are N → ∞ → 𝑁 N\to\infty limit of results from lattice SU ​ ( N ) SU 𝑁 \mathrm{SU}(N) gauge theory with Wilson gauge action  Alles:1998is .

$N\to \infty$	$\ell =1$	$\ell =2$	$\ell =3$
${\hat{W}}_{\ell }^{11}$	$\mathrm{0.124\hspace{0.17em}999\hspace{0.17em}9}$ p m 0.0000054	$\mathrm{0.005\hspace{0.17em}113\hspace{0.17em}4}$ p m 0.0000059	$\mathrm{0.000\hspace{0.17em}790\hspace{0.17em}5}$ p m 0.0000099
${\hat{W}}_{\ell }^{22}$	$\mathrm{0.342\hspace{0.17em}336}$ p m 0.000021	$-\mathrm{0.016\hspace{0.17em}775}$ p m 0.000027	$-\mathrm{0.000\hspace{0.17em}025}$ p m 0.000046
${\hat{W}}_{\ell }^{33}$	$\mathrm{0.576\hspace{0.17em}534}$ p m 0.000049	$-\mathrm{0.091\hspace{0.17em}185}$ p m 0.000063	$\mathrm{0.003\hspace{0.17em}044}$ p m 0.000098
${\hat{W}}_{\ell }^{44}$	$\mathrm{0.817\hspace{0.17em}004}$ p m 0.000080	$-\mathrm{0.223\hspace{0.17em}62}$ p m 0.00011	$\mathrm{0.022\hspace{0.17em}59}$ p m 0.00017
	PT Analytic	PT Analytic	PT Analytic
${\hat{W}}_{\ell }^{11}$	$\mathrm{0.125\hspace{0.17em}000\hspace{0.17em}00}$	$\mathrm{0.005\hspace{0.17em}106\hspace{0.17em}929\hspace{0.17em}7}$	$\mathrm{0.000\hspace{0.17em}794\hspace{0.17em}223}$ p m 0.000000019
${\hat{W}}_{\ell }^{22}$	$\mathrm{0.342\hspace{0.17em}327\hspace{0.17em}88}$	$-\mathrm{0.016\hspace{0.17em}814}$
${\hat{W}}_{\ell }^{33}$	$\mathrm{0.576\hspace{0.17em}298\hspace{0.17em}27}$	$-\mathrm{0.091\hspace{0.17em}07}$
${\hat{W}}_{\ell }^{44}$	$\mathrm{0.815\hspace{0.17em}370\hspace{0.17em}98}$	$-\mathrm{0.222\hspace{0.17em}87}$

Equations

$$Z$$

$$S[U]$$

$$z_{\mu \nu }$$

$$n_{\mu \nu }$$

$${ϵ}_{\mu \nu }$$

$$W_{\mu \nu }(R,T)$$

$$U_{\mu }=e^{-ig{A}_{\mu }}{\Gamma }_{\mu },$$

$${\Gamma }_{\mu }{\Gamma }_{\nu }={\Gamma }_{\nu }{\Gamma }_{\mu }{z}_{\mu \nu }^{\ast }.$$

$$F[G]$$

$$U_{\mu }\to U_{\mu }^{\prime }=G{U}_{\mu }{G}^{†},$$

$$G$$

$$\Theta$$

$$Y$$

$$G(t)$$

$$\frac{dw}{dt}$$

$$U_{\mu }(t)$$

$$Z_{\mathrm{HMD}}$$

$$H[P,U]$$

$$P_{0,\mu }$$

$$\mathrm{Prob}({\eta }_{\mu }^a)d{\eta }^a$$

$${\dot{U}}_{\mu }$$

$${\dot{P}}_{\mu }$$

$$F_{\mu }$$

$$S_{\mu }$$

$$P_{0,\mu }=c_1{P}_{\mu }+\sqrt{1-c_1^2}{\eta }_{\mu },$$

$$\{\begin{array}{rl}{P}_{\mu }& =c_1{P}_{0,\mu }+\sqrt{1-c_1^2}{\eta }_{\mu },\\ w& =c_1{w}_0,\\ U_{\mu }& =U_{0,\mu },\end{array}$$

$$\{\begin{array}{rl}{\bar{P}}_{\mu }& =P_{\mu }+(\Delta t/2)F_{\mu }[U],\\ U_{\mu }^{\prime }& =\exp \left[i{\bar{P}}_{\mu }\Delta t\right]U_{\mu },\\ P_{\mu }^{\prime }& ={\bar{P}}_{\mu }+(\Delta t/2)F_{\mu }[U^{\prime }],\end{array}$$

$$\{\begin{array}{rl}{w}_1& =w+\Delta t\alpha \Theta [U^{\prime }],\\ P_{1,\mu }& =\exp \left[i{w}_1\Delta t\right]P_{\mu }^{\prime }\exp \left[-i{w}_1\Delta t\right],\\ U_{1,\mu }& =\exp \left[i{w}_1\Delta t\right]U_{\mu }^{\prime }\exp \left[-i{w}_1\Delta t\right],\end{array}$$

$${\dot{U}}_{\mu }$$

$${\dot{P}}_{\mu }$$

$$\dot{w}$$

$$D_{\mu }w$$

$$⟨{\zeta }_{\mu }^a(t){\zeta }_{\nu }^b(s)⟩$$

$$C$$

$$A=\sum _{k=0}^{\infty }{g}^k{A}^{(k)},B$$

$$C^{(k)}=(A\star B{)}^{(k)}$$

$$P_{\mu }$$

$$U_{\mu }$$

$$w$$

$$\{\begin{array}{rl}{\bar{P}}_{\mu }^{(k)}& =P_{\mu }^{(k)}+(\Delta t^{\prime }/2)F_{\mu }^{(k)}[U],\\ {U_{\mu }^{(k)}}^{\prime }& ={\left(\exp \left[i{\bar{P}}_{\mu }\Delta t^{\prime }\right]\star U_{\mu }\right)}^{(k)},\\ {P_{\mu }^{(k)}}^{\prime }& ={\bar{P}}_{\mu }^{(k)}+(\Delta t^{\prime }/2)F_{\mu }^{(k)}[U^{\prime }],\end{array}$$

Full text

††preprint:

IFT-UAM/CSIC-19-18

FTUAM-19-4

HUPD-1903

1 Introduction

Gauge theories in the limit of infinite number of colours ($N\longrightarrow\infty$) are very interesting theoretically tHooft:1973alw . They are simpler than their finite $N$ counterparts, but share most of the fascinating properties of the latter. Nonetheless, their understanding remains a challenge. The relevance of this goal is also given by the fact that they seem a point of contact with other approaches Maldacena:1997re . A good deal of their difficulty lies in the non-perturbative character of most of its properties. The standard first principles approach to this kind of problems is the lattice formulation of quantum field theory Wilson:1974sk . However, in contrast with what happens in perturbation theory, the large $N$ limit involves an extrapolation and seems harder than the finite $N$ study. Only recently this line of approach has led to trustworthy computations (for a review see  Lucini:2012gg ), pioneered by the works of Teper and collaborators Lucini:2003zr ; Lucini:2004my ; Lucini:2005vg .

Fortunately, there is a very specific simplification which emerges when studying large $N$ gauge theories on the lattice. This is the so-called Eguchi-Kawai (EK) reduction Eguchi:1982nm . According to this result finite volume corrections are subleading in the large $N$ limit. Hence, there is the possibility that the dynamics of the large $N$ gauge theories is captured by a matrix model. There are several proposals that have been put forward to transform this possibility into a reality Bhanot:1982sh ; Gross:1982at ; Narayanan:2003fc ; Kovtun:2007py ; Unsal:2008ch . Here, we will be concerned with one of the early proposals called Twisted Eguchi-Kawai model (TEK for short) GonzalezArroyo:1982ub ; GonzalezArroyo:1982hz , and introduced by two of the present authors. Indeed, it was used shortly after EK proposal to compute the string tension at large $N$ GonzalezArroyo:1983pw ; Fabricius:1984un . Lately this has led to a calculation of this quantity with at least compatible precision to other methods GonzalezArroyo:2012fx . Several precise tests of the reduction mechanism have been obtained recently for other quantities Gonzalez-Arroyo:2014dua , providing a strong verification of the validity and usefulness of this approach.

The previous paragraph justifies our interest in the TEK model. The finite $N$ corrections are different for the matrix model and for the ordinary gauge theory, and their size and nature are very important from the practical point of view, in order to make this approach competitive computationally. Curiously, these corrections have a theoretically interesting interpretation in terms of the so-called non-commutative gauge theories Douglas:2001ba . Their Lagrangian and Feynman rules first appeared in the literature GonzalezArroyo:1983ac when looking for a continuous generalization of the TEK model, and before they emerged from the mathematical construction of non-commutative geometries Connes:1987ue . A more direct connection appears as a result of Morita duality on the non-commutative torus. This shows that ordinary gauge theories with twisted boundary conditions (TBC) a la ’t Hooft tHooft:1979rtg are particular cases of non-commutative field theories. The TEK model is nothing but the volume reduced version of a gauge theory with twisted boundary conditions on the lattice. TBC are characterized by a collection of integer-valued fluxes, and their appropriate choice has been found to have a crucial impact upon the size of the corrections and the absence of phase transitions.

The last ingredient entering in this work is perturbation theory. Although our major interest when using the lattice approach was in studying the non-perturbative aspects of the theory, there is also interest in understanding the theory from the perturbative side. At the least it offers us a method to analytically determine certain observables and estimating the size and nature of its $N$ dependence. In this spirit a recent perturbative calculation of Wilson loops in lattice gauge theories with twisted boundary conditions (including the TEK model) was addressed Perez:2017jyq .

In addition, a new set of ideas has focused in explaining the usefulness of perturbation theory in understanding also non-perturbative contributions. This takes its most extreme form in the concept of resurgence (see for example Dunne:2015eaa ; Cherman:2014ofa and references therein). At the least, as was already known, the large order behaviour of the perturbative coefficients is associated with non-perturbative aspects such as the action of other saddles. Very interesting results have been obtained in this spirit for large $N$ matrix theories Cherman:2014ofa ; Marino:2012zq . For example, it has been identified that the large $N$ phase transition of the Gross-Witten-Wadia unitary matrix model Gross:1980he ; Wadia:2012fr is governed by non-trivial saddles and the action can be reconstructed via the resurgence on the asymptotic expansion about the vacuum Buividovich:2015oju .

For large $N$ gauge theories in 4 dimensions there are many interesting aspects to be studied. Perturbation theory is dominated by planar diagrams, whose number does not grow factorially. This could suggest that the perturbative series is convergent within a given radius. Furthermore, the instanton action at fixed value of ’t Hooft coupling diverges. This implies that the corresponding singularity in the Borel plane moves to infinity, suggesting at least Borel summability. However, the expected singularity associated to infrared renormalons tHooft:1977xjm ; Beneke:1998ui does not move away and would induce a factorial growth of the coefficients similar to the finite $N$ case. An analysis of the perturbative expansion at high orders should settle this point. Another issue to be studied is the interplay between the order of perturbation theory and the value of $N$. The higher rate of growth in the number of non-planar diagrams would suggest that their contribution could overcome the $1/N^{2}$ suppression at sufficiently high orders. Something of this kind emerges when analyzing the simultaneous expansion in powers of ’t Hooft coupling and $1/N^{2}$ as seen in Ref. Marino:2012zq . The behaviour of the reduced model could however be quite different.

To close the circle, recently a new method has arisen that allows the numerical computation of high order coefficients in the perturbative expansion. It goes under the name of Numerical Stochastic perturbation theory (NSPT) and was pioneered by the works of Di Renzo and collaborators DiRenzo:1993hs ; DiRenzo:1994av ; DiRenzo:1994sy ; DiRenzo:1995qc ; Burgio:1997hc ; DiRenzo:2004hhl (for earlier developments on stochastic quantization Parisi:1980ys and stochastic perturbation theory, see also Damgaard:1987rr ; BOOK_SQ1 ; BOOK_SQ2 ). Recently, this methodology has led to remarkable results about high order coefficients in SU(3) gauge theories DiRenzo:2000ua ; Rakow:2005yn ; Bauer:2011ws ; Horsley:2012ra ; Bali:2013pla ; Horsley:2013pra ; Bali:2014fea ; DelDebbio:2018ftu ; DelDebbio:2018vhr . As the memory size required is proportional to the highest perturbation order involved in the computation in NSPT and to the volume, the authors of Bali:2013pla ; Bali:2014fea ; DelDebbio:2018ftu have used twisted boundary conditions aware that they lead to reduced finite volume dependence.

Stochastic quantization is based on the Langevin equation, and the first NSPT studies used the perturbative expansion of this equation. However, for lattice theories various other Monte Carlo algorithms, such as Kramers, hybrid molecular dynamics (HMD), and hybrid Monte Carlo (HMC) algorithms, have been developed and used. Recently, extensive studies have been carried out using NSPT algorithms DallaBrida:2017pex ; DallaBrida:2017tru based on them. In particular, the HMD based algorithm for NSPT is easy to implement by perturbatively expanding the HMD/HMC codes used for non-perturbative simulations. In this case, one can benefit from various numerical integration schemes for the molecular dynamics part to reduce the finite step size error of the integration. Thus, the NSPT algorithm based on the HMD/HMC algorithm could open the way to efficiently estimating very high order coefficients.

Our work is a first attempt to apply this methodology to matrix models. For the time being, our work is mostly exploratory but a necessary step before any attempt of a higher order and larger $N$ study. Nonetheless, apart from showing that the method works with an incredibly high precision for the low-order coefficients, we also present results which extend to higher order the previous analytical results Perez:2017jyq . In addition, our work provides interesting information about the probability distribution of the perturbative coefficients. In particular, we have studied the dependence of the cumulants of these distributions with respect to the different parameters: the value of $N$ and of the fluxes, the size of the loops and the order of perturbation theory. These results allow us to determine the necessary computational requirements for any further extension of these studies.

The layout of the paper is the following. In section 2, after introducing the TEK model, we explain the application of the NSPT algorithm based on the HMD algorithm to the TEK model. In Section 3 we present the numerical results for the perturbative coefficients of the Wilson loops. We employ $(N,k)=(16,1),(49,2),(121,3)$ for $\mathrm{SU}(N)$ and flux parameter $k$ of the TEK model. The coefficients are computed up to four-loop ($O(g^{8})$) and the first two-loop coefficients are compared to the analytic values. The probability distribution of these estimates is investigated and used to explore the feasibility for extending our results to higher orders. We estimate the numerical computational cost of NSPT for the TEK model at large values of $N$ and large perturbative orders in section 4. Possible improvements on the algorithm are also discussed. In the last section we summarize the main results of the paper. Two technical points are included in appendices.

2 TEK model and NSPT

In this section we start by briefly introducing the TEK model together with the gauge fixing method. Then we recall the Hybrid Molecular Dynamics (HMD) algorithm for nonperturbative simulations of the TEK model. This algorithm is then perturbatively expanded to derive the equation of motion for the NSPT algorithm.

The computational cost is estimated in terms of the highest order $N_{\mathrm{trunc}}$ in the perturbative expansion involved in the algorithm and the matrix size $N$ of $\mathrm{SU}(N)$.

2.1 TEK model

The partition function of the TEK model in four-dimensions is defined by

[TABLE]

where $U_{\mu}$ are $\mathrm{SU}(N)$ matrices. We choose the symmetric twist characterized by $N=\hat{L}^{2}$ and $z_{\mu\nu}$ given by

[TABLE]

where $k$ is an integer coprime with $\hat{L}$ ($=\sqrt{N}$). The bare coupling constant $g$ is defined through $\beta=1/g^{2}$.

The Wilson loop operator for an $R\times T$ rectangle in the $\mu$–$\nu$ plane is defined by

[TABLE]

In this paper we restrict ourselves to square loops $R=T$.

The analytic perturbative expansion of the observables proceeds by expanding the link matrices around the classical vacuum as

[TABLE]

where the matrices $\Gamma_{\mu}$ define the classical vacuum of eq. (2) and have the following property

[TABLE]

A particular solution (up to multiplication by a phase) is singled out by an appropriately chosen gauge fixing condition, accompanied by the corresponding ghost term.

In order to stabilize the runaway trajectories, NSPT also requires the use of the stochastic gauge fixing method Zwanziger:1981kg . Here we are using the functional for the Landau gauge condition, from which we can extract the stochastic gauge fixing contribution in NSPT DiRenzo:1993hs ; DiRenzo:1994av . The gauge fixing functional $F[G]$ is defined by

[TABLE]

where $G\in\mathrm{SU}(N)$ is the gauge transformation matrix. The Landau gauge condition is achieved by maximizing the functional, by application of several iterations of the form:

[TABLE]

with

[TABLE]

where $\alpha$ is a parameter. When this iteration process is transformed into a continuous evolution equation with a fictitious time $t$, we obtain the following differential equation:

[TABLE]

where $w(t)$ is an auxiliary Hermitian-traceless matrix tracing the steepest descent trajectory. This evolution equation can be combined with the NSPT process to implement the stochastic gauge fixing.

2.2 NSPT for the TEK model

NSPT is based on stochastic quantization Parisi:1980ys . This amounts to writing down a stochastic differential equation of the Langevin type for the field variables of a target system, which relaxes to the equilibrium probability density at large times Damgaard:1987rr ; BOOK_SQ1 ; BOOK_SQ2 . The NSPT method is obtained by expanding these field variables in powers of the coupling constant and casting the original Langevin equation into a tower of equations, one for each power. Observables that are analytic functions of the field variables acquire a corresponding perturbative expansion. Each term in the expansion becomes an stochastic variable whose mean value at large times gives us the corresponding coefficient that we are after.

In this paper we employ the HMD based NSPT, which has been introduced in refs. DallaBrida:2017tru ; DallaBrida:2017pex . This is easy to implement through modifications of existing codes of nonperturbative HMD or HMC algorithms, and is preferable to systematically improve the molecular dynamics (MD) part, reducing the error arising from the finite MD step size.

To explain the method in our case, we will begin by reviewing the generalized HMD (GHMD) algorithm for nonperturbative simulations. Then we will apply the perturbative expansion to the algorithm, to derive the HMD based NSPT algorithm.

The HMD partition function is introduced by adding the canonical momentum variable $P_{\mu}$ conjugate to $U_{\mu}$ to the original partition function eq. (1).

[TABLE]

The classical dynamics for $(P_{\mu},U_{\mu})$ with the Hamiltonian $H[P,U]$ reproduces the microcanonical ensemble at fixed energy $E=H[P,U]$. Subsequent refreshment of the variable $P_{\mu}$ with the Gaussian distributions generates the corresponding canonical ensemble $e^{-H[P,U]}$. The marginal distribution for $U_{\mu}$ becomes the $e^{-S[U]}$ distribution that we are looking for.

The HMD algorithm proceeds as described in Algorithm 1, where the equations in Step 2 are discretized in time. This discretization violates the exact energy conservation, which distorts the distribution shape. In order to ensure that the discretized Markov chain leads to the correct probability distribution the Monte Carlo process must satisfy time reversibility and area preservation in the MD evolution, for which the leapfrog type scheme is normally used. For nonperturbative simulations, the Metropolis test can be inserted between Step 1 and 2 to compensate the violation of energy conservation, yielding the HMC algorithm.

The full momentum refreshment in Step 1 will cause random walking in phase space so that the autocorrelation times becomes longer for the ensemble. To relax random walking behaviour, the generalized HMC algorithm (GHMC) has been introduced in ref. Kennedy:2000ju . For the GHMC algorithm the momentum is partially refreshed in Step 1 as

[TABLE]

instead of eq. (17), where $c_{1}$ is a mixing parameter. $P_{\mu}$ in the right hand-side is the solution of Step 2. A momentum reflection step is added after the rejection of the Metropolis test in the GHMC algorithm, while this step is absent in the GHMD algorithm. We can further include the gauge fixing process eq. (13) into the MD evolution equation.

With the leapfrog algorithm and the partial momentum refreshment, the updating algorithm for one step having a small $\Delta t$ interval, is given by

[TABLE]

where $(P_{0,\mu},U_{0,\mu},w_{0})$ is the initial state and $(P_{1,\mu},U_{1,\mu},w_{1})$ is the final state. Eq. (25) corresponds to the leapfrog evolution for $\Delta t$, for which various higher-order schemes are available. For simplicity we will explain the second order leapfrog scheme only. The gauge fixing evolution (12)–(14) is interleaved into the MD evolution Rossi:1987hv ; Davies:1987vs and the transformation in eq. (26) does not affect the gauge invariant observables, however, we introduce this to explain the stochastic gauge fixing for NSPT Zwanziger:1981kg ; DiRenzo:1993hs ; DiRenzo:1994av . Taking the $\Delta t\to 0$ limit with $c_{1}=e^{-\gamma\Delta t}$, as shown in ref. DallaBrida:2017tru , this evolution reduces to

[TABLE]

where $D_{\mu}w$ is defined as

[TABLE]

and $\zeta_{\mu}$ is a random noise satisfying

[TABLE]

The terms with $w$ act as the gauge damping force. Eq. (27) corresponds to the Kramers equation and in the limit $\gamma\to+\infty$ ($c_{1}=0$) to the Langevin equation with gauge damping force.

In the GHMD algorithm Step 2 is obtained by repeating eqs. (25) and (26) until the total evolution time becomes a fixed value, and then applying eq. (24) to implement Step 1. When the evolution time in Step 2 reduces to $\Delta t$, the method becomes Langevin ($c_{1}=0$) or Kramers ($c_{1}\neq 0$) algorithm.

Having explained the GHMD algorithm, we now describe the corresponding NSPT algorithm. This follows by expanding $(P_{\mu},U_{\mu},w)$ and the MD equation (22) or eqs. (24)–(26) in a power series in $g$. We will now describe the perturbative expansion for eqs. (24)–(26), because it is sufficient to write the simulation program. In order to simplify notation, we first define the $\star$-product as the convolution product of two perturbative series, given by

[TABLE]

where $A,B,C$ are matrices and $A^{(k)},B^{(k)},C^{(k)}$ are the coefficient matrices.

The perturbative expansion for $(P_{\mu},U_{\mu},w)$ is defined by

[TABLE]

where $U_{\mu}^{(0)}=\Gamma_{\mu}$ is kept fixed in NSPT as it is the perturbative vacuum. Rescaling the fictitious time and the gauge fixing parameter as $t^{\prime}=t/g$ and $\alpha^{\prime}=g^{2}\alpha$, and substituting eqs. (33)–(35) into eqs. (24)–(26), we obtain

[TABLE]

for each perturbative order $k=1,2,\dots$. The details of the perturbative expansion of the matrix exponential $\exp[A]=I+g\left(\exp[A]\right)^{(1)}+g^{2}\left(\exp[A]\right)^{(2)}+\dots$ with $A=gA^{(1)}+g^{2}A^{(2)}+\dots$ are explained in Appendix A. The perturbative expressions for $F^{(k)}_{\mu}$ and $\Theta^{(k)}$ are extracted from eqs. (21)–(22), and (10)–(11), respectively:

[TABLE]

In the limit $\Delta t^{\prime}\to 0$, the equation of motion should conserve the energy eq. (16) order by order in perturbation theory. This energy conservation can be monitored during the simulation.

The Wilson loop operator eq. (4) is expanded similarly:

[TABLE]

where we used a shorthand notation for matrix power with the $\star$-product such as $(A\star A\star A)\equiv A^{\star 3}$. The expectation values at large times of our $W^{(k)}_{\mu\nu}(R,T)$ yield the coefficients of the perturbative expansion of Wilson loops. For the latter we take the notation given in ref. Perez:2017jyq where the first two coefficients have been computed analytically. If we consider for example an $R\times T$ Wilson loop in the $\mu$–$\nu$ plane, the relation is as follows

[TABLE]

Notice that, since the perturbative expansion is in powers of ’t Hooft coupling $\lambda=g^{2}N$, the expectation value of $W^{(k)}_{\mu\nu}$ for odd values of $k$ has to vanish.

A technical point is the necessity to correct for deviations from the unitarity constraint induced by the numerical round-off error from the finite precision of computer arithmetic. While imposing the traceless-Hermitian character of $(P_{\mu}^{(k)},w^{(k)})$ is easily done, the conditions on $\{U_{\mu}^{(k)}\}$ following from the unitarity of the link matrices can be imposed by the matrix logarithm scheme for reunitarization DiRenzo:2004hhl . The details of the group projection is described in Appendix B.

Let us conclude by summarizing our NSPT algorithm. As explained earlier, the finite time step induces a distortion in the probability distribution, which in NSPT cannot be corrected by a Metropolis test. Hence, it is preferable to employ a higher order integration scheme in the MD evolution. Although, for simplicity, we used the simple (second order) leapfrog scheme to explain the HMD based NSPT algorithm, in practice we employed the fourth order leapfrog scheme for eq. (36), while the evolution of $w$ is not changed. The properties of second and fourth order Omelyan-Mryglod-Folk schemes OMF ; Takaishi:2005tz , and second order Runge-Kutta scheme in NSPT have been investigated in DallaBrida:2017pex ; DallaBrida:2017tru . We summarize the NSPT algorithm used in this paper in algorithm 2. This corresponds to the Kramers type NSPT algorithm (called KSPT in DallaBrida:2017pex ). The specific parameters are fixed to $\gamma=0.5$ for the momentum refreshing parameter $c_{1}=\mathrm{e}^{-\gamma\Delta t^{\prime}}$ and $\alpha^{\prime}=2$ for the gauge damping parameter. The trajectories are performed over a fixed time $t^{\prime}=1$, with a perturbative reunitarization step at the end. The Wilson loop coefficients are computed at the end of every trajectory.

2.2.1 Computational estimates

The maximum power of $g$ that is studied $N_{\mathrm{trunc}}$ is limited by the total computer memory available. For the $\mathrm{SU}(N)$ TEK model, the total memory requirement is of $O(N_{\mathrm{trunc}}N^{2})$. Most of the computational time is spent in matrix multiplication, whose cost is of $O(N^{3})$. The cost of the convolutional $\star$-product is of $O(N_{\mathrm{trunc}}^{2})$ with a naive implementation. The evaluation of the perturbative matrix function requires one more factor of $N_{\mathrm{trunc}}$, yielding the cost of $O(N_{\mathrm{trunc}}^{3})$ (see Appendix A). The cost of the reunitarization we implemented is of $O(N_{\mathrm{trunc}}^{4})$ (see Appendix B).

The energy difference in a trajectory $\Delta H$ is non-zero for a finite MD step size $\Delta t$. The relation between the energy conservation violation $\Delta H$ and $\Delta t$ depends on the MD integration scheme. We assume that algorithm 2 with 4th order leapfrog scheme yields $\Delta H^{(k)}\sim\Delta t^{4}$ at each perturbation order DallaBrida:2017pex , where $\Delta H^{(k)}$ is the energy conservation violation for the perturbative coefficient of the total energy in NSPT. Since the energy is proportional to $N^{2}$, to achieve a constant energy conservation violation, the number of steps $N_{\mathrm{MD}}$, therefore, should scale as $N_{\mathrm{MD}}=\sqrt{N}$.

The total computational cost of the NSPT algorithm truncated at $N_{\mathrm{trunc}}$ for the $\mathrm{SU}(N)$ TEK model, then, scales as $O(N_{\mathrm{trunc}}^{3}N^{3}N_{\mathrm{MD}})=O(N_{\mathrm{trunc}}^{3}N^{7/2})$ for the MD part and $O(N_{\mathrm{trunc}}^{4}N^{3})$ for the reunitarization part. Our estimates of total CPU time do not take into account autocorrelation times. For that we need experimental studies on the statistical properties of the Monte Carlo data, which depends on models and algorithms. As far as the parameters we investigated, no sizable autocorrelation and parameter dependence are observed and the autocorrelation length is of $O(1)$ in units of the trajectory length $t^{\prime}=1$.

3 Numerical results

In this section we show the results for the perturbative coefficients of the Wilson loops with NSPT, and investigate the statistical properties of the distribution. We have implemented the NSPT algorithm as explained in the previous section. We accumulate the statistics for the perturbative coefficients of Wilson loops up to order $g^{8}$ ($N_{\mathrm{trunc}}=8$). The mean value gives us the corresponding perturbative quantity, but the variance and higher cumulants allow us to estimate the required statistics to achieve a given error in these coefficients. Fortunately, we have analytic results for the one and two-loop coefficients to test our results at these orders, but we can extend these results two more orders in powers of $\lambda$. Furthermore, we monitor the coefficients at odd-orders $O(g^{2\ell+1})$ which should be zero for Wilson loops. This provides a test of the cancellation of the non-loop effect in NSPT. Indeed, our results are consistent with zero within the two standard deviation. Therefore, we will concentrate in giving the results of even-order coefficients.

The parameters and statistics are shown in table 1. $N_{\mathrm{MD}}$ is the number of MD steps for unit trajectory $t^{\prime}=1$. Statistical errors are estimated with the jackknife method after binning in 1000 trajectory samples. We discard the first $\geq 1000$ trajectories to account for thermalization. Since we employ the 4th order leapfrog scheme for the MD integrator, a ${\Delta t^{\prime}}^{4}=1/N_{\mathrm{MD}}^{4}$ dependence is expected in observables DallaBrida:2017pex . We also study the limit of vanishing step size for some representative cases.

Given the pilot nature of our study we have concentrated in studying cases which have been analyzed in detail in the analytic calculations of ref. Perez:2017jyq . Thus we concentrated in the three values of $N=16,49,121$. The lowest values are not so interesting from the point of view of approximating large $N$ Yang-Mills theory at infinite volume. Corrections are expected to be large. However, only for a value as low as $N=16$ can one see clearly two effects of the twisted boundary conditions: the breakdown of CP and of cubic-rotational invariance. The first phenomenon manifests itself in non-vanishing imaginary parts of Wilson loops. We included this case to test this phenomenon in our NSPT determinations. On the other hand having $N=49$ and $121$ allows us to test the dependence on $N$ of both the physical and computational parameters. We also tested several values of the flux parameter $k$. For sufficiently large values of $k/\hat{L}$ the dependence is quite small. For large values of $N$ and small values of $k$ the dependence can be quite strong even in the perturbative calculation. Non-perturbatively this restriction is even more important to preserve a remnant of center symmetry GonzalezArroyo:2010ss that validates reduction in the limit $N\to\infty$.

3.1 Comparison with analytic calculations

As mentioned earlier the presence of the twist breaks part of the symmetries of the standard lattice symmetries with periodic boundary conditions. In particular, part of the cubic group is broken. This translates into a dependence of the coefficients of the Wilson loop perturbative expansion on the plane in which the Wilson loop lies. Due to the remaining symmetry for our choice of twist, there are two sets of planes which we label as $S_{1}=\{(\mu,\nu):(1,2),(2,3),(3,4),(4,1)\}$ and $S_{2}=\{(\mu,\nu):(1,3),(2,4)\}$ Perez:2017jyq , such that the results for all $\mu-\nu$ planes contained in each $S_{i}$ should be the same. The results for all planes in $S_{1}$ need not be equal to those contained in $S_{2}$. To increase the statistics we can average all the planes within each set.

Our results for $\hat{W}^{11}_{\ell}$ and $\hat{W}^{33}_{\ell}$ at $(N,k)=(16,1)$ are shown in figures 1 and 2, respectively. The mean values are plotted as a function of $1/N_{\mathrm{MD}}^{4}$ and extrapolated linearly to vanishing step size. The linear dependence in this variable is the expectation for our 4th order leapfrog scheme for the MD integrator, and our data are consistent with this expectation. The data for the $S_{1}$ and $S_{2}$ planes are plotted and extrapolated separately (purple dashed : $S_{1}$, green dotted : $S_{2}$). The analytic results, depicted by horizontal lines in the same plots, only predict differences among planes for $\hat{W}_{2}$. At all orders, the cubic symmetry breaking in the real parts is found to be small and of the order of the errors of our calculation. As mentioned earlier the imaginary parts beyond the leading order are non-zero and different for the two families of planes. Our results reproduce both features and match nicely with the analytic results for the first two coefficients.

All our results for $(N,k)=(16,1)$ and $(49,2)$, extrapolated linearly in $1/N_{\mathrm{MD}}^{4}$ to zero, are collected in Tables 2 and 3. We also tabulate the analytic values from ref. Perez:2017jyq and $\chi^{2}/\mathrm{DoF}$ from the extrapolating fit. The violation of CP and cubic invariance due to the twist is well seen for the $N=16$ case at the two loop level and the differences between $S_{1}$ and $S_{2}$ averages are consistent with the analytic values. For the $(N,k)=(49,2)$ case, the results of our NSPT analysis are also consistent with analytic values. However, the violation of CP and cubic invariance is too small to be seen in comparison with the statistical errors. Since the violation of CP and cubic invariance disappears in the large $N$ limit, we conclude that in practice one would not be able see any effect of this symmetry breaking for even larger values of $N$. Hence, hereafter, we show the coefficients averaged over all $\mu$–$\nu$ planes ($S=S_{1}+S_{2}$) in order to study the $k$ dependence and $N$ dependence. We emphasize the small errors of our coefficient determinations, typically of order $10^{-5}$. This gives about 3 to 5 significant digits in some determinations.

We next investigate the dependence on the flux parameter $k$ at $N=49$. For that purpose we studied the $k=1$ and $k=3$ results at fixed value of $N_{\mathrm{MD}}=32$. The results for the coefficients for the plaquette and $4\times 4$ Wilson loop are displayed in figure 3. The results for $k=2$ are given both at $N_{\mathrm{MD}}=32$ as well as extrapolated $1/N_{\mathrm{MD}}\to 0$ (open triangles). From the data we conclude that $N_{\mathrm{MD}}=32$ is sufficient for the current analysis as the statistical error and the systematic error from the finite MD time step are comparable. The horizontal dash lines in $\hat{W}^{RR}_{1}$ and the filled squares in $\hat{W}^{RR}_{2}$ correspond to the analytic values. For the plaquette, the observed $k$-dependence is compatible with statistical errors. As expected, the differences are small but well beyond statistical errors for the $4\times 4$ loop. This is also the case for the third order coefficient $\hat{W}^{RR}_{3}$ for which we have no analytic results. Anyhow, no dramatic $k$ dependence is seen for the third and fourth order coefficients.

Our results averaged over all planes are displayed in table 4. The choice of values $(N,k)=(16,1),(49,2),(121,3)$ are selected so as to keep $k/\hat{L}$ at $\sim 0.27$ ($k/\hat{L}=1/4,2/7,3/11$ for each data set). The results at the one and two-loop level are consistent with the analytic values within two standard deviations. We also display (NPT fit) an estimate of the three-loop coefficients for $N=16$ and $49$ obtained by fitting nonperturbative expectation values to a third order polynomial in $\lambda$ with the first coefficients fixed to the analytic values Perez:2017jyq . The results are roughly consistent with our NSPT determinations, which are expected to be more reliable.

Since the main interest of the TEK model is its large $N$ limit, which coincides with Yang-Mills at infinite $N$ and infinite volume, we give the infinite $N$ extrapolation of our results. As mentioned previously finite $N$ corrections are expected to grow for large values of $R/\hat{L}$. Thus, we exclude the data of $(N,k)=(16,1)$ and base our analysis on $(N,k)=(121,3)$ and $(49,2)$ only and assuming a $1/N^{2}$ dependence of the coefficients. The results are displayed in figure 4 and the numerical values of the extrapolated coefficients are tabulated in table 5. In the figure we only show the plaquette and $4\times 4$ Wilson loop up to three-loop order. The error includes only the statistical one. The values extrapolated to $N\to\infty$ are shown as filled circles, while the horizontal solid lines are the corresponding analytic values DiGiacomo:1981lcx ; Alles:1993dn ; Alles:1998is ; Perez:2017jyq ; privcomm . We see a reasonable good agreement within errors. The deviations for the $4\times 4$ loop might include also a systematic error due to the relatively large values of $4/\hat{L}$. For a more precise control of the systematic errors in the extrapolated results it would be necessary to have at least 3 sufficiently large values of $N$, which is left out of this exploratory work.

3.2 Analysis of the distribution

For the purpose of determining the statistical requirements involved in a precise determination of the perturbative coefficients we performed an analysis of the probability distribution of the corresponding NSPT estimates. As a quantitative measure we computed the cumulants up to the 5th cumulant. The $N$ dependence of the variance (second cumulant) is important since the statistical error at fixed statistics is proportional to the square root of this variance. The non-zero value of higher order cumulants measures the deviation from the normal distribution. This is of practical importance since, as seen in ref. Alfieri:2000ce , the distribution for perturbative coefficients at higher order in NSPT could show strong deviations from a normal distribution, the so-called “Pepe effect”. We show the histograms for $\hat{W}^{11}_{\ell}$ and $\hat{W}^{33}_{\ell}$ in figures 5 and 6, respectively. We have fixed $N_{\mathrm{MD}}=32$ for better comparison. The distributions are plotted in logarithmic scale for which the normal distribution corresponds to a parabola. The dotted lines are the best fit to a normal distribution. One can see deviations in the tail for some plots for $N=16$, but it disappears for larger $N$. Hence, the problem of “Pepe effect” is absent in the TEK model for sufficiently large $N$, similarly to what happens for pure SU(3) lattice gauge theory. Since the 3rd–5th cumulants have minor effect or are not detectable on the distribution, we focus on the variance only in the following analysis.

We have studied the dependence of the variance of $\hat{W}_{\ell}^{RR}$ on $R$, $\ell$ and $N$. Notice that we can also use the results for half-integer $\ell$ in this analysis. The $R$ dependence is illustrated in figure 7, and is seen to go linearly with $R^{4}$. On the other hand the dependence on $j=2\ell$, displayed in figure 8 shows an exponential behaviour for sufficiently large $j$. Altogether, an ansatz of the form

[TABLE]

with $c=0.0013$ (solid lines in these figures) describes our data quite well. For large $N$ this corresponds to a dependence on the ratio of the loop size with respect to the effective size of the box $R/\hat{L}$ to the fourth power, similar to the observed subleading corrections in perturbation theory.

The $1/N^{2}$ dependence of the variance follows from the perturbative realization of factorization $\langle W^{2}\rangle=\langle W\rangle^{2}+O(1/N^{2})$. Similar arguments lead to the prediction that the $p$-th cumulant scales as $O(1/N^{2p-2})$. This by itself explains why the distribution tends to normal at large $N$. We cannot check this behaviour of the higher cumulants since their values are too small.

The exponential dependence of the variance on the order $4^{\ell}$ coincides precisely with the characteristic growth in the number of planar diagrams Brezin:1977sv ; Koplik:1977pf ; Marino:2012zq , given by the Catalan number. In the absence of renormalons, a similar growth is also expected for the mean values $\hat{W}_{\ell}$, but it is hard to check this behaviour with our results reaching only up to order $\ell=4$. Indeed, as observed in ref. Bali:2014fea , extremely large order perturbative coefficients are needed to detect the expected factorial behavior for the pure SU(3) lattice gauge theory. Hence, it would be interesting to extend our results to higher order to explore this point. The computational requirements will be studied in the next subsection.

4 Estimate of the computational cost at high order and large $N$

Here we will evaluate the computer requirements to extend our calculation $\hat{W}_{\ell}$ up to order $\ell=N_{\mathrm{trunc}}/2$. The precision in the determination is proportional to the standard deviation of the corresponding probability distribution. Hence, if we want to keep this precision fixed our data sample should grow with the square, the variance. We previously estimated that for each sample point the computational cost goes as $O(N^{7/2}N_{\mathrm{trunc}}^{3})$ for the MD part and as $O(N^{3}N_{\mathrm{trunc}}^{4})$ for the reunitarization part. Using eq. (159), we can then give an estimate of the total computational cost as follows:

[TABLE]

or

[TABLE]

This estimate can be reduced by improving the algorithm. To obtain perturbative coefficients at very higher order, the exponential scaling behaviour, $2^{N_{\mathrm{trunc}}}$, has to be relaxed by reducing the variance using a method such as a kind of reweighting method in NSPT. Furthermore, the factor $N_{\mathrm{trunc}}^{4}$ comes from our implementation of the reunitarization algorithm, for which a better algorithm could be applicable. One possibility is to use the perturbative Gram-Schmidt algorithm for the reorthogonalization followed by the perturbative subtraction of the U(1) phase. The estimation of the U(1) phase still requires the perturbative trace-log computation whose cost is of $O(N_{\mathrm{trunc}}^{3})$. At least one can reduce the factor from $N_{\mathrm{trunc}}^{4}$ to $N_{\mathrm{trunc}}^{3}$. In that case eq. (160) would dominate the total cost.

The cost of $O(N_{\mathrm{trunc}}^{2})$ in the convolutional product of two series could be ameliorated by using the fast Fourier transformation (FFT) algorithm, which reduces the cost to $O(N_{\mathrm{trunc}}\log N_{\mathrm{trunc}})$. Improvement in this direction can be used in future studies.

5 Summary

In this paper we have developed and applied an HMD-based NSPT algorithm to compute the perturbative expansion of Wilson loops in the TEK model. At large $N$ these coefficients coincide with those of large volume Yang-Mills theory, hence their interest. The low order coefficients are very precise and match with the values obtained by the analytic calculation of ref. Perez:2017jyq . We have been able to extend the perturbative series two more orders, up to $O(\lambda^{4})$.

We also studied the statistical properties of the coefficients in NSPT. We found that their distribution tends to normal at large $N$ and there is no Pepe effect. Furthermore, we have determined the dependence of the corresponding variance on $N$ and the order of perturbation theory. This has allowed us to estimate the computational requirements necessary to extend our results to higher orders. This would allow extending studies of the type described in ref. Marino:2012zq to Yang-Mills theory in the limit of infinite number of colours.

Acknowledgements.

We thank Margarita García Pérez for discussions and for providing us with the analytic values of some perturbative coefficients, and Alberto Ramos for his encouraging comments on the draft. A.G.-A. acknowledges financial support from the MINECO/FEDER grant FPA2015-68541-P and the MINECO Centro de Excelencia Severo Ochoa Program SEV-2016-0597. K.-I.I. and M.O. are supported by JSPS KAKENHI Grant Numbers 16K05326 and 17K05417, respectively. K.-I.I. and I.K. acknowledge for financial support by Priority Issue 9 to be tackled by using Post K Computer. The numerical computation for this paper were performed on the subsystem A of the ITO supercomputer system at Research Institute for Information Technology, Kyushu University, and workstations at INSAM (Institute for Nonlinear Sciences and Applied Mathematics), Hiroshima University.

Appendix A Perturbative matrix functions

In this appendix we describe the details of evaluating the perturbative coefficients of a matrix function with a perturbative series as the argument. After introducing the generic algorithm for perturbative matrix functions, we explain the cases with matrix exponential and matrix logarithm explicitly.

We consider a matrix function $f(F)$ whose argument $F$ is a matrix. We assume that $f(F)$ has the following Taylor expansion form:

[TABLE]

We also assume that the perturbative expansion of $F$ is

[TABLE]

where $g$ is the coupling constant and $F^{(k)}$ are the coefficient matrices. Note that the leading term of $F$ is $O(g)$. We truncate the perturbative expansion at $N_{\mathrm{trunc}}$-th order.

[TABLE]

The perturbative expansion of $f(F)$ with eq. (164) is also truncated at $N_{\mathrm{trunc}}$ as

[TABLE]

We would like to know the coefficient $\left[f(F)\right]^{(k)}$ in terms of $c_{k}$’s and $F^{(k)}$’s. A recursion relation to evaluate the coefficient $\left[f(F)\right]^{(k)}$ can be obtained by expanding the Horner’s method for polynomials as follows.

As $F$ is truncated at $N_{\mathrm{trunc}}$-th order, the Taylor expansion is also truncated at $N_{\mathrm{trunc}}$-th order,

[TABLE]

where the matrix $G$ is introduced for the truncated form. $G$ can be evaluated via the Horner’s method:

[TABLE]

This form leads to the following recursion relation:

[TABLE]

where $G_{j}$ are working area. As $F$ has the series form eq. (163), $G_{j}$ also has the series form:

[TABLE]

Substituting the series form (164) for $F$ and (172) for $G_{j}$ into the recursion relation, we obtain

[TABLE]

Here we discard higher order terms $O(g^{k})$ with $k>N_{\mathrm{trunc}}$.

Consequently the perturbative coefficient $[f(F)]^{(k)}$ of $f(F)$ can be obtained as algorithm 3. The computational cost scales with $N_{\mathrm{trunc}}^{3}$.

A.1 Matrix exponential for updating $U$

Here we consider the case of $f(F)=\exp(F)=\exp(iP\Delta\tau)$ with $P=\sum_{k=1}^{\infty}\beta^{-k/2}P^{(k)}$. The Taylor expansion for $f(F)$ is

[TABLE]

Thus the coefficients are $c_{k}=(i\Delta\tau)^{k}/{k!}$. We identify $F=P$, $F^{(k)}=P^{(k)}$, $d_{j}=i\Delta\tau/j$, and $[\exp(iP\Delta\tau)]^{(k)}=G_{1}^{(k)}$ with $g=\beta^{-1/2}$ for the recursion algorithm 3 (omitting the last step lines 9–11 as $c_{0}=1$).

A.2 Matrix logarithm for reunitarization

Next we consider the case of $f(F)=\ln(I+F)=\ln(U)$ with $U=I+\sum_{k=1}^{\infty}\beta^{-k/2}U^{(k)}=I+F$. The Taylor expansion for $f(F)$ is

[TABLE]

Thus the coefficients are $c_{0}=0$, $c_{\ell}=(-1)^{\ell+1}/\ell$. As $F$ is expanded as $F=\sum_{k=1}^{\infty}\beta^{-k/2}U^{(k)}$, we can identify $F^{(k)}=U^{(k)}$, and $[\ln(U)]^{(k)}=G_{1}^{(k)}$ with $g=\beta^{-1/2}$ and $d_{j}=-(j-1)/j$ for the recursion algorithm 3.

Appendix B Perturbative reunitarization

Reunitarization of $U$ is applied perturbatively on $U^{(k)}$ via the perturbative expansion of the matrix logarithm:

[TABLE]

where $A$ should be anti-Hermitian and traceless for $U$ to be $\mathrm{SU}(N)$. From the perturbative expansion of the matrix logarithm and $A=\sum_{k=1}^{\infty}g^{k}A^{(k)}$, we have

[TABLE]

where $\left[\ln[U]\right]^{(k)}$ is described in subsection A.2. The $\mathrm{SU}(N)$ condition on the perturbative coefficients $A^{(k)}$ is

[TABLE]

By inspecting algorithm 3, we can find that the coefficient $\left[\ln[U]\right]^{(k)}=A^{(k)}$ has the following dependency on $U^{(k)}$:

[TABLE]

where $X[U^{(1)},U^{(2)},\dots,U^{(k-1)}]\equiv\left[\ln[U]\right]^{(k)}-U^{(k)}$. Thus the $\mathrm{SU}(N)$ condition on $A^{(k)}$ can be guaranteed by applying algorithm 4 on $U^{(k)}$. The computational cost of algorithm 4 is $O(N_{\mathrm{trunc}}^{4})$.

In order to apply algorithm 4 to the TEK model, we have to reunitarize the matrix $V^{(k)}_{\mu}=U^{(k)}_{\mu}\Gamma^{{\dagger}}_{\mu}$ as the perturbative vacuum is $U^{(0)}_{\mu}=\Gamma_{\mu}$. After reunitarizing $V_{\mu}^{(k)}$, $U^{(k)}_{\mu}=\Gamma_{\mu}V^{(k)}_{\mu}$ is computed.