Universal Renormalons in Principal Chiral Models

TL;DR

This paper numerically verifies the presence of renormalons in two-dimensional sigma models, demonstrating their role in nonperturbative physics and asymptotic series behavior similar to gauge theories.

Contribution

First numerical verification of renormalon effects in 2D sigma models, linking perturbative ambiguities to nonperturbative phenomena in these models.

Findings

High order coefficients are independent of model rank.

Ratios of consecutive coefficients follow renormalon asymptotics.

Individual coefficients decrease strongly, yet renormalon behavior is evident.

Abstract

Perturbative expansions in many physical systems yield 'only' asymptotic series which are not even Borel resummable. Interestingly, the corresponding ambiguities point to nonperturbative physics. We numerically verify this renormalon mechanism for the first time in two-dimensional sigma models, that, like four-dimensional gauge theories, are asymptotically free and generate a strong scale through dimensional transmutation. We perturbatively expand the energy through a numerical version of stochastic quantization. In contrast to the first energy coefficients, the high order coefficients are independent on the rank of the model. Technically, they require a sophisticated analysis of finite volume effects and the continuum limit of the discretized model. Although the individual coefficients do not grow factorially (yet), but rather decrease strongly, the ratio of consecutive coefficients clearly obey the renormalon asymptotics.

Tables

Table 1. Table 1: Lattice setups for the Langevin runs. The simulations use five different stochastic time step sizes ϵ = 0.01 , 0.00875 , 0.0075 , 0.00625 , 0.005 italic-ϵ 0.01 0.00875 0.0075 0.00625 0.005 \epsilon=0.01,0.00875,0.0075,0.00625,0.005 , except for the geometries printed in bold face, where ϵ = 0.01 , 0.0075 , 0.005 italic-ϵ 0.01 0.0075 0.005 \epsilon=0.01,0.0075,0.005 .

Order $M$	Rank $N$	Lattice extensions $L$ $(V=L\times L$)
$20\left({\beta }^{-10}\right)$	$12$	$8,12,16,20,24,32,\mathrm{𝟒𝟖}$
$40\left({\beta }^{-20}\right)$	$3,4,5,6$	$8,12,16,20,24,32,\mathrm{𝟒𝟖}$

Table 2. Table 2: Range of τ 𝜏 \tau where the simulations are thermalized and unitarity violations are small.

$V$	$\tau$ range	Remarks
$16\times 16$	$100\text{–}500$	Only $N=3$ and w/o $ϵ=0.01,0.0625$
$20\times 20$	$\mathrm{\hspace{0.33em}\hspace{0.33em}75}\text{–}300$	$N=6$ w/o $ϵ=0.0075$
$24\times 24$	$150\text{–}400$
$\begin{array}{c}32\times 32\\ 48\times 48\end{array}$	$\begin{array}{c}170\text{–}500\\ 220\text{–}500\end{array}\}$	No unitarity violation found

Equations

$$E=\sum _{n=0}^{\infty }{c}_n{\alpha }^n(\alpha \to 0)\text{with }{c}_n\stackrel{n\to \infty }{\sim }{\gamma }^n{n}^{\kappa }n!.$$

$$\mathcal{S}$$

$$s_{\mu }(x)$$

$$a^{2}E=1-\frac{1}{4N^{2}V}\frac{\partial \log Z}{\partial \beta }=1-\frac{1}{4NV}⟨\sum _{x,\mu }{s}_{\mu }(x)⟩,$$

$$U(x)$$

$$a^{2}E$$

$$a{\Lambda }_L$$

$$=(\beta /{\beta }_0{)}^{{\beta }_1/{\beta }_0^2}\exp (-\beta /{\beta }_0),$$

$${\beta }_0=\frac{1}{8\pi },{\beta }_1=\frac{1}{128{\pi }^2}.$$

$$E_n$$

$$n^{\ast }\approx 16\pi \approx 50.$$

$$r_n:=\frac{E_n}{E_{n-1}n}$$

$$a^{2}E=\sum _{n=0}^{\infty }{E}_n{\beta }^{-n}-\frac{1}{L^2}\sum _{k=0}{F}_k{\beta }^{-k},$$

$$U=\text{mathbbm}1+U_1{\beta }^{-\frac{1}{2}}+U_2{\beta }^{-1}+\cdots +U_{2M}{\beta }^{-M},$$

$$U+U^{\prime }$$

$$U\cdot U^{\prime }$$

$$\frac{\partial \varphi (x,\tau )}{\partial \tau }=-\frac{\delta \mathcal{S}[\varphi (x,\tau )]}{\delta \varphi (x,\tau )}+\bar{\eta }(x,\tau ),$$

$${⟨\bar{\eta }(x,\tau )⟩}_{\bar{\eta }}$$

$${⟨\bar{\eta }(x,\tau )\bar{\eta }(x^{\prime },{\tau }^{\prime })⟩}_{\bar{\eta }}$$

$$\varphi (x,{\tau }_{n+1})=\varphi (x,{\tau }_n)-F_x[\varphi ,\eta ],$$

$$F_x[\varphi ,\eta ]:=ϵ\frac{\delta \mathcal{S}[\varphi (x,\tau )]}{\delta \varphi (x,{\tau }_n)}-\sqrt{ϵ}\eta (x,{\tau }_n)$$

$${⟨\eta (x_i,{\tau }_n)\eta (x_j,{\tau }_m)⟩}_{\eta }=2{\delta }_{x_i,x_j}{\delta }_{{\tau }_n,{\tau }_m}.$$

$$\bar{\mathcal{S}}[\varphi ]=\mathcal{S}[\varphi ]+ϵS_1[\varphi ]+\mathcal{O}({ϵ}^2).$$

$$U(x,{\tau }_{n+1})=\exp \left(-i\sum _a{T}^a{F}_x^a\right)U(x,{\tau }_n),$$

$$F_x^a=ϵ{\nabla }_x^a\mathcal{S}[U]+\sqrt{ϵ}{\eta }^a$$

$$f(e^{i\sum _a{T}^a{\omega }^a}U)=f(U)+\sum _a{\omega }^a{\nabla }^{a}f(U)+\mathcal{O}({\omega }^2).$$

$${⟨\eta (a,x_i,{\tau }_i)\eta (b,x_j,{\tau }_j)⟩}_{\eta }=2{\delta }_{x_i{x}_j}{\delta }_{{\tau }_i{\tau }_j}{\delta }_{ab}.$$

$$\begin{split}\nabla_{x}^{a}\mathcal{S}[U]&=-i\beta N\sum\limits_{\mu}\operatorname{Tr}\Big{(}T^{a}U(x)U^{\dagger}(x+\mu)\\ &\qquad\qquad-U(x+\mu)U^{\dagger}(x)T^{a}\\ &\qquad\qquad-U(x-\mu)U^{\dagger}(x)T^{a}\\ &\qquad\qquad+T^{a}U(x)U^{\dagger}(x-\mu)\Big{)}\end{split}$$

$$T_{ij}^a{T}_{kl}^b=\frac{1}{2}\left({\delta }_{il}{\delta }_{jk}-\frac{1}{N}{\delta }_{ij}{\delta }_{kl}\right)$$

$$\sum\limits_{a}\nabla_{x}^{a}\mathcal{S}[U]T^{a}=-i\frac{\beta N}{2}\sum\limits_{\mu}\Big{(}M_{\mu}(x)-M_{\mu}^{\dagger}(x)\\ -M_{\mu}(x-\mu)+M_{\mu}^{\dagger}(x-\mu)\Big{)},$$

$$U=e^{iA}$$

$$A_k\to \frac{1}{2}\left(A_k+A_k^{†}\right),A_k\to A_k-\frac{\mathrm{Tr}(A_k)}{N}\cdot \text{mathbbm}{1}_n.$$

$${\Delta }_k=\{\begin{array}{ll}|\sqrt{N}-\Vert V_k{\Vert }_F|/N& k=0\\ \Vert V_k{\Vert }_F/N& k>0\end{array},$$

Full text

Perturbation theory, the expansion in a small parameter, is a straightforward approach to many physical systems, both classical and quantum. The high-order behavior of the perturbative expansion, however, can be very elaborate: it may not be a convergent, rather an asymptotic series. This situation already occurs for relatively simple quantum mechanical systems, like the anharmonic oscillator expanded around zero anharmonicity, and even for ordinary integrals (as toy models for path integrals) Le Guillou and Zinn-Justin (1990). In quantum field theories (QFTs) the asymptotic nature of an expansion in powers of the coupling constant $\alpha$ hints at a physical instability, as pointed out by Dyson for quantum electrodynamics Dyson (1952).

The typical dependence of an observable $E$ on $\alpha$ contains factorially growing coefficients,

[TABLE]

Straightforward minimization and the application of Stirling’s formula reveal that the summand with the smallest magnitude comes at order $n\approx 1/\gamma\alpha$ and its value $\exp(-1/\gamma\alpha)$ is a proxy for the limited accuracy of the asymptotic expansion. If the coefficients $c_{n}$ have an alternating sign, the series (1) can be Borel resummed and this way a unique value can be assigned to the observable $E$. For sign coherent series, on the other hand, Borel resummation comes with an (imaginary) ambiguity proportional (in leading order) to $\exp(-1/\gamma\alpha)$. In practical calculations this is often not an issue, thanks to the smallness of the expansion parameter $\alpha$. However, in asymptotically free systems — such as four-dimensional non-Abelian gauge and two-dimensional nonlinear sigma models — the long-range regime is always strongly coupled and the ambiguity can severely limit the ability to make physically meaningful statements about $E$. These effects go under the name of renormalons Beneke (1999). Operator product expansion (OPE) offers another view on them as it involves nonperturbative condensates and strong scales with a similar dependence on $\alpha$ (see, e.g., Beneke (1999) and Eq. (7) below).

From experience with various models, the emergence of sign-coherent asymptotic expansions is connected to vacuum degeneracy, as labeled by topological quantum numbers. What makes this subject so fascinating is that the nonperturbative nature of the perturbative ambiguity seems to be connected to the system’s nonperturbative classical tunneling solutions with its typical factors of $\exp(-1/\alpha)$. This concerns instantons as stable topological configurations and their superpositions, but also unstable saddles, as is the case in the models at hand Cherman et al. (2015). There is mounting evidence from resurgence theory calculations that, taken together, information from asymptotic perturbative series and from topology leads to a cancellation of nonperturbative ambiguities and could potentially be used to construct a nonperturbative continuum formulation of QFTs Dunne and Ünsal (2015a, b); Cherman et al. (2014); Dunne and Ünsal (2013, 2012).

So far, a rigorous mathematical proof of the resurgence picture does not exist for most theories of interest and the calculations rely on certain assumptions (or special features of supersymmetric theories). While these assumptions are physically and mathematically well motivated, it is nevertheless important to check their validity. Moreover, most studies of the renormalon mechanism are based on the summation of a special class of Feynman diagrams. The only known ab initio determination of high-order perturbative coefficients is provided through numerical simulations on space-time lattices. The method of choice is not the common Monte Carlo framework, rather stochastic quantization (Langevin dynamics) Parisi and Wu (1981) and the numerical version thereof, combined with numerical perturbation theory to numerical stochastic perturbation theory (NSPT) Di Renzo et al. (1994); Di Renzo and Scorzato (2004). It has the great advantage that its effort grows only quadratically in the expansion order, not factorially as in diagrammatic perturbation theory. With this tool renormalons in four-dimensional $SU(3)$ Yang-Mills theory have been clearly demonstrated in two observables: the lattice action (plaquette) Bauer et al. (2012); Bali et al. (2014) and the energy of static sources (Polyakov loop) Bali et al. (2013). High expansion orders (up to 35 in $1/\beta$), extrapolations in volume and Langevin time, and other sophisticated lattice methods had to be used to improve on previous studies that could not find renormalons. Very recently, quarks have been included in that framework Del Debbio et al. (2018).

Renormalons and resurgence in sigma models have been investigated before; see for example Volin (2010); Dunne and Ünsal (2012, 2013); Cherman et al. (2014, 2015); Fujimori et al. (2016, 2019); Ishikawa et al. (2019); Marino and Reis (2019). There are however no NSPT studies of renormalons for theories other than quantum chromodynamics. To provide cross-checks of the universality of the renormalon picture and resurgence it is important to have reliable first-principle results for a large variety of different theories. In this paper, we present first numerical results for the perturbative coefficients of the energy density of 1+1 dimensional principal chiral models $PC(N)$, with special emphasis on the universality in the rank $N$. In these models, the degrees of freedom are $SU(N)$ group valued fields, and the Euclidean action is nothing but the obvious kinetic term $\frac{1}{g^{2}}\int\!d^{2}x\,\text{tr}\,\partial_{\mu}U(x)\partial_{\mu}U^{\dagger}(x)$. Superficially this looks like the action of a free theory, but we emphasize that the constraint $U\in SU(N)$ introduces couplings between the field components and makes the $PC(N)$ models highly nontrivial. On lattices with spacing $a$ the derivatives translate into nearest neigbor interactions, and the lattice action reads

[TABLE]

Like their $O(N)$ 111Since $SU(2)$ is a three-sphere, there is an isomorphism between the models $PC(2)$ and $O(4)$. and $CP(N-1)$ cousins, these sigma models are asymptotically free 222Correspondingly, the continuum limit $a\to 0$ is achieved by $\beta\to\infty$. and generate a mass and strong scale through quantum fluctuations.

From a statistical physics analogy, the energy density $E$ is related to the $\beta$-derivative of the partition function, in our convention

[TABLE]

where $V$ stands for the number of lattice sites. Note that the energy density has mass dimension $d=2$.

As $\beta N$ is related to the inverse of $g^{2}$ [see Eq. (2)], both the fields and the observable are expanded in powers of $\beta^{-1/2}$,

[TABLE]

where we have immediately used that $E$ only contains integer powers of $1/\beta$ 333Note that in our approach the energy is formally expanded in half-integer powers, just like the fields, but every second of these coefficients is found to vanish, i.e., to be consistent with zero.. Actually, the first few terms of this weak coupling expansion are known analytically Brihaye and Rossi (1984) and will be used as benchmarks for our numerical results. For more details on the expansion within NSPT, we refer to Appendix A.

To develop an expectation for the renormalon behavior of the energy expansion, we first of all notice that, although their homotopy groups are trivial, $PC(N)$ models contain nonperturbative saddles Uhlenbeck (1989) (unitons), which may cure the ambiguity of a sign coherent perturbative expansion Cherman et al. (2014, 2015). Second, we invoke from the large-$N$ expansion 444Where the ’t Hooft coupling kept fixed is $g^{2}N$ and thus $1/\beta$. the two-loop relation between the lattice spacing (inverse cutoff scale), the generated strong scale $\Lambda_{L}$ (proportional to the mass), and the bare lattice coupling $\beta$,555The coefficients $\beta_{0,1}$ parametrize the beta function of the coupling $\alpha=\beta^{-1}$ with the cutoff: $a^{-1}\partial_{a^{-1}}\alpha=-\beta_{0}\alpha^{2}-\beta_{1}\alpha^{3}+\mathcal{O}(\alpha^{4})$; see e.g., Rossi and Vicari (1994); Gonzalez-Arroyo and Okawa (2018)

[TABLE]

The OPE relates the constant $\gamma$ in Eq. (1) to the running of the coupling $\beta=1/\alpha$ and the energy dimension of the observable $E$: since our observable has mass dimension two, its perturbative part should, in leading order, receive nonperturbative corrections of the form $\exp(-2\beta/\beta_{0})$. We therefore anticipate the energy coefficients to behave like

[TABLE]

That the energy possesses factorially growing sign coherent perturbative coefficients has first been demonstrated in Fateev et al. (1994). Applying the general arguments from below Eq. (1) to the expansion coefficients $E_{n}$ themselves (e.g., by formally setting $\alpha=1$) using $\gamma=\beta_{0}/2$, the $E_{n}$ are expected to start growing around order

[TABLE]

Their ratios divided by the order,

[TABLE]

should approach a constant. Note that these leading-order statements should hold independently of the rank $N$. Three-loop corrections are related to the (regularization dependent) beta function coefficient $\beta_{2}$, which contains order $\mathcal{O}(1/N^{2})$ terms in the lattice scheme Rossi and Vicari (1994).

We apply NSPT to calculate the expansion coefficients $E_{n}$ of the energy density on symmetric two-dimensional lattices with the same number $L$ of sites in the spatial and Euclidean time direction, $V=L\times L$, and periodic boundary conditions in both of them. To study finite size effects and the assumed universality of Eq. (12), we consider different lattice geometries and a variety of different ranks $N$. In particular, we calculate the expansion coefficients up to $E_{10}$ for $N=12$ and even up to $E_{20}$ for $N=3,4,5,6$; see Appendix B for details. The simulations for different $N$ are completely independent.

In Langevin simulations the finite stochastic time step $\epsilon$ introduces a systematic error and it is necessary to perform the extrapolation $\epsilon\to 0$. For the numerical integration of the Langevin equation we utilize the Runge-Kutta algorithm from Bali et al. (2013), which is exact up to terms proportional to $\epsilon^{2}$. We use up to five different values of $\epsilon$ in our simulations and perform the $\epsilon\to 0$ limit by fitting a function linear in $\epsilon^{2}$ to our data.

A remarkable feature of NSPT simulations is that neither the lattice spacing $a$ nor the coupling $\beta$ enters the calculations explicitly. All computations are done directly with the expansion coefficients. Therefore, it is not possible to assign a physical volume to our lattice and it is not straightforward to go to the infinite volume limit. Indeed, from OPE arguments very large finite size effects are expected, even on the largest lattices that are achievable in present day simulations Bali et al. (2013, 2014). Fortunately, the OPE enables us to infer the functional form of the finite volume dependence, which makes it possible to extrapolate our results to infinite volume.

For the results presented in this work we used the following equation to take finite volume effects into account:

[TABLE]

where the coefficients $F_{k}$ are polynomials of order $k$ in $\ln(L)$. The details of the infinite volume extrapolation are somewhat lengthy, and we postpone a more detailed discussion of the extrapolation and its systematics to Appendix C.

Figure 1 shows a synopsis of our numerical coefficients $E_{n}$ for the $PC(6)$ model, i.e., a fixed $N$. They perfectly meet the analytically known formulas for $E_{1,2,3}$ and have been determined very precisely over many orders of magnitude. This figure also demonstrates the importance of performing the extrapolation to $\epsilon=0$ and to infinite volume: for high orders the corrections to the expansion coefficients are extremely large and the finite volume results are off by orders of magnitude.

It is interesting to note that the coefficients $E_{n}$ seem to fall off exponentially with $n$. Up to the expansion orders we consider, the asymptotic nature of the expansion is completely hidden by the large value of $n^{*}\approx 50$, only after which the coefficients grow factorially; see Eq. (11). With our current numerical setup, it is not feasible to calculate expansion coefficients up to such high orders to directly observe this growth of the $E_{n}$. This is at variance with four-dimensional gauge theories, where $\beta_{0}=11$ and for the plaquette one computes a much smaller $n^{*}\approx 4/11$, such that those coefficients grow from the start Bali et al. (2014).

However, the derivation of $n^{*}$ assumes that the asymptotic behavior of Eq. (10) of the expansion coefficients sets in before the coefficients start to grow; thus, the ratios $r_{n}$ of consecutive $E_{n}$ should approach the constant in Eq. (12) already for $n<n^{*}$. These ratios are therefore a more sensitive signal for renormalons.

To calculate the ratios $r_{n}$, we first perform the extrapolation to $\epsilon=0$ and then use Eq. (13) to take finite size effects into account. The results are plotted in Fig. 2. We find that the asymptotic behavior sets in somewhere around order $n=15$ independent of the rank $N$ in accordance with the renormalon picture.

While our data clearly show the expected asymptotic behavior, the error bars on the ratios are relatively large. It is tempting to look for a plateau in the ratios $r_{n}$ to perform a fit to a constant and get a better estimate for the asymptotic behavior of the expansion coefficients. The problem with this approach is that subsequent ratios $r_{n}$ and $r_{n+1}$ are strongly correlated, since the coefficient $E_{n}$ is used in the calculation of both of them. The simulations and fits for different rank $N$ are, however, independent. In Fig. 3 we show the results of fitting a constant to $r_{n}(N)$ for fixed $n$. The fit essentially averages over $N$ while also taking the individual errors into account. The results for such a $N$-average with Gaussian error propagation are shown for comparison. For large expansion orders $n\gtrsim 15$, after the asymptotic behavior sets in, the ratios $r_{n}$ should no longer depend on $N$ and combining data for different ranks is justified. The plots in Fig. 3 are in very good agreement with the prediction from Eq. (12) with relative errors that are smaller than $10\%$.

To summarize, we have determined the perturbative coefficients of the energy in particular two-dimensional sigma models by a suitable lattice technique – implicitly summing up factorially many diagrams – including dedicated continuum and infinite volume extrapolations. Using the first few analytically known coefficients as benchmarks, we have determined up to 20 high-order coefficients which spread over many orders of magnitude. Their ratios (divided by the order) clearly approach a constant consistent with the renormalon picture. The latter is based mainly on the coupling dependence of the strong scale (and the mass dimension of the quantity), thus reflecting nonperturbative physics. Our ab initio results validate the universality of the renormalon/resurgence picture, i.e., the independence on the rank of the sigma model. The method utilized in this study can easily be extended to twisted boundary conditions in one or even two directions, which ought to reduce the volume dependence. Likewise, asymptotic series in other observables and/or other sigma models can be studied to shed more light on the resurgence picture.

The NSPT code used in this work Puhr and Simeth (2019) as well as the original data and the data analysis scripts Puhr and Bruckmann (2019) are available on GitHub.

Acknowledgments – We thank Gerald Dunne, Mithat Ünsal and in particular Gunnar Bali for useful discussions, and Marcos Mariño for pointing us to reference Volin (2010). We are indebted to Jakob Simeth for his contributions to the NSPT code used for our simulations. The computations were performed on a local cluster at the University of Regensburg and on the Linux cluster of the LRZ in Garching. GNU parallel Tange (2018) was used to speed up the data analysis. We acknowledge support from DFG (Contracts No. BR 2872/6-2 and No. BR 2872/8-1).

Appendix A NSPT for the Principal Chiral Model

The framework of NSPT applied in this work is not new, but not well-known outside the lattice community. For the convenience of the reader, we therefore give a brief introduction to NSPT in this section. A more thorough exposition can be found in the referenced original works. NSPT was first developed in the context of QCD in Di Renzo et al. (1995); Di Renzo et al. (1994) and a review can be found in Di Renzo and Scorzato (2004).

A.0.1 Numerical Perturbation Theory

The idea of numerical perturbation theory is to formally perform a weak coupling expansion of the lattice fields $U$ in powers of $\beta^{-\frac{1}{2}}$ up to order $\beta^{-M}$,

[TABLE]

where $\beta$ is the lattice coupling. In this work we only consider expansions around the vacuum, where the leading term $U_{0}$ is given by the unit matrix on all lattice sites.

Algebraic operations with these truncated series are straightforward,

[TABLE]

Once addition and multiplication are defined, any analytic function of the fields $U$ can be evaluated by inserting the field expansion into the power series of the function.

The most expensive operation is the multiplication of fields, which requires $\mathcal{O}(M^{2})$ multiplications of coefficients 666This is only true for a naive implementation of the convolution (16). Using, e.g., a fast Fourier transform algorithm for the convolution would require only $\mathcal{O}(M\log_{2}M)$ multiplications.. The numerical cost of an NSTP simulation is consequently roughly proportional to $M^{2}$. Compared to diagrammatic perturbation theory, where the number of Feynman diagrams to be taken into account typically grows like $\mathcal{O}(M!)$, this is very efficient.

A.0.2 Stochastic Quantization on the Lattice

An obvious precondition for the use of numerical perturbation theory is that all the functions involved in the computations can be expanded in powers of $\beta^{-\frac{1}{2}}$. Most state of the art lattice simulations use Monte Carlo methods based on the Metropolis-Hastings algorithm to sample the configuration space. Metropolis-Hastings has a big disadvantage from the point of view of numerical perturbation theory: it employs an accept-reject step, where a proposed configuration update is accepted with a certain probability depending on the change in the action. An accept-reject step cannot be formulated as a function of powers of $\beta^{-\frac{1}{2}}$ and is therefore not suitable for numerical perturbation theory.

One alternative to Metropolis-Hastings type algorithms is stochastic quantization Damgaard and Hüffel (1987), which is based on the Langevin equation. Stochastic quantization introduces a new dimension, usually called the stochastic time $\tau$. The evolution of the fields in $\tau$ is governed by the Langevin equation. For a scalar field $\phi$ with action $\mathcal{S}[\phi]$, the Langevin equation is given by

[TABLE]

where $\bar{\eta}(x,\tau)$ is a Gaussian noise term with the properties

[TABLE]

In the equations above, $\left\langle\cdots\right\rangle_{\bar{\eta}}$ denotes the average over the noise $\bar{\eta}$.

For numerical calculations, the partial differential equation (17) has to be replaced by a finite-difference equation. Making $\tau$ discrete with a step size $d\tau=\epsilon$ leads to

[TABLE]

where $-F_{x}[\phi,\eta]$ is a discretization of the $\tau-$integral of the right-hand side of Eq. (17). A choice with an error of order $\epsilon$ is the Euler method, 777It seems counterintuitive that the Euler method with a term proportional to $\sqrt{\epsilon}$ gives an error of order $\epsilon$. The origin of the square root term is the discretization of Eq. (19). Note that the dimension of the Dirac $\delta$-function $[\delta(\tau-\tau^{\prime})]=[\frac{1}{\tau}]=[\frac{1}{\epsilon}]$, whereas the Kronecker-delta in Eq. (22) is dimensionless. The correct discretization of the $\delta$-function is therefore $\delta(\tau-\tau^{\prime})\to\frac{1}{\epsilon}\delta_{\tau_{i},\tau_{j}}$. If we want to keep the variance of the noise (22) independent of the step size $\epsilon$, we have to multiply every occurrence of $\eta$ by $\epsilon^{-1/2}$, i.e., $\bar{\eta}\to\frac{\eta}{\sqrt{\epsilon}}$. In the Euler discretization this gives rise to the term $\epsilon\frac{\eta}{\sqrt{\epsilon}}=\sqrt{\epsilon}\,\eta$.

[TABLE]

with the Gaussian noise $\eta$, which has zero mean and obeys the discretized version of Eq. (19)

[TABLE]

It can be shown that for sufficiently large $\tau_{n}$ the field configurations $\{\phi(x,\tau_{n})\}$ produced by the discretized Langevin with $F_{x}$ given by Eq. (21) are distributed with probability density $\exp(-\bar{\mathcal{S}}[\phi])$ Batrouni et al. (1985), where the equilibrium action $\bar{\mathcal{S}}$ differs from the continuum action $\mathcal{S}$ by terms of $\mathcal{O}(\epsilon)$,

[TABLE]

Evidently, discrete Langevin simulations suffer from a systematic error, since they sample configurations with respect to a probability distribution that differs from the desired one. This is a disadvantage of the Langevin ansatz and makes it necessary to run simulations for several different values of $\epsilon$ and to extrapolate the results to $\epsilon=0$. On the other hand, the big advantage of stochastic quantization from the point of view of numerical perturbation theory is the absence of an accept-reject step. Stochastic quantization in combination with numerical perturbation theory goes by the name of NSPT and can be used to numerically compute expansion coefficients on the lattice.

For completeness, we mention a recently developed new formulation of NSPT, where the Langevin equation is replaced by an stochastic molecular dynamics (SMD) algorithm Dalla Brida et al. (2017); Dalla Brida and Lüscher (2017). The SMD based NSPT has potential advantages over the Langevin formulation, like smaller autocorrelation times. In this work, we stick to the established Langevin method, which has proven itself in practice in large scale NSPT simulations in the context of QCD.

A.0.3 Discrete Langevin for Constrained Systems

So far, we have only discussed the simple case of an unconstrained scalar field $\phi$. The generalization to fields which are elements of a Lie group is straightforward Batrouni et al. (1985); Catterall et al. (1991). For $U\in SU(N)$, the discretized Langevin equation reads

[TABLE]

where $T^{a}$ stands for the generators of $SU(N)$ (normalized to $\operatorname{Tr}(T^{a}T^{b})=\delta_{a,b}/2$) and the generalization of the Euler term is given by

[TABLE]

with the noise $\eta^{a}=\eta(a,x_{i},\tau_{i})T^{a}$ and the Lie derivative $\nabla_{x}^{a}$ defined via

[TABLE]

Note that $\nabla_{x}^{a}$ acts only on the field $U(x)$ at site $x$. The noise $\eta(a,x_{i},\tau_{i})$ is Gaussian with zero mean and variance,

[TABLE]

For the action (2) of the $PC(N)$ model, the Lie derivative is given by

[TABLE]

Exploiting the property

[TABLE]

of the $SU(N)$ generators, the sum over Lie derivatives in the exponent of Eq. (24) can be performed analytically,

[TABLE]

where M_{\mu}(x):=U(x)U^{\dagger}(x+\mu)-\frac{1}{N}\operatorname{Tr}\big{(}U(x)U^{\dagger}(x+\mu)\big{)}. Note that $M_{\mu}(x)-M_{\mu}^{\dagger}(x)$ is an anti-Hermitian traceless matrix and consequently \exp\big{(}i\sum_{a}\nabla_{x}^{a}\mathcal{S}[U]T^{a}\big{)}\in SU(N), as it should be to ensure that the updated field is also an element of the group.

The use of the Euler integrator (25) again leads to systematic errors of order $\mathcal{O}(\epsilon)$. Runge-Kutta type algorithms for discrete Langevin updates of $SU(N)$-valued fields have been discussed in the literature Batrouni et al. (1985); Torrero and Bali (2008); Bali et al. (2013). In this work we employ the algorithm from Bali et al. (2013), which is of order $\mathcal{O}(\epsilon^{2})$. We find that the higher numerical cost in comparison to the Euler method is offset by the possibility to use larger values of $\epsilon$ in the simulations.

Having defined a discrete Langevin equation for fields $U$, it seems straightforward to apply numerical perturbation theory to this framework. There is, however, a subtle issue. Plugging the expansion (14) into Eq. (28), it becomes clear that the force starts with terms of order $\beta^{\frac{1}{2}}$, whereas the noise term in Eq. (21) is of order $\mathcal{O}(1)$ and the fields are expansions in $\beta^{-\frac{1}{2}}$ 888We use the Euler method as an example here, but the situation is the same for Runge-Kutta algorithms.. Evidently, this leads to inconsistencies. The solution is to redefine the time step $\epsilon\to\epsilon^{\prime}:=\beta\epsilon$. By absorbing $\beta$ in the time step, the force term starts with order $\beta^{-\frac{1}{2}}$. Note that the noise is now also of order $\beta^{-\frac{1}{2}}$ and that the leading-order term of the fields is not changed during a Langevin update. Since the leading-order term is the point around which the fields are expanded — in our case the vacuum — this is what we should expect.

Appendix B Simulation Details

B.0.1 Lattice Setups

For our simulations, we use only symmetric lattices with $V=L\times L$ and a lattice spacing $a$ that is the same in all directions. All simulations are performed with at least three different values of $\epsilon$ for every lattice geometry considered; see Table 1 for details about the lattice setups. As discussed in the previous section, the noise enters the NSPT Langevin equation only at order $\beta^{-\frac{1}{2}}$. It takes several Langevin sweeps before the noise can affect the higher order terms of the field expansion. Moreover, from Eqs. (15) and (16), it is obvious that higher order coefficients have no influence on lower order terms in NSPT calculations. As long as the field coefficients of a given order $l$ are not thermalized, it makes therefore no sense to update terms of order $m>l$. The initial condition of our simulations corresponds to a “cold start”, where all fields are equal to $\mathbbm{1}$. To initialize the higher orders, we start with an expansion up to order two (in $\beta^{-\frac{1}{2}}$) and perform $n_{\text{init}}=500$ Langevin sweeps. Then we increase the expansion order by one and run the simulation for another $n_{\text{init}}$ time steps. This is repeated until we have reached the desired expansion order $M$, after which we perform a final $n_{\text{init}}+1/\epsilon$ initialization sweeps.

We emphasize that our lattice fields only fulfill the constraint $U\in SU(N)$ up to the current expansion order. When we go to a higher order, the unitarization is no longer valid and we have to unitarize the fields up to the new order. It is not straightforward to enforce $U\in SU(N)$, because it is in general not easy to see how the constraints $\det(U)=1$ and $U^{\dagger}U=\mathbbm{1}$ translate to conditions for the expansion coefficients $U_{k}$. A well-known trick is to work with fields $A\in\mathfrak{su}(N)$ from the algebra instead. Any element of $SU(N)$ can be written as

[TABLE]

for a suitable $A$. For $U\in SU(N)$ the field $A=\sum A_{k}\beta^{-\frac{k}{2}}$ has to fulfill $A=A^{\dagger}$ and $\operatorname{Tr}(A)=0$, which can be achieved by making every $A_{k}$ individually Hermitian and traceless,

[TABLE]

To go back and forth between the fields $U$ and $A$ one utilizes the series expansion of $\ln(1+x)$ and $\exp(x)$, respectively. Taking advantage of the expansion around the vacuum, which fixes $U_{0}=\mathbbm{1}$ and $A_{0}=0$, considerably simplifies the calculations. By taking the route over the auxiliary fields $A$, normalizing $U$ becomes straightforward, albeit numerically expensive.

Once the lattice is initialized, we perform $n_{\text{up}}=500/\epsilon$ Langevin sweeps with measurements every $n_{\text{m}}=100$ sweeps. Note that we chose the same stochastic time duration for all our simulations, since we expect the autocorrelation to depend on the stochastic time, not the number of discrete Langevin time steps.

B.0.2 Constraint Violation in Langevin Runs

In practical calculations, it is important to keep in mind that the NSPT algorithm described in the last section respects the constraint $U\in SU(N)$ only if expressed in infinite precision arithmetic. Round-off and numerical errors can lead to a violation of unitarity.

During our runs, we periodically check if the fields still fulfill the constraints. Calculating $A$ every time would be very expensive and considerably slow down the simulations. We therefore check, order by order in the expansion, if the fields $U$ are still unitary. To this end, we calculate

[TABLE]

with $V_{k}$ the expansion coefficients of unitarity,

[TABLE]

where $\|\cdot\|_{F}$ stands for the Frobenius matrix norm. The expectation value $\left\langle\Delta_{k}\right\rangle$ (over all lattice sites) is taken as an estimate of the violation of the constraint $U\in SU(N)$.

Monitoring the Langevin histories shows that the initialization is in general not sufficient to thermalize the highest order coefficients of the energy density. We do not observe any significant violation of unitarity for the lattices with $V\geq 32\times 32$, but the small lattices with $V\lesssim 20\times 20$ show large $\left\langle\Delta_{k}\right\rangle$ even for moderate Langevin times $\tau$. The deviation from unitarity strongly depends on the lattice volume and on the rank $N$. In general, all other parameters being equal, $\left\langle\Delta_{k}\right\rangle$ is larger for higher rank $N$ and higher expansion order $k$.

For the final computation of the expectation values $\left\langle E_{n}\right\rangle$, we make sure to choose only data from a stochastic time interval where the simulation is equilibrated and unitarity violations are small. An example for $N=6$ is shown in Fig. 4. Remember that the $E_{n}$ are coefficients for an expansion in $\beta^{-1}$, whereas the index on the $\Delta_{k}$ corresponds to an expansion in $\beta^{-\frac{1}{2}}$.

To be consistent, for a given lattice size, we choose the same $\tau$ interval for all orders and all parameter sets. In Table 2 we summarize which data we used in the final fits. For convenience, we simply call this data set setup I from now on.

The choice of which data to include is somewhat arbitrary and one should verify that it does not influence the final results. As a cross-check, we consider two additional datasets. The first consists of all the setups in Table 1 and for the second one we use the configurations of Table 2. In both cases, we take data from the stochastic time interval $\tau\in[250,500]$. For later use, we will refer to these as setup II and setup III, respectively. Remarkably, as will be discussed in more detail in the next section, within uncertainties the results for the ratios $r_{n}$ for all three datasets overlap. The reason is probably that for large lattices with $V\gtrsim 24\times 24$, which seem to have a higher weight in the infinite volume extrapolation, unitarity violations are very small.

B.0.3 Extrapolation $\epsilon\to 0$

The Runge-Kutta integrator employed in the simulations is of second order. The extrapolation to $\epsilon=0$ can therefore be performed by fitting a function linear in $\epsilon^{2}$ to our data. The systematics of neglecting higher order terms in the fit function are estimated to be

[TABLE]

where $\epsilon^{\star}$ is the largest $\epsilon$ value used for the given lattice setup which is not larger than the median of the $\epsilon$ values available. Results for a run with $N=3$ and $L=32$ are shown in Fig. 5. In our estimates for the statistical uncertainties of the expansion coefficients we take the autocorrelation into account. The integrated autocorrelation time is estimated using the automatic windowing algorithm from (Madras and Sokal, 1988, Appendix C) with $c=5$. Where available, we compare the NSPT coefficients with results from analytic perturbation theory and find good agreement. The final uncertainty for the expansion coefficients extrapolated to $\epsilon=0$ is obtained by adding the uncertainties from the fit and the systematics in quadrature.

Appendix C Finite volume effects

The extrapolation to infinite volume is based on OPE methods and we closely follow reference Bali et al. (2014). The energy density is symmetric under an exchange of $a\leftrightarrow-a$ and the same is also true for our lattice action. Additionally, we only consider symmetric lattices with equal extend in all directions. It follows that the finite volume effects can only depend on even powers of $L$ and can be parametrized as

[TABLE]

where $E_{n}^{\infty}$ is the infinite volume result.

C.0.1 OPE of the Energy Density

The lattice regularization introduces two distinct scales: the inverse lattice spacing $1/a$ and the inverse of the linear lattice extent $1/(La)$. In the infinite volume limit $aL\gg a$, the two scales are well separated and it makes sense to apply an OPE. The Wilson coefficients then depend on the hard modes of scale $1/a$, whereas the soft modes of scale $1/(La)$ can be described in terms of expectation values of local operators. Formally, we can then expand the expectation value of the perturbative energy density as

[TABLE]

where

[TABLE]

with the infinite volume expansion coefficients $E_{n}^{\infty}$. The expectation value $\left\langle O_{2}\right\rangle$ has to be proportional to $1/(La)^{2}$ on dimensional grounds and we write it as

[TABLE]

Finally, absorbing a constant factor into the definition of the $f_{k}$, the Wilson coefficient is

[TABLE]

With the help of the equations above, we can write $\left\langle E_{\text{pert}}\right\rangle(L)$ in two different ways,

[TABLE]

and

[TABLE]

Note that we are using a different index convention than Bali et al. (2014) here.

We are interested in calculating the coefficients $E_{n}^{\infty}$ from the finite volume results $E_{n}(L)$. To this end, we can use the $\beta$-function of the $PC(N)$ model Rossi and Vicari (1994); Gonzalez-Arroyo and Okawa (2018),

[TABLE]

to expand $\alpha^{k}(1/(La))$ around $\alpha^{k}(1/a)$ in Eq. (41). Comparing (41) and (42) order by order in $\alpha^{k}(1/a)$ finally yields expressions for the functions $F_{n}(L)$ in terms of the coefficients $f_{k},c_{i}$ and the $\beta$-function. In any volume, the leading order of the expansion — the expansion point — is fixed and therefore $f_{0}=F_{0}(L)=0$. The first few nontrivial functions read:

[TABLE]

[TABLE]

and

[TABLE]

In general, $F_{k}(L)$ is a polynomial of degree $(k-1)$ in $\ln(L)$ with coefficients depending on $\{f_{l}\}_{l\leq k}$, $\{c_{j}\}_{j\leq(k-1)}$, and on the $\beta$-function (via the coefficients $\beta_{0},\beta_{1},\cdots$).

A simple closed expression to generate $F_{k}(L)$ for given $k$ is not known. We use SymPy Meurer et al. (2017) to explicitly compute the expansion of Eq. (42) up to order $\alpha^{20}$, from which the $F_{k}(L)$ can be read off.

C.0.2 Fitting the Volume Dependence

Knowing the functional form of the finite volume effects, we can perform fits to our data to extract the infinite volume coefficients. In its most generic form, the fit function we use is given by

[TABLE]

where $F_{n}$ depends on the unknown fit parameters $\{f_{k}\}$ and $\{c_{i}\}$.

The $f_{k}$ and $c_{i}$ couple finite size effects for different expansion orders. We follow the ansatz of Bali et al. (2014) and perform a simultaneous fit to all expansion coefficients. If the fields in the NSPT are expanded up to order $M$, the coefficients $c_{M-1}$ and $f_{M}$ solely appear in $F_{M}$ and only as constant terms in the polynomial in $\ln(L)$. Therefore, from the perspective of the fit, one of them is redundant and we set $c_{M-1}=0$. The coefficients up to $E_{3}^{\infty}$ are known from perturbation theory and are used as input values in the fits. All in all, we end up with a fit function with $3M-5$ free parameters: $(M-3)$ from the infinite volume coefficients $E_{n}^{\infty}$, $M$ parameters $f_{k}$, and $(M-2)$ unknown $c_{i}$.

In all our fits we neglect higher order terms in the $\beta$-function and set $\beta_{i>1}=0$. The leading coefficients are set to their known (and regularization independent) values. The systematic error of truncating the $\beta$-function is estimated by performing separate fits where $\beta_{1}$ is also set to zero. In our fits we observe a behavior that has also been noticed in Bali et al. (2014): in the functions $F_{n}$, the terms containing $c_{i}f_{j-i}$ for fixed $j$ and different $i$ are hard to distinguish in the fitting procedure. The reason can be understood by looking at Eq. (42): the running of the terms $c_{i}f_{j-i}\alpha(1/(La))^{i}\alpha(1/a)^{j-i}$ is very similar for fixed $j$, especially if $j$ is small. This introduces strong correlations between such terms in the fit and ultimately leads to large uncertainties in the infinite volume expansion coefficients. Moreover, when we include systematic errors for the extrapolation $\epsilon\to 0$ in the uncertainties for $E_{n}(L)$, the fits yield unrealistically small $\chi^{2}$ values. Our estimate for the systematic errors is probably too conservative and does not put strong constraints on the fit parameters. As a cross-check, we repeat the fits with only the statistical uncertainties included.

The vastly different scales of the $E_{n}$ for different expansion orders make it hard to compare the coefficients directly. Instead, we consider the ratios $r_{n}$ obtained with the infinite volume coefficients $E_{n}^{\infty}$ from the fits. Uncertainties for the ratios are calculated with Gaussian error propagation from the uncertainties in $E_{n}^{\infty}$.

The plots in Fig. 6 clearly show that we get unacceptably large uncertainties if the coefficients $c_{i}$ are included in our fits. More data and in particular data for much larger lattice volumes would be needed to be able to capture the difference between the $f_{k}c_{i}$ terms.

Since we cannot resolve the $c_{i}$ anyhow, we leave them out of the fitting procedure entirely. Setting $c_{i}=0$ drastically reduces the parameters in our fit function, while keeping its general functional form — a polynomial in $\ln(L)$ — intact. This works astonishingly well and the results are plotted in Fig. 7. Again, we find very small $\chi^{2}$ values of order $10^{-1}$ for the fits with the systematic errors included. With the exception of the fits for $N=6$ and $N=12$, even if we ignore the systematic errors the $\chi^{2}$ values are of order $1$. All four fits yield, with overlapping error bars, the same final results for the ratios. Moreover, the independent fits for different $N$ lie on top of each other after the asymptotic behavior sets in at $n\sim 15$.

The results plotted in Figs. 6 and 7 are obtained with the data from setup I. It is essential to check how strong the choice of the $\tau$ interval and the lattice configurations used for the fits influence the final results. In Figs. 8 and 9 we show the results for the ratios obtained with setup II and setup III, respectively. Finally, in Fig. 10 we show the final results after averaging over the rank $N$ by fitting a constant to $r_{n}(N)$ for fixed $n$ for the setups and fit variants we considered. There is very good agreement between all results. We can therefore conclude that our extrapolation to infinite volume is very robust and systematics from neglecting higher order terms in the beta function are negligible. For the results in the main part of this work, we use the fits obtained with finite $\beta_{1}$ and with systematic errors for the $\epsilon\to 0$ extrapolation. We do not include systematic error estimates for the truncation of the $\beta$-function, since setting $\beta_{1}=0$ in the infinite volume extrapolation does not change the results within uncertainties.