Diagrammatic Monte-Carlo for weak-coupling expansion of non-Abelian lattice field theories: large-N U(N)xU(N) principal chiral model

TL;DR

This paper introduces a Diagrammatic Monte-Carlo method for non-Abelian lattice field theories in the large-N limit, enabling direct study of weak-coupling expansions and non-perturbative phenomena like mass gaps.

Contribution

It develops numerical tools for simulating large-N non-Abelian lattice theories using weak-coupling expansions and demonstrates their effectiveness on the principal chiral model.

Findings

Convergence of double series to non-perturbative mass gap.

Good agreement with traditional Monte-Carlo for energy and transition temperature.

Feasibility of applying the method to planar QCD.

Abstract

We develop numerical tools for Diagrammatic Monte-Carlo simulations of non-Abelian lattice field theories in the t'Hooft large-N limit based on the weak-coupling expansion. First we note that the path integral measure of such theories contributes a bare mass term in the effective action which is proportional to the bare coupling constant. This mass term renders the perturbative expansion infrared-finite and allows to study it directly in the large-N and infinite-volume limits using the Diagrammatic Monte-Carlo approach. On the exactly solvable example of a large-N O(N) sigma model in D=2 dimensions we show that this infrared-finite weak-coupling expansion contains, in addition to powers of bare coupling, also powers of its logarithm, reminiscent of re-summed perturbation theory in thermal field theory and resurgent trans-series without exponential terms. We numerically demonstrate the convergence of these double series to the manifestly non-perturbative dynamical mass gap. We then develop a Diagrammatic Monte-Carlo algorithm for sampling planar diagrams in the large-N matrix field theory, and apply it to study this infrared-finite weak-coupling expansion for large-N U(N)xU(N) nonlinear sigma model (principal chiral model) in D=2. We sample up to 12 leading orders of the weak-coupling expansion, which is the practical limit set by the increasingly strong sign problem at high orders. Comparing Diagrammatic Monte-Carlo with conventional Monte-Carlo simulations extrapolated to infinite N, we find a good agreement for the energy density as well as for the critical temperature of the "deconfinement" transition. Finally, we comment on the applicability of our approach to planar QCD at zero and finite density.

Equations

$$⟨\mathcal{O}⟩\left(\lambda \right)=\sum _i{e}^{-S_i/\lambda }\sum _{l=0}^{m_i}\log {\left(\lambda \right)}^l\sum _{p=0}^{+\infty }{c}_{i,l,p}{\lambda }^p,$$

$$\mathcal{Z}=\underset{U\left(N\right)}{\int }d{g}_x\exp \left(-\frac{N}{\lambda }\sum _{<x,y>}2\mathrm{Re}\mathrm{Tr}\left(g_x^{†}{g}_y\right)\right),$$

$$g_x=\frac{1+i\alpha {\varphi }_x}{1-i\alpha {\varphi }_x}=1+2\sum _{k=1}^{+\infty }{\left(i\alpha \right)}^k{\varphi }_x^k,$$

$$\underset{U\left(N\right)}{\int }d{g}_x=\int d{\varphi }_x\det {\left(1+{\alpha }^2{\varphi }_x^2\right)}^{-N}=$$

$$=\int d{\varphi }_x\exp \left(-N\mathrm{Tr}\ln \left(1+{\alpha }^2{\varphi }_x^2\right)\right)=$$

$$=\int d{\varphi }_x\exp \left(-N\left({\alpha }^2\mathrm{Tr}{\varphi }_x^2+\frac{{\alpha }^4}{2}\mathrm{Tr}{\varphi }_x^4+\dots \right)\right),$$

$$S_0={\lambda }^{-1}\sum _{x,y}{\Delta }_{xy}\mathrm{Tr}\left(g_x^{†}{g}_y\right)=$$

$$=4{\lambda }^{-1}\sum _{k,l=1}^{+\infty }{\left(-i\alpha \right)}^k{\left(i\alpha \right)}^l\sum _{x,y}{\Delta }_{xy}\mathrm{Tr}\left({\varphi }_x^k{\varphi }_y^l\right)=$$

$$=4{\lambda }^{-1}\sum _{\begin{array}{c}k,l=1\\ k+l=2n\end{array}}^{+\infty }{\left(-1\right)}^{\frac{k-l}{2}}{\alpha }^{k+l}\sum _{x,y}{\Delta }_{xy}\mathrm{Tr}\left({\varphi }_x^k{\varphi }_y^l\right),$$

$${\Delta }_{x,y}=2D{\delta }_{x,y}-\sum _{\mu =0}^{D-1}{\delta }_{x,y-\hat{\mu }}-\sum _{\mu =0}^{D-1}{\delta }_{x,y+\hat{\mu }},$$

$$S\left[\varphi \right]=\frac{1}{2}\sum _{x,y}\left({\Delta }_{xy}+\frac{\lambda }{4}{\delta }_{xy}\right)\mathrm{Tr}\left({\varphi }_x{\varphi }_y\right)+$$

$$+\sum _{k=2}^{+\infty }{\left(-1\right)}^{k-1}\frac{{\lambda }^k}{k{8}^k}\sum _x\mathrm{Tr}{\varphi }_x^{2k}+$$

$$+\sum _{\begin{array}{c}k,l=1\\ k+l=2n,k+l>2\end{array}}^{+\infty }{\left(-1\right)}^{\frac{k-l}{2}}\frac{4{\lambda }^{\frac{k+l-2}{2}}}{8^{\frac{k+l}{2}}}\sum _{x,y}{\Delta }_{xy}\mathrm{Tr}\left({\varphi }_x^k{\varphi }_y^l\right),$$

$$\mathcal{Z}=\int d{\varphi }_x\exp \left(-NS\left[\varphi \right]\right).$$

$$\frac{\lambda }{8}\sum _x{\varphi }_x^2=\frac{m_0^2}{2}\sum _x{\varphi }_x^2$$

$$S_F\left[\varphi \right]=\frac{1}{2}\sum _{x,y}\left({\Delta }_{xy}+\frac{\lambda }{4}{\delta }_{xy}\right)\mathrm{Tr}\left({\varphi }_x{\varphi }_y\right)\equiv$$

$$\equiv \frac{1}{2}\sum _{x,y}\left({\Delta }_{xy}+m_0^2{\delta }_{xy}\right)\mathrm{Tr}\left({\varphi }_x{\varphi }_y\right)\equiv$$

$$\equiv \frac{1}{2}\sum _{x,y}{\left({\left(G^0\right)}^{-1}\right)}_{xy}\mathrm{Tr}\left({\varphi }_x{\varphi }_y\right)$$

$$S_I\left[{\varphi }_x\right]=\sum _{v=1}^{+\infty }{\left(-\frac{\lambda }{8}\right)}^v{\xi }^v{S}_I^{\left(2v+2\right)}\left[{\varphi }_x\right],$$

$$S_I^{\left(2v+2\right)}\left[{\varphi }_x\right]=\frac{m_0^2}{2\left(v+1\right)}\sum _x\mathrm{Tr}{\varphi }_x^{2v+2}+$$

$$+\frac{1}{2}\sum _{l=1}^{2v+1}{\left(-1\right)}^{l-1}\sum _{x,y}{\Delta }_{xy}\mathrm{Tr}\left({\varphi }_x^{2v+2-l}{\varphi }_y^l\right),$$

$$⟨\mathcal{O}\left[{\varphi }_x\right]⟩=\underset{\xi \to 1}{\lim }\underset{M\to \infty }{\lim }\sum _{m=0}^M{\xi }^m{\left(-\frac{\lambda }{8}\right)}^m{\mathcal{O}}_m\left(m_0\right),$$

$$W_K\sim \int d^2{q}_1\dots d^2{q}_{f-1}$$

$$\frac{V_1\dots V_v}{\left(\Delta \left(Q_1\right)+m_0^2\right)\dots \left(\Delta \left(Q_l\right)+m_0^2\right)},$$

$$\Delta \left(Q\right)=\sum _{\mu =0}^{D-1}4{\sin }^2\left(\frac{Q_{\mu }}{2}\right)$$

$$m^2\sim \exp \left(-\frac{1}{{\beta }_0\lambda }\right),$$

$$\exp \left(-\frac{1}{{\beta }_0\lambda }\right)=\exp \left(-{\beta }_0^{-1}{e}^{-z}\right)=$$

$$=\sum _{k=0}^{+\infty }\frac{{\left(-{\beta }_0\right)}^{-k}}{k!}{e}^{-kz}=\sum _{l=0}^{+\infty }\frac{{\left(-z\right)}^l}{l!}\sum _{k=0}^{+\infty }\frac{{\beta }_0^{-k}{k}^l}{k!}.$$

$$\mathcal{Z}=\int \mathcal{D}{n}_x\exp \left(-\frac{N}{2\lambda }\sum _{x,y,a}{\Delta }_{xy}{n}_{ax}{n}_{ay}\right),$$

$$\lambda I_0\left(m\right)=1,I_0\left(m\right)\equiv \int \frac{d^{D}p}{{\left(2\pi \right)}^D}\frac{1}{\Delta \left(p\right)+m^2}.$$

$$\displaystyle I_{0}\left(m\right)=\left\{\begin{array}[]{ll}\frac{1}{2m}\,\left(1+\frac{m^{2}}{4}\right)^{-1/2},&\hbox{D = 1;}\\ -\frac{1}{4\pi}\,\log\left(\frac{m^{2}}{32}\right)\left(1+O\left(m^{2}\right)\right),&\hbox{D = 2;}\\ A+Bm+\cdots,&\hbox{D = 3 ,}\end{array}\right.$$

$$\displaystyle m^{2}\left(\lambda\right)=\left\{\begin{array}[]{ll}\sqrt{\lambda^{2}+4}-2,&\hbox{D = 1;}\\ 32\exp{\left(-\frac{4\pi}{\lambda}\right)},&\hbox{D = 2;}\\ m=|B^{-1}|\left(A-1/\lambda\right),&\hbox{D = 3.}\end{array}\right.$$

$$n_{0x}=\frac{1-\frac{\lambda }{4}{\varphi }_x^2}{1+\frac{\lambda }{4}{\varphi }_x^2},n_{ix}=\frac{\sqrt{\lambda }{\varphi }_{ix}}{1+\frac{\lambda }{4}{\varphi }_x^2},$$

$$\mathcal{D}{n}_x=\mathcal{D}{\varphi }_x{\left(1+\frac{\lambda }{4}{\varphi }_x^2\right)}^{-N}.$$

$$S\left[\varphi \right]=\frac{1}{2}\sum _{x,y}\left({\Delta }_{xy}+\frac{\lambda }{2}{\delta }_{xy}\right)\left({\varphi }_x\cdot {\varphi }_y\right)+$$

$$+\sum _{\begin{array}{c}k,l=0\\ k+l\ne 0\end{array}}^{+\infty }\frac{{\left(-1\right)}^{k+l}{\lambda }^{k+l}}{2\cdot 4^{k+l}}\sum _{x,y}{\Delta }_{xy}{\left({\varphi }_x^2\right)}^k{\left({\varphi }_y^2\right)}^l\left({\varphi }_x\cdot {\varphi }_y\right)+$$

$$+\sum _{k=2}^{+\infty }\frac{{\left(-1\right)}^{k-1}{\lambda }^k}{4^{k}k}\sum _x{\left({\varphi }_x^2\right)}^k,$$

Full text

I Introduction

An infamous fermionic sign problem in Monte-Carlo simulations based on sampling field configurations in the path integral is currently one of the main obstacles for systematic first-principle studies of the phase diagram and equation of state of dense strongly interacting matter as described by Quantum Chromodynamics (QCD). This problem has motivated the search for alternative simulation algorithms for non-Abelian gauge theories, which would hopefully avoid the sign problem or make it milder.

One of the alternative first-principle simulation strategies is the so-called Diagrammatic Monte-Carlo approach Van Houcke et al. (2010), conventionally abbreviated as DiagMC, which stochastically samples the diagrams of the weak- or strong-coupling expansions, rather than field configurations. DiagMC algorithms, complemented with the “worm” algorithm techniques Prokof’ev and Svistunov (2001), turned out to be very helpful for eliminating sign problem in physical models with relatively simple weak- or strong-coupling expansions. The physical models and theories which were successfully studied using DiagMC range from strongly interacting fermionic systems Burovski et al. (2006); Prokof’ev and Svistunov (2007); van Houcke et al. (2012); Rossi et al. (2017); Rossi (2017); Tupitsyn and Prokof’ev (2017); Ayyar and Chandrasekharan (2017) relevant in condensed matter physics, via scalar field theories Wolff (2009); Korzec et al. (2011); Gattringer and Kloiber (2013); Davody (2013); Pfeffer and Pollet (2017), to $O\left(N\right)$ and $\mathbb{CP}^{N-1}$ nonlinear sigma models Wolff (2010a, b); Bruckmann et al. (2015); Rindlisbacher and de Forcrand (2017) and Abelian gauge theories with fermionic and bosonic matter fields Korzec and Wolff (2010); Gattringer et al. (2015); Delgado et al. (2013), which are more similar to QCD. This selection of references is by no means exhaustive.

This success has motivated various attempts to apply DiagMC to simulations of non-Abelian lattice gauge theories at finite density Miura et al. (2011); Unger and de Forcrand (2011); de Forcrand et al. (2014); Vairinhos and de Forcrand (2014); Unger (2014); Gattringer and Marchis (2016); de Forcrand et al. (2016) or the effective models thereof Fromm et al. (2013); Mercado et al. (2011); Gattringer (2011). Despite significant progress, these attempts have not so far resulted in a first-principle, systematically improvable simulation strategy. At the qualitative level, one of the most successful strategies is based on the strong-coupling expansion of lattice QCD, which captures such non-perturbative features of QCD as confinement, chiral symmetry breaking and the baryon and meson bound states already at the leading order Rossi and Wolff (1984); Wolff (1985). DiagMC simulations of lattice QCD with staggered fermions at infinitely strong coupling are free of the sign problem in the massless limit, and allow one to reproduce the expected qualitative features of the QCD phase diagram Unger and de Forcrand (2011); de Forcrand et al. (2014). Unfortunately, an extension of this approach to the finite values of gauge coupling meets conceptual difficulties, as the next orders of the strong-coupling expansion coming from the bosonic part of the action should be incorporated in the DiagMC algorithm manually de Forcrand et al. (2014), which becomes more and more complicated for higher orders Unger (2014). Ideally, one would also like to sample the strong-coupling expansion in powers of inverse gauge coupling using DiagMC algorithm. Since on finite-volume lattices the strong-coupling expansion is convergent for all values of coupling, simulating sufficiently high orders should allow to access even the weak-coupling scaling region where lattice theory approaches continuum QCD. Recent attempts to extend the DiagMC simulations to finite gauge coupling using the method of auxiliary Hubbard-Stratonovich variables Vairinhos and de Forcrand (2014) or the Abelian color cycles Gattringer and Marchis (2016) are very promising, but still not quite practical. Development of practical tools for generating high orders of strong-coupling expansion is thus currently an important unsolved problem.

However, in the infinite-volume limit strong coupling expansion is known to have only a finite radius of convergence (see e.g. Langelage et al. (2008)). One can thus expect that as one approaches the weak-coupling regime, the convergence of the strong-coupling series will increasingly stronger depend on the lattice size, which might well make the continuum and infinite-volume extrapolation difficult in practice, especially for the infrared-sensitive physical observables such as hadron masses and transport coefficients. Moreover, in the t’Hooft large-$N$ limit one even expects a phase transition separating the strong-coupling and the weak-coupling phases Gross and Witten (1980); Green and Samuel (1981); Campostrini et al. (1995) even on the finite-volume lattice. These properties of the strong-coupling expansion might significantly limit the ability of DiagMC to work directly in the thermodynamic limit, which is one of the most important advantages of DiagMC over conventional Monte-Carlo methods Rossi (2017); Rossi et al. (2017).

These considerations suggest that the extrapolation to the continuum and thermodynamic limit might be easier for DiagMC simulations based on the weak-coupling expansion. However, the weak-coupling expansion for asymptotically free non-Abelian field theories, such as four-dimensional non-Abelian gauge theories and two-dimensional nonlinear sigma models, has two closely related problems, which make a direct application of a DiagMC approach conceptually difficult at first sight. The first problem are infrared divergences due to massless gluon propagators, which cancel only in the final expressions for physical observables. The second problem is the non-sign-alternating factorial growth of the weak-coupling expansion coefficients, which is commonly referred to as the infrared renormalon Beneke (1999) and cannot be dealt with using the Borel resummation techniques. These factorial divergences in the perturbative expansions of physical observables originate in the renormalization-group running of the coupling constant, and have nothing to do with the purely combinatorial factorial growth of the number of Feynman diagrams with diagram order111The recent work Ayyar and Chandrasekharan (2017) suggests that renormalon divergences might be absent in two-dimensional asymptotically free fermionic theories.. In particular, infrared renormalon behaviour persists also in the large-$N$ limit David (1981, 1982); Fateev et al. (1994a, b); Beneke et al. (1998), in which the number of Feynman diagrams grows only exponentially with order Brézin et al. (1978). It is known that if one stops the running of the coupling at some artificially introduced infrared cutoff scale (e.g. by considering the theory in a finite volume), this factorial growth disappears Bali et al. (2014); Bali and Pineda (2016), thus in a sense it is a remainder of IR divergences in perturbative expansion.

Recently infrared renormalon divergences have been given a physical interpretation in terms of the action of unstable or complex-valued saddles of the path integral, along with a resummation prescription based on the mathematical notion of resurgent trans-series Dunne and Unsal (2013); Cherman et al. (2014); Dunne and Unsal (2014a, b, 2015); Behtash et al. (2015). Trans-series is a generalization of the ordinary power series to a more general form

[TABLE]

where $\langle\,\mathcal{O}\,\rangle\left(\lambda\right)$ is the expectation value of some physical observable as a function of the coupling constant $\lambda$, $c_{i,l,p}$ are the coefficients of perturbative expansion of path integral around the $i$-th saddle point of the path integral with the action $S_{i}$ and $l$ labels flat directions of the action in the vicinity of this saddle. It turns out that non-Borel-summable factorial divergences cancel between perturbative expansions around different saddle points $i$ in (1). Unfortunately, currently no method is available to generate the trans-series (1) for strongly coupled field theories in a systematic way. The resurgent analysis of perturbative expansions has been done so far almost exclusively in the regime where gauge theories or sigma models are artificially driven to the weak-coupling limit by a special compactification of space-time to $D=1+0$ dimensions with twisted boundary conditions Dunne and Unsal (2013); Cherman et al. (2014); Dunne and Unsal (2014a, b, 2015); Behtash et al. (2015).

The aim of this work is to develop DiagMC methods suitable for the bosonic sector of non-Abelian lattice field theories, which is currently a challenge for DiagMC approach. We develop two practically independent tools, which are finally combined together to perform practical simulations.

The first tool is the resummation prescription which renders the bare weak-coupling expansion in non-Abelian field theories infrared-finite and thus suitable for DiagMC simulations which sample Feynman diagrams. This prescription is introduced in Section II and tested in Section III on some exactly solvable examples. The basic idea is to work with massive bare propagators, where the bare mass term comes from the nontrivial path integral measure for non-Abelian groups and is proportional to the bare coupling. By analogy with Fujikawa’s approach to axial anomaly Fujikawa (1979), this bare mass serves as a seed for dynamical mass generation and conformal anomaly in non-Abelian field theories which are scale invariant at the classical level. A similar partial resummation of the perturbative series with the help of the bare mass term is also often used in finite-temperature quantum field theory Espinosa et al. (1993); Arnold and Zhai (1994), in particular, in the contexts of hard thermal loop approximation Braaten and Pisarski (1992) and screened perturbation theory Karsch et al. (1997). Upon this resummation, the weak-coupling expansion becomes formally similar to the trans-series (1), but contains only the powers of coupling and of the logs of coupling only, while the exponential factors $e^{-S_{i}/\lambda}$ seem to be absent. This expansion appears to be manifestly infrared-finite and free of factorial divergences at least in the large-$N$ limit. An exactly solvable example of the $O\left(N\right)$ nonlinear sigma model demonstrates that our expansion reproduces non-perturbative quantities such as the dynamically generated mass gap (see Subsection III.2).

The second tool, described in Section IV, is the general method for devising DiagMC algorithms from the full untruncated hierarchy of Schwinger-Dyson equations Buividovich (2011a, b), bypassing any explicit construction of a diagrammatic representation. Here the basic idea is to generate a series expansion by stochastic iterations of Schwinger-Dyson equations. We believe the method to be advantageous for sampling expansions for which the complicated form of diagrammatic rules makes the straightforward application of e.g. worm algorithm hardly possible, with strong-coupling expansion in non-Abelian gauge theories being one particular example. In this work we apply the method to sample the (re-summed) weak-coupling expansion of Section II, which for non-Abelian lattice field theories also has complicated form due to the multitude of higher-order interaction vertices.

As a practical application which takes advantage of both tools, in this work we consider the large-$N$ $U\left(N\right)\times U\left(N\right)$ principal chiral model on the lattice, which is defined by the following partition function:

[TABLE]

where integration in the path integral is performed over unitary $N\times N$ matrices $g_{x}$ living on the sites of the square two-dimensional lattice, $\lambda$ is the t’Hooft coupling constant which is kept fixed when taking the large-$N$ limit, and $\sum\limits_{<x,y>}$ denotes summation over all pairs of neighboring lattice sites $x$ and $y$. Numerical results for this model obtained by combining the methods of Sections II and IV are presented in Section V.

While we expect that our approach can be extended to large-$N$ pure gauge theory without conceptual difficulties, in this proof-of-concept study we prefer to work with the principal chiral model (2). The reason is that while this model exhibits non-perturbative features very similar to those of non-Abelian gauge theories, including dynamical mass gap generation, dimensional transmutation and infrared renormalons Fateev et al. (1994a, b), it has a much simpler structure of the perturbative expansion than non-Abelian gauge theories. Correspondingly, the DiagMC algorithm is also simpler to implement. We discuss the extension of our approach to large-$N$ gauge theories and illustrate the expected structure of infrared-finite weak-coupling expansion for this case in Section VI.

II Infrared-finite weak-coupling expansion for the $U\left(N\right)\times U\left(N\right)$ principal chiral model

The first step in the construction of lattice perturbation theory is to parameterize the small fluctuations of lattice fields around some perturbative vacuum state. For non-Abelian lattice fields $g_{x}\in SU\left(N\right)$ the most popular parametrization is the exponential mapping from the space of traceless Hermitian matrices $\phi_{x}$ to $SU\left(N\right)$ group: $g_{x}=\exp{\left(i\alpha\phi_{x}\right)}$, where $\alpha$ is related to the coupling constant and $g_{x}=1$ is the perturbative vacuum. While all the results in this paper could be also obtained for exponential mapping, for our purposes it will be more convenient to use the Cayley map

[TABLE]

where $\phi_{x}$ are again Hermitian matrices and the value of $\alpha$ will be specified later. Since in this work we will be working in the large-$N$ limit in which the $U\left(N\right)$ and $SU\left(N\right)$ groups are indistinguishable, we omit here the zero trace condition for $\phi_{x}$ and work with $U\left(N\right)$ group. This will significantly simplify the derivation of Schwinger-Dyson equations in terms of $\phi_{x}$ fields, see Subsection IV.1.

One of the advantages of the Cayley map is the particularly simple form of the $U\left(N\right)$ Haar measure (see Appendix B.1 for the derivation):

[TABLE]

where the integral on the right-hand side is over all Hermitian matrices $\phi_{x}$ and $\ldots$ represent the terms of order $O\left(\alpha^{6}\phi_{x}^{6}\right)$ or higher. Thus the Cayley map provides a unique one-to-one mapping between the $U\left(N\right)$ group manifold and the space of all Hermitian matrices. Strictly speaking, one should exclude a single point $g_{x}=-1$ from this mapping, which, however, has zero measure in the path integral. In contrast, for exponential mapping the integration over $\phi_{x}$ should be restricted to a certain region within the space of Hermitian matrices in order to avoid multiple covering of $U\left(N\right)$ group. In addition, for exponential mapping the exponentiated Jacobian contains double-trace terms Buividovich (2015), which make the planar perturbation theory technically somewhat more complicated (see Appendix B.2). The Cayley map was also used in the context of lattice QCD in the Landau gauge von Smekal et al. (2007, 2008); Sternbeck and von Smekal (2010); von Smekal and Bischoff (2012) and for studying unitary matrix models Mizoguchi (2005).

The classical action of the principal chiral model (2) can be rewritten in terms of the fields $\phi_{x}$ as

[TABLE]

where we have introduced the lattice Laplacian operator in $D$-dimensional Euclidean space

[TABLE]

and used the unitarity of $g_{x}$ to add the diagonal terms ${\rm Tr}\,\left(g^{{\dagger}}_{x}g_{x}\right)=N$ to the action. Here $\mu=0,1$ labels the two lattice coordinates, and $\hat{\mu}$ denotes the unit lattice vector in the direction $\mu$. We see now that in order to get the canonical normalization of the kinetic term in the action $S_{0}=\frac{1}{2}\sum\limits_{x,y}\phi_{x}\Delta_{x,y}\phi_{y}+\ldots$, where $\ldots$ denotes the terms with more than two fields $\phi$, we have to set $\alpha^{2}=\frac{\lambda}{8}$. We also note that the classical action (II) is invariant under the shifts of the field variable $\phi_{x}\rightarrow\phi_{x}+cI$, where $I$ is the $N\times N$ identity matrix. This symmetry is a consequence of the scale invariance of the classical action of the principal chiral model.

Rewriting the path integral in the partition function (2) in terms of the new fields $\phi_{x}$, we also need to include the Jacobian (II). Adding the expansion of the Jacobian (II) in powers of $\lambda$ and $\phi$ to the classical action (II), we obtain the full action for the fields $\phi_{x}$:

[TABLE]

so that the partition function (2) reads

[TABLE]

The observation which is central for this work is that the Jacobian (II) contributes a term

[TABLE]

to the quadratic part

[TABLE]

of the action (II), which endows the field $\phi_{x}$ with the bare mass $m_{0}^{2}=\frac{\lambda}{4}$. Thus the effect of the Jacobian (II) is to break the shift symmetry $\phi_{x}\rightarrow\phi_{x}+c\,I$ of the classical action, serving as a seed for dynamical scale generation and conformal anomaly in the principal chiral model (2).

A conventional approach in lattice perturbation theory, which constructs weak-coupling expansion as formal series in powers of coupling $\lambda$, is to treat the quadratic term (9) in the action (II) as a perturbation Brihaye and Rossi (1984) on top of the kinetic term. Perturbative expansion then involves massless propagators $\left(\Delta^{-1}\right)_{xy}$ coming from the scale-invariant classical part of the action, which leads to infrared divergences in the contributions of individual diagrams. These infrared divergences, however, cancel in the final results for physical observables Brihaye and Rossi (1984). Furthermore, the resulting power series in the coupling constant $\lambda$ contain factorially divergent terms, some of which have non-alternating signs and thus cannot be removed by Borel resummation techniques. The only recently proposed way to deal with these divergences is to use the mathematical apparatus of resurgent trans-series Dunne and Unsal (2013); Cherman et al. (2014); Dunne and Unsal (2014a, b, 2015); Behtash et al. (2015).

In this work we consider a different approach to the weak-coupling expansion based on the action (II), somewhat similar in spirit to the screened perturbation theory Espinosa et al. (1993); Arnold and Zhai (1994); Karsch et al. (1997) and to hard thermal loop perturbation theory Braaten and Pisarski (1992); Rebhan (1994). Namely, we include the mass $m_{0}^{2}=\frac{\lambda}{4}$ coming from (9) into the bare lattice propagators $G^{0}_{xy}=\left(\Delta+m_{0}^{2}\right)^{-1}_{xy}$ which are used to construct the weights of Feynman diagrams. This corresponds to a trivial resummation of an infinite number of chain-like diagrams with “two-leg” vertices. Since this prescription puts the coupling $\lambda$ both in the propagators and vertices, we need to define the formal power counting scheme for our expansion. To this end we group together the terms with the products of $\left(2v+2\right)$ fields, $v=1,2,\ldots+\infty$, which correspond to interaction vertices with $\left(2v+2\right)$ legs, and multiply them with the powers $\xi^{v}$ of an auxiliary parameter $\xi$ which should be set to $\xi=1$ to obtain the physical result:

[TABLE]

where in the last equation we have used our definition $m_{0}^{2}=\frac{\lambda}{4}$ in order to formally eliminate the coupling $\lambda$ from interaction vertices $S_{I}^{\left(2v+2\right)}\left[\phi_{x}\right]$. This construction is similar to the “shifted action” of Rossi et al. (2016).

Perturbative expansion of any physical observable $\langle\,\mathcal{O}\left[\phi_{x}\right]\,\rangle$ is then organized as a formal power series expansion in the auxiliary parameter $\xi$ in (11):

[TABLE]

where $M$ is the maximal order of the expansion. This means that we ignore the relation between the mass $m_{0}^{2}=\frac{\lambda}{4}$ and the coupling constant $\lambda$ and treat the terms $S_{I}^{\left(2v+2\right)}\left[\phi_{x}\right]$ as bare vertices with $2v+2$ legs, counting only those powers $\lambda^{v}$ of the t’Hooft coupling constant which come from the vertex pre-factors in (11). The bare mass $m_{0}^{2}$ in the expansion coefficients in (13) is substituted with $\lambda/4$ at the very end of the calculation.

Let us discuss the general structure of an expansion (13) for the free energy of the two-dimensional principal chiral model (2) in the large-$N$ limit, in which only connected planar Feynman diagrams with no external legs contribute. Consider such a planar diagram with $f-1$ independent internal loop momenta $q_{1},\ldots,q_{f-1}$, and with $l$ bare propagators and $v$ vertices. The kinematic weight of such a diagram can be written as

[TABLE]

where $Q_{1}\ldots Q_{l}$ are momenta flowing through each diagram line,

[TABLE]

is the lattice Laplacian (6) in momentum space which behaves as $\Delta\left(Q\right)\sim Q^{2}$ at small momenta and $V_{1}\ldots V_{v}$ are the weights of the vertices which in general depend on $q_{1}\ldots q_{f-1}$. The momenta $Q_{1}\ldots Q_{l}$ are in general not independent and can be expressed as some linear combinations of $q_{1},\ldots,q_{f-1}$. All momenta belong to the Brillouin zone $Q_{\mu}\in\left[-\pi,\pi\right]$ of the square lattice. Thus any possible ultraviolet divergences in the integral (II) are regulated by the lattice ultraviolet cutoff $\Lambda_{UV}\sim 1$, and we only have to care about infrared divergences. Each vertex in (II) contributes either a power of $m_{0}^{2}$ or a square of some combination of momenta $q_{1}\ldots q_{f-1}$, coming from the first or the second terms in (II), respectively. If our planar Feynman diagram is considered as a planar graph drawn on the sphere, the number $f$ is the number of faces of this planar graph. We now take into account the identity $f-l+v=2$ for the Euler characteristic of a planar graph. Applying now the standard dimensional analysis, we see that the above relation between $f$, $l$ and $v$ implies that $W_{K}$ can only contain the terms proportional to $\Lambda_{UV}^{2}$, $m_{0}^{2}$, $m_{0}^{4}/\Lambda_{UV}^{2}$ and so on, probably times logarithmic terms $\sim\log\left(m_{0}^{2}\right)$. In Subsection III.2 we will explicitly illustrate the importance of these logarithmic terms. No negative powers of $m_{0}$ can ever appear.

Remembering now the relation $m_{0}^{2}=\frac{\lambda}{4}$, and taking into account the powers of coupling associated with vertices, we immediately conclude that the weights of Feynman diagrams in our perturbation theory with the bare mass $m_{0}\sim\lambda$ are proportional to positive powers of $\lambda$, times possible logarithmic terms $\log\left(\lambda\right)^{m}$, $\log\left(\log\left(\lambda\right)\right)^{m}$ ($m>0$) and so on. From (II) we see that the combinatorial pre-factor of vertex with $2v+2$ legs grows at most linearly with $v$. Since the number of planar diagrams which contribute to the expansion (13) at order $m$ grows at most exponentially with $m$, and the weight of any diagram is both UV and IR finite, contributions which grow factorially with $m$ cannot appear in the expansion (13) in $D=2$ dimensions.

The numerical results which we present below suggest that (13) is a convergent weak-coupling expansion. At first sight this contradicts the common lore that perturbative expansion cannot capture the non-perturbative physics of the model (2). Let us remember, however, that our expansion is no longer a strict power series in $\lambda$, which is known to be factorially divergent. Rather, the mass term $m_{0}^{2}\sim\lambda$ in bare propagators allows for logarithms of coupling in the series, in close analogy with resummation of infrared divergences in hard thermal loop perturbation theory Rebhan (1994). It is easy to convince oneself that the formal expansion which contains powers of logs of coupling along with the powers of coupling can incorporate the non-perturbative scaling of the dynamical mass gap

[TABLE]

where $\beta_{0}$ is the first coefficient in the perturbative $\beta$-function. Indeed, we can rewrite the function (16) as a convergent expansion in powers of $z=\log\left(\lambda\right)$:

[TABLE]

In what follows we apply this resummation prescription to construct the infrared-finite weak-coupling expansion for several models: exactly solvable large-$N$ $O\left(N\right)$ sigma model in $D=1,2,3$ dimensions, exactly solvable large-$N$ principal chiral model on one-dimensional lattices with $L_{0}=2,\,3,\,4$ lattice sites, and, most importantly, the principal chiral model in $D=2$ dimensions at zero and finite temperature.

III Exactly solvable examples and counterexamples

III.1 Large-$N$ $O\left(N\right)$ sigma model in $D=1,\,2,\,3$ dimensions

Large-$N$ $O\left(N\right)$ sigma model provides one of the simplest examples of a non-perturbative quantum field theory which can be exactly solved by using saddle point approximation. In $D=2$ dimensions this model has dynamically generated non-perturbative mass gap, similarly to the more complicated principal chiral models and non-Abelian gauge theories. On the other hand, for this model one can also explicitly construct the formal weak-coupling perturbative expansion, which is dominated by the “cactus”-type diagrams, schematically shown on Fig. 1. These properties make this model an ideal testing ground for the infrared-finite weak-coupling expansion described in Section II.

The partition function of the $O\left(N\right)$ sigma model is given by the path integral over $N$-component unit vectors $n_{a\,x}$, $a=0\ldots N-1$ attached to the sites of $D$-dimensional square lattice, which we label by $x$:

[TABLE]

where $\Delta_{x,y}$ is again the lattice Laplacian (6). The well-known saddle-point solution shows that the two-point correlation function $G_{xy}=\langle\,n_{a\,x}\,n_{a\,y}\,\rangle$ is proportional to the propagator of a free massive scalar field, with the mass $m$ being the solution of the gap equation:

[TABLE]

In $D=1,\,2,\,3$ dimensions the function $I_{0}\left(m\right)$ is

[TABLE]

where in the last line $A=0.252731$ and $B=-0.0795775$. If the mass gap $m$ is small enough, the corresponding solutions for $m$ read

[TABLE]

The goal of this Section is to construct the infrared-finite weak-coupling expansion following the prescription of Section II. In order to parameterize the fluctuations of the $n_{a\,x}$ field around the classical vacuum with $n_{0\,x}=1$, we use the stereographic projection from the unit sphere in $N$-dimensions to $N-1$-dimensional real space of fields $\phi_{i\,x}$, $i=1,\ldots,N-1$, which is the analogue of the Cayley map (3):

[TABLE]

where $\phi_{x}^{2}\equiv\sum\limits_{i}\phi_{i\,x}\phi_{i\,x}$. This is a unique map between the unit sphere embedded in $N$-dimensional space and the $N-1$-dimensional real space, in which only a single point with $n_{0\,x}=-1$, $n_{i\,x}=0$ is excluded. The integration measure on the unit sphere can be written in terms of the stereographic coordinates $\phi_{i\,x}$ as Buividovich (2015) (see also Appendix B.3):

[TABLE]

Following the construction of Section II, we now include this Jacobian in the action for the fields $\phi_{i\,x}$, which reads:

[TABLE]

where $\phi_{x}\cdot\phi_{y}\equiv\sum\limits_{i}\phi_{i\,x}\cdot\phi_{i\,y}$. We have also taken into account that the terms of the form $\sum\limits_{x,y}\Delta_{xy}\left(\phi_{x}^{2}\right)^{k}\left(\phi_{y}^{2}\right)^{l}$ vanish at large $N$ due to factorization property $\langle\,\left(\phi_{x}\cdot\phi_{y}\right)\,\left(\phi_{z}\cdot\phi_{t}\right)\,\rangle=\langle\,\phi_{x}\cdot\phi_{y}\,\rangle\langle\,\phi_{z}\cdot\phi_{t}\,\rangle$ and translational invariance. In full analogy with the derivation of the action (II) in Section II, the Jacobian of the transformation from compact to non-compact field variables results in the bare “mass term” $m_{0}^{2}=\lambda/2$.

A convenient way to arrive at the perturbative expansion for the action (III.1) is to iterate the corresponding Schwinger-Dyson equations

[TABLE]

which in the large-$N$ limit involve only the two-point function

[TABLE]

The form of equation (III.1) suggests that the solution has the form of the free scalar field propagator with the mass $m$, times some wave function renormalization factor $z$:

[TABLE]

Substituting this ansatz in (III.1), we obtain the following equations for $m$ and $z$:

[TABLE]

where we have again introduced the auxiliary parameter $\xi$, to be set to $\xi=1$, in order to define the formal power counting scheme. Substituting $\xi=1$ and $m_{0}^{2}=\lambda/2$, we immediately arrive at the exact solutions (27) for the mass $m$. For $z$ the exact solution is simply $z=2$.

On the other hand, we can use the equations (III.1) for an iterative calculation of the coefficients of the formal weak-coupling expansion

[TABLE]

which is also the formal power series in the auxiliary parameter $\xi$ in (III.1). In practice, we continue the iterations up to some finite order $M$, and calculate $m^{2}$ and $z$ by summing up all the terms in the series (III.1) with powers of $\xi$ less than or equal to $M$. After such a truncation of series, we substitute $\xi=1$ and $m_{0}\rightarrow\lambda/2$.

In $D=1$ dimensions, our prescriptions yields the series with summands of the form $\lambda^{k}\left(8+\lambda\right)^{-M-2}$ and $\lambda^{k+1/2}\left(8+\lambda\right)^{-M-3/2}$, $k=1,2,\ldots$, times some integer-valued coefficients which depend on $M$. The expansion thus behaves regularly in the limit $\lambda\rightarrow 0$, where the mass goes to zero linearly in $\lambda$.

In $D=2$ dimensions, we get the expansion which involves both powers of $\lambda$ and $\log\left(\lambda\right)$ (see also Buividovich (2015)):

[TABLE]

and similarly for $z\left(\lambda,M\right)$. As discussed above, the expansion involving the logs of $\lambda$ can reproduce the non-perturbative scaling (27) of the mass gap in $D=2$. Finally, in $D=3$ dimensions we obtain the expansion in positive powers of $\sqrt{\lambda}$.

On Fig. 2 we compare the convergence of the expansions in $D=1,\,2,\,3$ dimensions, plotting the relative error

[TABLE]

of the truncated series (III.1) with respect to the exact result $m\left(\lambda\right)$ given by (27). The values of the coupling $\lambda$ are fixed in such a way that the exact non-perturbative mass gap is equal to $m=0.4$ for all dimensions $D=1,\,2,\,3$. The series (III.1) exhibit the fastest convergence at $D=1$. At $D=2$, the convergence is not monotonic, and also slower than in $D=1$. In $D=3$ the convergence is slowest, and from the convergence plot in Fig. 2 it is difficult to say whether the series converge or not, even with up to $500$ terms in the expansion. These results suggest that for the large-$N$ $O\left(N\right)$ sigma model $D=2$ seems to be the critical dimension separating the convergent and non-convergent expansions.

III.2 Non-perturbative mass gap in the 2D $O\left(N\right)$ sigma model

We now turn to the physically most interesting example of the two-dimensional $O\left(N\right)$ sigma model in the large-$N$ limit, where the mass gap (27) exhibits non-perturbative scaling $m^{2}\sim\exp{\left(-\frac{1}{\beta_{0}\lambda}\right)}$ similar to the one found also in two-dimensional principal chiral model and in non-Abelian gauge theories. On Fig. 3 we illustrate the convergence of the series (III.1) for the mass $m$ and the renormalization factor $z$ to the exact results (27) and $z=2$ for different values of the coupling constant $\lambda$. To make the convergence at large truncation orders $M$ more obvious, on Fig. 4 we also show the relative error (37) of the truncated series (III.1) as a function of $M$ in logarithmic scale for several values of $\lambda$. While at large values of $\lambda$ the series converge with the precision of $10^{-16}$ at $M\sim 500$, the convergence becomes slower at smaller values of $\lambda$ - and, correspondingly, at smaller values of the mass $m$.

An interesting feature of the dependence of $m\left(\lambda,M\right)$ on $M$ is that it is not monotonous, but rather has several extrema with respect to $M$. We denote the smallest value of $M$ at which the dependence of $m\left(\lambda;M\right)$ on $M$ has a minimum as $M^{\star}_{1}$, the position of the next maximum of this dependence - as $M^{\star}_{2}$, and so on. These positions depend on $\lambda$ and shift towards larger $M$ at smaller $\lambda$. We can think of the extrema positions $M^{\star}_{1},M^{\star}_{2},\ldots$ as of some characteristic scales which determine the convergence rate of the truncated series. In particular, at $M>M^{\star}_{1}$ the convergence to the exact value is rather fast. Remarkably, we have found that $M^{\star}_{1}$ as a function of $\lambda$ with good precision appears to be proportional to the square of correlation length (in lattice units):

[TABLE]

For illustration, on Fig. 5 we plot $\log\left(M^{\star}_{1,2}\left(\lambda\right)\right)$ versus $1/\lambda$, together with the linear fits of the form $\log\left(M^{\star}\left(\lambda\right)\right)={\rm const}-\beta/\lambda$. For $M^{\star}_{1}$ the best fit yields the coefficient $\beta=12.49$ which is very close to $4\pi=12.57$. For $M^{\star}_{2}$ we get the best fit value $\beta=9.33$, which is close to $3\pi=9.42$ and suggests the scaling $M^{\star}_{2}\sim l_{c}^{3/2}\sim m^{-3/2}$. This might imply that at sufficiently small $\lambda$ the two extrema might merge. However, for $M^{\star}_{2}$ we have fewer data points, and those exhibit some tendency for faster growth at larger values of $1/\lambda$, thus the scaling exponent for $M^{\star}_{2}$ might be different at asymptotically small $\lambda$.

The data presented on Figs. 4 and 5 provide a clear numerical evidence that the infrared-finite weak-coupling expansion (III.1) indeed incorporates the non-perturbative dynamically generated mass gap $m^{2}=32\exp{\left(-\frac{4\pi}{\lambda}\right)}$. If we were to sample the series (III.1) using DiagMC algorithm, the linear scaling of $M^{\star}_{1}\left(\lambda\right)$ with correlation length would also result in the expectable critical slowing down of simulations near the continuum limit.

The series which we obtain from the field-theoretical weak-coupling expansion (III.1) based on Schwinger-Dyson equations (III.1) are in this respect drastically different from the simplest formal expansion (II) of the exact answer (27) in powers of logs. We have explicitly checked that for this case the positions of the extrema of $m$ with respect to $M$ are independent of $\lambda$. In this respect our infrared-finite weak-coupling expansion (III.1) seems to be closer to the notion of trans-series.

Let us stress, however, that the series (III.1) are not trans-series yet, despite formal similarity with (1). In a mathematically strict sense, resurgent trans-series (1) should give a “minimal” representation of a function with optimal convergence. For the large-$N$ $O\left(N\right)$ sigma model, this optimal trans-series representation for the dynamically generated mass is given by a single term $m^{2}=32\exp{\left(-\frac{4\pi}{\lambda}\right)}$ (see the exact solution (27)). Thus the true trans-series for this model consist of only one term, which is enough for convergence! The expansion (III.1) also converges to the exact result (27), but at a non-optimal rate 222We thank Ovidiu Costin for pointing out these properties of trans-series to us at the “Resurgence 2016” workshop in Lisbon..

III.3 Exactly solvable counter-example: finite chiral chain models

The $U\left(N\right)\times U\left(N\right)$ principal chiral model (2) can be solved exactly in $D=1$ dimensions for lattice sizes $L_{0}=2,\,3,\,4$ and $L_{0}=+\infty$ Rossi et al. (1998). The first and the last cases can be reduced to the Gross-Witten unitary matrix model Gross and Witten (1980). For other two cases, solutions can be found in Rossi et al. (1998). For all values of $L_{0}$ the one-dimensional principal chiral model exhibits a third-order Gross-Witten phase transition at some critical coupling $\lambda=\lambda_{c}\left(L_{0}\right)$, which separates the strong-coupling and the weak-coupling regimes. In particular, for the mean link $\langle\,\frac{1}{N}{\rm Tr}\,\left(g^{{\dagger}}_{1}g_{0}\right)\,\rangle$ the exact results in the weak-coupling regime $\lambda<\lambda_{c}$ are:

[TABLE]

For $L_{0}=4$ the weak-coupling solution is only available in the implicit form:

[TABLE]

where $\delta$ is the solution of the equation $\frac{8}{\pi\lambda}\left(E\left(1-\delta^{2}\right)-\delta^{2}F\left(1-\delta^{2}\right)\right)=1$ and $E$ and $F$ are the complete elliptic integrals of second and first type, respectively.

These analytic expressions offer the possibility to check the convergence of the infrared-finite weak-coupling expansion of Section II for small lattice sizes. Of course, for finite lattices the dimensional analysis following the equation (II) would not work, and loop integrals will only produce some rational functions of $\lambda$ instead of logs. Nevertheless it is interesting to check whether these rational functions provide a good approximation to the exact results (39), (III.3) and (41). To this end we have performed exact calculations of the coefficients of the weak-coupling expansion (13), generating them recursively up to some maximal order $M$ from the Schwinger-Dyson equations (55), to be presented in Subsection IV.1. Such an exact calculation is possible for $L_{0}=2$ up to $M=16$ and for $L_{0}=4$ up to $M=8$. Explicit expression for the mean link $\langle\,\frac{1}{N}{\rm Tr}\,\left(g^{{\dagger}}_{1}g_{0}\right)\,\rangle$ in terms of the correlators of $\phi_{x}$ matrices and the prescription for truncating the expansion at some finite order $M$ are given in Subsection V.1.

On Fig. 6 we illustrate the relative error of the truncated weak-coupling expansion (13) for the mean link $\langle\,\frac{1}{N}{\rm Tr}\,\left(g^{{\dagger}}_{1}g_{0}\right)\,\rangle_{M}$ as a function of inverse truncation order $1/M$ for $L_{0}=2,\,3,\,4$, fixing $\lambda/\lambda_{c}=0.5$ and $\lambda/\lambda_{c}=0.8$. The relative error is defined in full analogy with (37). At the first sight the convergence seems to be rather fast, however, a detailed inspection of numerical results reveals a tiny discrepancy of order of $0.1\%$ which cannot be ascribed to numerical round-off errors, and which seems to stabilize at some finite value as $M$ is increased (at least for $L_{0}=2$). An analytic calculation of the few first orders revealed that the weak-coupling expansion (13) correctly reproduces only the three leading coefficients of the standard perturbative expansion in powers of $\lambda$, which is convergent for finite $L_{0}$.

It seems thus that for finite-size one-dimensional lattices we observe the misleading convergence of the re-summed perturbative expansion to some unphysical branch of the solutions of Schwinger-Dyson equations. This phenomenon might happen even in zero-dimensional models Rossi et al. (2016) and was first noticed in Bold DiagMC simulations of interacting fermions Kozik et al. (2015). Remembering the fact that for the large-$N$ $O\left(N\right)$ sigma model in the infinite space the expansion (13) converged with much higher precision than the discrepancy which we observe here for finite-size chiral chains (compare Figs. 6 and 4), we conclude that infinite lattice size and the presence of $\log\left(\lambda\right)$ terms seem to be crucial for the convergence to the correct physical result. An important extension of our analysis, which we relegate for future work, would be to extend the sufficient conditions of Rossi et al. (2016) for the convergence to the physical result to the weak-coupling expansion (13).

IV Schwinger-Dyson equations and DiagMC algorithm for $U\left(N\right)\times U\left(N\right)$ principal chiral model

The aim of this Section is to construct a DiagMC algorithm for sampling the weak-coupling expansion (13). While in this work we will only consider the weak-coupling expansion (13) based on the action (II) of the large-$N$ principal chiral model (2), our construction of the DiagMC algorithm is quite general and can be applied to both to the weak and strong-coupling expansions in general large-$N$ matrix field theories. We will briefly outline the extension to non-Abelian gauge theories in Section VI.

IV.1 Schwinger-Dyson equations for the large-$N$ $U\left(N\right)\times U\left(N\right)$ principal chiral model with Cayley parameterization of $U\left(N\right)$ group

As in Prokof’ev and Svistunov (2007); van Houcke et al. (2012); Tupitsyn and Prokof’ev (2017); Davody (2013); Pfeffer and Pollet (2017), the starting point for the construction of our DiagMC algorithm are the Schwinger-Dyson equations, which we derive for the action (II) of the principal chiral model (2). In contrast to the works Davody (2013); Pfeffer and Pollet (2017) on the bosonic DiagMC, here we consider the full untruncated hierarchy of Schwinger-Dyson equations for multi-trace, in general disconnected, correlators of the $\phi_{x}$ fields defined in (3):

[TABLE]

To derive the Schwinger-Dyson equations for the field correlators (IV.1), let’s consider the full derivative of the form

[TABLE]

where we have explicitly separated the action (II) into the quadratic and interacting parts given by (II) and (11), (II) to simplify the derivation and notation. We now expand the derivative according to the Leibniz product rule, contract matrix indices and finally contract the result with the bare propagator $G^{0}_{x^{1}_{1}\,x}=\left(\Delta+m_{0}^{2}\right)^{-1}_{xy}$ coming from the quadratic part $S_{F}\left[\phi_{x}\right]$ of the action (II). This yields the following infinite system of equations for the multi-trace observables (IV.1):

[TABLE]

[TABLE]

[TABLE]

where we have used red and violet colors to highlight the field arguments $x$ and the single-trace terms which are different on the left- and the right-hand sides. By $\frac{\partial S_{I}}{\partial\phi_{x}}$ we have schematically denoted the terms which are produced by applying the derivative operator in (43) to the interacting part (11) of the action (II):

[TABLE]

These terms contain three or more $\phi$ fields and should be inserted into the single-trace expectation values (IV.1) within the equations (44), (IV.1) and (IV.1): $\langle\,\left[\frac{\partial S_{I}\left[\phi\right]}{\partial\phi_{x}}\,y\right]\,\rangle=\langle\,\frac{1}{N}{\rm Tr}\,\left(\frac{\partial S_{I}\left[\phi\right]}{\partial\phi_{x}}\phi_{y}\right)\,\rangle$, and so on. For the sake of brevity, we omit the matrix indices of $\phi_{x}$, so that all products of $\phi_{x}$ are understood as matrix products. In particular, this implies that $\phi_{x}$ and $\phi_{y}$ do not commute for $x\neq y$.

In the t’Hooft large-$N$ limit, the last summands in the equations (IV.1) and (IV.1) can be omitted, since they are suppressed by a factor $1/N^{2}$ and all multi-trace correlators (IV.1) are of order $N^{0}$. Since in this limit the correlators (IV.1) with odd numbers of $\phi_{x}$ fields within some of the traces, such as $\langle\,\left[x\right]\left[y\right]\,\rangle\equiv\langle\,\frac{1}{N}{\rm Tr}\,\left(\phi_{x}\right)\frac{1}{N}{\rm Tr}\,\left(\phi_{y}\right)\,\rangle$, are suppressed as $1/N^{2}$, we can also limit the equations to the space of multi-trace correlator which contain an even number of $\phi$ fields within each trace. Furthermore, in the t’Hooft large-$N$ limit the equations (44), (IV.1) and (IV.1) without the $1/N^{2}$-suppressed terms admit the large-$N$ factorized solution

[TABLE]

with single-trace correlators $\langle\,\left[x_{1},\ldots,x_{n}\right]\,\rangle$ obeying the much simpler non-linear equation

[TABLE]

The form of the Schwinger-Dyson equation (44) for the two-point correlator does not change.

Schwinger-Dyson equations for large-$N$ quantum field theories are almost always considered in the nonlinear form (IV.1), which is much more compact than the linear equations (IV.1). One of the well-known examples of such non-linear Schwinger-Dyson equations in the large-$N$ limit are the Migdal-Makeenko loop equations for non-Abelian gauge theories Makeenko and Migdal (1981); Eguchi and Kawai (1982). Likewise, even at finite $N$ one can also arrive at a more compact non-linear representations of the Schwinger-Dyson equations (44), (IV.1) and (IV.1) by rewriting them in terms of connected or one-particle-irreducible correlators. This strategy is numerically very efficient for approximate solutions based on some truncations of the full set of equations (especially for correlators with small numbers of fields), and is commonly used for numerical solutions of QCD Schwinger-Dyson equations Alkofer and von Smekal (2001); Berges (2004) and in condensed matter physics Prokof’ev and Svistunov (1998); Van Houcke et al. (2010); Pfeffer and Pollet (2017). However, Schwinger-Dyson equations are always linear equations when written in terms of disconnected correlators of the form (IV.1). This linear form is very convenient for the construction of DiagMC algorithms, which we consider in the next Subsection.

To proceed towards the DiagMC algorithm, we go to the momentum space, and define the momentum-space field operators and correlators as

[TABLE]

Note that we have assumed the finite lattice volume $V$, which will be convenient for further analysis. In particular, this allows us to interpret all the delta-functions with momentum-valued arguments as Kronecker deltas on the discrete set of lattice momenta. For brevity, we will denote $\delta\left(p\right)\equiv\delta_{p,0}$. We will see in what follows that the algorithm performance is practically volume-independent, thus one can make a smooth transition from discrete to continuous momenta if necessary.

It is also convenient to introduce the momentum-space vertex functions

[TABLE]

so that the Fourier transform of the derivative (IV.1) of the interacting part of the action can be compactly represented as

[TABLE]

An important property of the vertex functions (IV.1) is that they are positive for all values of $q_{1},\ldots,q_{2v+1}$. While we don’t have an explicit proof of this statement, in our Monte-Carlo simulations with over $10^{14}$ Metropolis updates in total we have explicitly checked that the values of vertex function were always positive. Furthermore, in Appendix A we demonstrate that if all the momenta in (IV.1) are small and the lattice Laplacian can be approximated as $\Delta\left(p\right)=p^{2}$, the vertex functions take the simple form

[TABLE]

which manifestly demonstrates their positivity.

Furthermore, in order to define the power counting scheme to be used in our DiagMC simulations, we expand the momentum-space correlators of $\phi$ fields in formal power series in our auxiliary parameter $\xi$:

[TABLE]

In terms of the correlators (IV.1) of the momentum-space fields $\phi_{p}$, the large-$N$ limit of the Schwinger-Dyson equations (44), (IV.1) and (IV.1) read:

[TABLE]

where we have again used red and violet colors to highlight the individual momenta, order labels and momentum subsequences in correlators (IV.1) which are different on the right- and left-hand sides. We have also introduced the bare propagator in momentum space

[TABLE]

It is easy to convince oneself that recursive solution of the Schwinger-Dyson equations (55) generates the conventional weak-coupling expansion, with the auxiliary parameter $\xi$ formally serving as a coupling. For example, to find the order $m=2$ contribution to the two-point correlator $\langle\,\left[p_{1},p_{2}\right]\,\rangle$, we use the first equation in (55) to express it in terms of four-point correlator of order $m=1$ and the six-point correlator of order $m=0$. The four-point correlator at order $m=1$ can be related to two-point correlator at order $m=1$ and the six-point correlator at $m=0$ using the last equation in (55). Finally, the two-point correlator at order $m=1$ is again related to four-point correlator at $m=0$ by virtue of the first equation in (55). This chain of relations is schematically illustrated on Fig. 7.

It is instructive to explicitly calculate the first perturbative correction to the two-point correlator, which can be readily found from the Schwinger-Dyson equations (55) for the two-point correlator:

[TABLE]

where

[TABLE]

In particular, $\Sigma_{1}\left(p=0\right)=-2$. If we interpret this quantity as the first-order correction to the self-energy, then at zero momentum this correction, being multiplied by $\left(-\lambda/8\right)$ and $\xi=1$, exactly cancels the bare mass term $\lambda/4$! Of course, this cancellation requires re-summation of the infinite chain of bubble diagrams, and at any finite order of the expansion all diagram weights are still infrared-finite.

IV.2 Metropolis algorithm for stochastic solution of linear equations

As discussed in the previous Subsection, Schwinger-Dyson equations in terms of disconnected field correlators always form an infinite system of linear equations, which can be written in the following general form

[TABLE]

where $A\left(X|Y\right)$ is some theory-dependent linear operator, $b\left(X\right)$ represents the “source” terms in the Schwinger-Dyson equations (typically, the contact delta-function term in the lowest-order equation, like in (44) and (55), and $X$ and $Y$ are generalized indices incorporating all arguments of field correlators which enter the Schwinger-Dyson equations. For the Schwinger-Dyson equations (55) discussed in the previous Subsection, $X$ and $Y$ are sequences of momenta $\left\{\left\{p^{1}_{1},\ldots,p^{1}_{n_{1}}\right\},\ldots,\left\{p^{r}_{1},\ldots,p^{r}_{n_{r}}\right\}\right\}$ which parameterize the multi-trace correlators (IV.1) with an additional order-counting integer variable $m$. Note that the set of $X$ indices is in most cases infinite, which excludes the application of standard numerical linear algebra methods to the equation (59) (see however Anderson and Kruczenski (2016) for an attempt in this direction). Bringing the Schwinger-Dyson equations into the form (59) can be subtle for field theories for which the free (quadratic) part of the action has flat directions or vanishes, but these subtleties are not relevant for the present work.

The solution $\phi\left(X\right)$ of (59) can be written as formal geometric series in $A$:

[TABLE]

If the diagrammatic expansion which is obtained by iterating Schwinger-Dyson equations is truncated at some finite order, there is only a finite number of terms in these series, which correspond to a finite number of diagrams contributing at a given expansion order.

The key idea of the DiagMC algorithm which we devise in this work is to evaluate the series (IV.2) by stochastically sampling the sequences of generalized indices $\mathcal{S}=\left\{X_{n},\ldots,X_{0}\right\}$ with arbitrary $n$ in (IV.2) with probability

[TABLE]

where $\mathcal{N}_{w}$ is the appropriate normalization factor.

The solution $\phi\left(X\right)$ to the system (59) can be found by histogramming the last generalized index $X_{n}$ in the sequence $\left\{X_{n},\ldots,X_{0}\right\}$, where each occurrence of $X_{n}$ is weighted with the sign $\sigma\left(X_{n},\ldots,X_{0}\right)={\rm sign}\left(A\left(X_{n}|X_{n-1}\right)\right)\ldots{\rm sign}\left(A\left(X_{1}|X_{0}\right)\right){\rm sign}\left(b\left(X_{0}\right)\right)$. This sign reweighting puts certain limitations on the use of the method, since we are effectively sampling the series in which the linear operator $A\left(X|Y\right)$ is replaced by $|A\left(X|Y\right)|$ and which typically have a smaller radius of convergence. In the worst case, the expectation value of the reweighting sign $\sigma\left(X_{n},\ldots,X_{0}\right)$ can decrease exponentially or faster with $n$, thus limiting the maximal expansion order which can be sampled by the algorithm.

Let us also note that since for Schwinger-Dyson equations the generalized indices $X$ take infinitely many values, in practice histogramming of all possible values of $X$ is not possible, but also not necessary. E.g. if one is interested only in the momentum-space two-point correlator $\langle\,\frac{1}{N}{\rm Tr}\,\left(\phi_{p_{1}}\,\phi_{p_{2}}\right)\,\rangle=\langle\,\left[p_{1},p_{2}\right]\,\rangle$ which corresponds to $X$ of the form $X=\left\{\left\{p_{1},p_{2}\right\}\right\}$, one can histogram only the momenta $p_{1}$ and $p_{2}$ whenever $X$ contains only two of them. Thus no field correlators should be explicitly saved in computer memory, and the memory required by the algorithm does not directly scale with system volume (however, it does scale with the maximal expansion order, which might be correlated with the volume via required statistical error).

In order to simplify the notation in what follows, let us also define the quantities

[TABLE]

In this work we propose to sample the sequences $\mathcal{S}$ with probability (61) by using the Metropolis-Hastings algorithm with the following updates:

Algorithm IV.IV.2.1

Add index:

With probability $p_{+}$ a new generalized index $X_{n+1}$ is added to the sequence $\left\{X_{n},\ldots,X_{0}\right\}$, where the probability distribution of $X_{n+1}$ is $\pi\left(X_{n+1}|X_{n}\right)=|A\left(X_{n+1}|X_{n}\right)|/\mathcal{N}\left(X_{n}\right)$. The sign variable $\sigma$ is multiplied by ${\rm sign}\left(A\left(X_{n+1}|X_{n}\right)\right)$.

Remove index:

If the sequence contains more than one generalized index, with probability $1-p_{+}$ the last generalized index $X_{n}$ is removed from the sequence, which thus changes from $\left\{X_{n},X_{n-1},\ldots,X_{0}\right\}$ to $\left\{X_{n-1},\ldots,X_{0}\right\}$. The sign variable $\sigma$ is multiplied by ${\rm sign}\left(A\left(X_{n}|X_{n-1}\right)\right)$.

Restart:

If the sequence contains only one generalized index $X_{0}$, it is replaced by $X_{0}^{\prime}$ with probability $1-p_{+}$. The probability distribution of $X_{0}^{\prime}$ is $\pi\left(X_{0}^{\prime}\right)=|b\left(X_{0}^{\prime}\right)|/\mathcal{N}_{b}$. The sign variable $\sigma$ is set to ${\rm sign}\left(b\left(X_{0}^{\prime}\right)\right)$.

These updates should be accepted or rejected with the probability Hastings (1970)

[TABLE]

Explicit calculation yields the following acceptance probabilities for the three updates introduced above:

[TABLE]

After some algebraic manipulations the overall average acceptance rate can be expressed as

[TABLE]

where in the last line the expectation value is taken with respect to the stationary probability distribution of $X_{n}$. For the optimal performance of the algorithm it is desirable to choose $p_{+}$ (which is the only free parameter in our algorithm) in such a way that the acceptance rate is maximally close to unity. To this end we first perform Metropolis updates with some initial value of $p_{+}$, estimate the expectation values $\langle\,\min\left(p_{+}^{\prime},\left(1-p_{+}^{\prime}\right)\mathcal{N}\left(X_{n}\right)\right)\,\rangle$ in (IV.2) for several trial values $p_{+}^{\prime}$ and then select the value for which the acceptance is maximal. This process is repeated several times until the values of the average acceptance and $p_{+}$ stabilize.

Sometimes the transformations $X_{n}\rightarrow X_{n+1}$ can be joined into groups of very similar transformations. For example, if we have to generate a random loop momentum when sampling the Feynman diagrams in momentum space, in the notation of Algorithm IV.2 generating every single possible momentum would be a separate transformation, but it is more convenient to describe them as a single group of transformations. Denoting possible outcomes of such transformations as $X_{n+1}\in\mathcal{T}\left(X_{n}\right)$, where $\mathcal{T}\left(X_{n}\right)$ labels the group of transformations, we can express the probability distribution of the new element $X_{n+1}$ in the “Add index” update in (IV.2) as the product of the conditional probability $\pi\left(X_{n+1}|\mathcal{T}\left(X_{n}\right)\right)$ to generate the new element $X_{n+1}$ by one of the transformations in $\mathcal{T}\left(X_{n}\right)$, times the overall probability $\pi\left(\mathcal{T}|X_{n}\right)$ to perform any transformation from $\mathcal{T}$:

[TABLE]

Such grouping of transformations will be very convenient for a compact description of our DiagMC algorithm in what follows. Sometimes we will also omit the arguments of $\mathcal{N}_{\mathcal{T}}\left(X_{n}\right)$ and $\mathcal{T}\left(X_{n}\right)$ for the sake of brevity.

Finally, let us calculate the normalization factor $\mathcal{N}_{w}$ relating the (sign weighted) histogram of the last generalized index $X_{n}$ in the sequences $\mathcal{S}=\left\{X_{n},\ldots,X_{0}\right\}$ and the solution $\phi\left(X\right)$ of the linear equations (59). Combining (61) and (59), we obtain a linear equation for $\mathcal{N}_{w}$:

[TABLE]

where $w\left(X\right)$ is the probability distribution of the last element $X_{n}$ in the sequences $\mathcal{S}$, with any sign $\sigma$. The solution of this equation can be written as

[TABLE]

where again the expectation value is taken with respect to the stationary probability distribution of the Metropolis algorithm IV.2.

The general implementation of the Metropolis algorithm IV.2 for which the user-defined matrix elements $A\left(X|Y\right)$ and the source vector $b\left(X\right)$ should be provided is available as a part of GitHub repository Buividovich (2017).

IV.3 DiagMC algorithm for the large-$N$ $U\left(N\right)\times U\left(N\right)$ principal chiral model

In this Subsection we adapt the general Metropolis algorithm described in the previous Subsection IV.2 for solving the linear Schwinger-Dyson equations (55) obtained from the action (II) of the large-$N$ principal chiral model (2). The generalized indices $X$ will take values in the space of partitioned sequences of momenta

[TABLE]

additionally labelled with an expansion order $m=0,1,2,\ldots$. These sequences will be generated with the probability proportional to the correlators $\langle\,\left[p^{1}_{1},\ldots,p^{1}_{n_{1}}\right]\ldots\left[p^{r}_{1},\ldots,p^{r}_{n_{r}}\right]\,\rangle_{m}$ defined in (IV.1), up to the sign re-weighting with the sign variable $\sigma=\pm 1$, introduced in algorithm IV.2.

The system of equations (55) already has the form (59). The only non-zero elements of the source vector $b\left(X\right)$ correspond the delta-function term (in the right-hand side of line (55a)) in the first, inhomogeneous Schwinger-Dyson equation in (55):

[TABLE]

From equations (55) we can also directly read off the following groups of nonzero elements of the $A\left(X|Y\right)$ matrix:

[TABLE]

Here we have used red and violet colors to highlight the individual momenta and order labels and momentum subsequences which are different for the first and the second indices of the $A\left(X|Y\right)$ matrix.

In the above representation we write down all nonzero elements of the matrix $A\left(X|Y\right)$ for the fixed first index $X$. However, for the “Add index” update in the Algorithm IV.2 we need to know all nonzero elements $A\left(X_{n+1}|X_{n}\right)$ for the fixed second index $X_{n}$. In order to obtain a suitable representation of the matrix $A\left(X|Y\right)$, we re-label the matrix elements (71), assuming that the second index $Y$ in $A\left(X|Y\right)$ has the most general form $Y=\left\{\left\{p^{1}_{1},\ldots,p^{1}_{n_{1}}\right\},\ldots,\left\{p^{r}_{1},\ldots,p^{r}_{n_{r}}\right\}\right\}_{m}$, and express the momenta in the first index $X$ in terms of these momenta $p^{a}_{i}$. This gives

[TABLE]

where the alphabetic labels of matrix elements are in one-to-one correspondence with alphabetic labels in (71).

The matrix elements (72) - (72) form sets of elements parameterized by the momentum variable $p$. We denote the normalization factor $\mathcal{N}_{\mathcal{T}}\left(X_{n}\right)$ in (IV.2) for these sets as

[TABLE]

where $\Sigma_{0}$ has an obvious interpretation of the one-loop self-energy. Correspondingly, the probability distribution of the new momenta $p$ is

[TABLE]

At the first sight it might seem that the Metropolis algorithm (IV.2) would require keeping in memory all the generalized indices $X_{0},\ldots,X_{n}$ which enter the series (IV.2). However, this is not necessary in practice, since the matrix $A\left(X|Y\right)$ for the Schwinger-Dyson equations (55) is very sparse, as one can see from the equations (72). Also, only a few momenta in the generalized index $X_{n}=\left\{\left\{p^{1}_{1},\ldots,p^{1}_{n_{1}}\right\},\ldots,\left\{p^{r}_{1},\ldots,p^{r}_{n_{r}}\right\}\right\}$ are different between the first and the second generalized indices of the matrix $A\left(X_{n+1}|X_{n}\right)$. Thus only a few momenta will be changed when generating the new generalized index $X_{n+1}$ from $X_{n}$. Therefore it is sufficient to keep in computer memory only the topmost sequence $X_{n}$, as well as the sequence of transformations $X_{n}\rightarrow X_{n+1}$ which led to the current sequence, along with a few parameters which characterize these transformations. In fact, only the information which is necessary to recover $X_{n}$ from $X_{n+1}$ in the “Remove index” update in Algorithm IV.2 should be stored.

We note also that the transformations $X_{n}\rightarrow X_{n+1}$ which correspond to matrix elements (72), (71), (71), (72) involve only the first momentum subsequence $\left\{p^{1}_{1},\ldots,p^{1}_{n_{1}}\right\}$ in (69), and the transformation corresponding to (72) involves also the second subsequence $\left\{p^{2}_{1},\ldots,p^{2}_{n_{2}}\right\}$. This naturally endows the partitioned sequences (69) with the stack structure, in the algorithmic sense. Correspondingly, we will refer to the subsequence $\left\{p^{1}_{1},\ldots,p^{1}_{n_{1}}\right\}$ as the “topmost” or “first” subsequence, the subsequence $\left\{p^{2}_{1},\ldots,p^{2}_{n_{2}}\right\}$ - as the “second” sequence and so on.

Interestingly, the fact that the random sampling process operates on the stack implies the factorization of probabilities for subsequences in the stack Etessami and Yannakakis (2005); P.McKean (1967); Buividovich (2011a), which corresponds to factorization of single-trace correlators in the large-$N$ limit. If one proceeds keeping the $1/N^{2}$ terms in the equations (IV.1) and (IV.1), our partitioned momentum sequences will no longer have stack structure. In particular, the fourth line of the Schwinger-Dyson equations (IV.1) would correspond to a transformation where the topmost momentum subsequence is split into two subsequences, one of which is inserted into arbitrary subsequence in the partitioned sequence (69) Buividovich (2016). Obviously, the number of ways to do such a splitting and insertion grows polynomially both in the number of subsequences and their lengths. This polynomial growth of the number of terms in the Schwinger-Dyson equations results in the factorial growth of the number of non-planar Feynman diagrams with order Buividovich (2016, 2011b). In contrast, the matrix elements (72) do not depend on the number and lengths of subsequences, so that the algorithm samples a much smaller space of planar diagrams, the number of which grows only exponentially with expansion order. This is the reason why Diagrammatic Monte Carlo should be in general much more efficient in the large-$N$ limit.

Since the generalized indices $X$ take an infinite number of values, it is impossible but also not necessary to keep the matrix elements $A\left(X|Y\right)$ in computer memory. Rather, they can be implemented as a function which, given the topmost index $X_{n}$ in the series (IV.2), returns only nonzero values of $A\left(X_{n+1}|X_{n}\right)$ and the corresponding set of Metropolis proposals $X_{n}\rightarrow X_{n+1}$. The C code of the Metropolis algorithm IV.2 available at Buividovich (2017) relies exactly on such a “functional” implementation of the matrix $A\left(X_{n+1}|X_{n}\right)$.

With all these considerations in mind, from (72) we can finally deduce the following set of transformations $X_{n}\rightarrow X_{n+1}$ in the “Add index” update of Algorithm IV.2:

Transformation set IV.IV.3.1

Push a new pair of momenta to the stack,

from matrix elements (72):

• •

A new subsequence containing a pair of momenta {\color[rgb]{.5,0,.5}\left\{{\color[rgb]{1,0,0}p,-p}\right\}} is pushed to the top of the momentum stack. The probability distribution for {\color[rgb]{1,0,0}p} is \pi_{0}\left({\color[rgb]{1,0,0}p}\right), defined in (74).

• •

The order $m$ and the sign $\sigma$ do not change.

• •

The corresponding term in the Schwinger-Dyson equations is (55c).

• •

The overall weight for these transformations, summed over all {\color[rgb]{1,0,0}p}, is $\mathcal{N}_{\mathcal{T}}=\sum_{X_{n+1}\in\mathcal{T}}|A\left(X_{n+1}|X_{n}\right)|=\Sigma_{0}$, using the same notation as in (IV.2).

• •

To recover $X_{n}$ back from $X_{n+1}$, the topmost subsequence should be removed from the stack configuration $X_{n+1}$. No additional information should be saved for this backward step.

Add new momenta to the topmost subsequence,

from matrix elements (72), (72):

• •

A new pair of momenta is inserted into the topmost subsequence {\color[rgb]{.5,0,.5}\left\{p^{1}_{1},\ldots,p^{1}_{n_{1}}\right\}}, either a) at the beginning as {\color[rgb]{.5,0,.5}\left\{{\color[rgb]{1,0,0}p,-p,}\,p^{1}_{1},\ldots,p^{1}_{n_{1}}\right\}}, or b) at the beginning and at the end as {\color[rgb]{.5,0,.5}\left\{{\color[rgb]{1,0,0}p,}\,p^{1}_{1},\ldots,p^{1}_{n_{1}},{\color[rgb]{1,0,0}-p}\right\}}. The probability distribution for {\color[rgb]{1,0,0}p} is again \pi_{0}\left({\color[rgb]{1,0,0}p}\right).

• •

The order $m$ and the sign $\sigma$ do not change.

• •

The corresponding terms in the Schwinger-Dyson equations (55) are (55e), (55f).

• •

The overall weights for each of the transformations a) and b) are $\mathcal{N}_{\mathcal{T}}=\Sigma_{0}$.

• •

To recover $X_{n}$ back from $X_{n+1}$, either the two first momenta, or the first and the last momenta should be removed from the topmost subsequence in the stack. Correspondingly, one should remember whether the transformation a) or b) was chosen.

Merge two topmost subsequences,

from matrix elements (72):

• •

If there are two or more momentum subsequences in $X$, they are merged and interleaved with a pair of momenta {\color[rgb]{1,0,0}p,-p}, where {\color[rgb]{1,0,0}p} is generated with probability distribution \pi_{0}\left({\color[rgb]{1,0,0}p}\right):

[TABLE]

• •

The order $m$ and the sign $\sigma$ do not change.

• •

The corresponding term in the Schwinger-Dyson equations (55) is (55g).

• •

The overall weight for this transformation, summed over all {\color[rgb]{1,0,0}p}, is $\mathcal{N}_{\mathcal{T}}=\Sigma_{0}$.

• •

To recover $X_{n}$ back from $X_{n+1}$, one should split the topmost subsequence in the stack, which now has the form of the right-hand side of (75), at the position $n_{1}+1$, and remove the first momenta from the resulting subsequences. Correspondingly, the length {\color[rgb]{1,0,0}n_{1}} of the topmost subsequence in $X_{n}$ should be stored before the transformation $X_{n}\rightarrow X_{n+1}$.

Create $\left(2v+2\right)$-leg vertex,

from matrix elements (72):

• •

If there are four or more momenta in the topmost momentum subsequence in $X_{n}$, replace the first $2v+1$, $v=1\ldots n_{1}/2-1$ momenta by a single momentum equal to their sum:

[TABLE]

• •

The order {\color[rgb]{1,0,0}m} is increased by $v$, and the sign {\color[rgb]{1,0,0}\sigma} is multiplied by $-1$.

• •

The corresponding terms in the Schwinger-Dyson equations (55) are (55b), (55) and (55).

• •

The weight of this transformation for each particular $v$ is

[TABLE]

• •

To recover $X_{n}$ back from $X_{n+1}$, one should replace the first momentum in the topmost subsequence in the stack by the momenta {\color[rgb]{1,0,0}p^{1}_{1},\ldots,p^{1}_{2v+1}}. To this end, the vertex order {\color[rgb]{1,0,0}v} and the $2v$ momenta {\color[rgb]{1,0,0}p^{1}_{1},\ldots,p^{1}_{2v}} should be stored. The momentum {\color[rgb]{1,0,0}p^{1}_{2v+1}} can be recovered from momentum conservation constraint.

Restart,

from the source vector (70):

• •

If the “Restart” update is selected in Algorithm IV.2, the new value of $X_{0}$ has a single subsequence with just two momenta: X_{0}^{\prime}=\left\{{\color[rgb]{.5,0,.5}\left\{{\color[rgb]{1,0,0}p,-p}\right\}}\right\}, where {\color[rgb]{1,0,0}p} is distributed with probability \pi_{0}\left({\color[rgb]{1,0,0}p}\right).

• •

The order {\color[rgb]{1,0,0}m} is set to {\color[rgb]{1,0,0}m}=0, and the sign variable to {\color[rgb]{1,0,0}\sigma}=+1. This state is also the initial state for the algorithm.

• •

The corresponding term in the Schwinger-Dyson equations (55) is the term on the right-hand side of line (55a).

• •

The norm of the source vector is $\mathcal{N}_{b}=\sum\limits_{X_{0}^{\prime}}b\left(X_{0}^{\prime}\right)=\Sigma_{0}$, which can be used to recover the overall normalization factor for field correlators from (68).

The sequence of updates $X_{n}\rightarrow X_{n+1}$ in the Metropolis algorithm IV.2 with transformations IV.3 can be visually interpreted as a step-by-step process of drawing planar, in general disconnected Feynman diagrams with arbitrary number of external legs. The transformations IV.3 are used in the “Add index” update and correspond to adding a randomly chosen new element to the diagram. The first three transformations in the set IV.3 correspond to drawing free propagator lines with different positions of endpoints. The last transformation in the set IV.3 draws an interaction vertex with $2v+2$ legs by joining $2v+1$ external lines into a single external line. It is also this transformation which turns external lines into internal ones. The “Remove index” update in the Metropolis Algorithm IV.2 erases the diagram element which was drawn last, and can be thought of as an “Undo” operation. The initial state of the algorithm is just a single free propagator line, which is the simplest Feynman diagram one can ever draw. Simultaneous generation of random momenta distributed with the probability $\pi_{0}\left(p\right)$ proportional to the free propagator $G_{0}\left(p\right)$ implements the Monte-Carlo integration over internal loop momenta in Feynman diagrams. In this way the diagram weights are represented as products of propagators with freely generated momenta times vertex functions and propagators which contain some linear combinations thereof. Let us also note that the momentum conservation is automatically valid for all the transformations in (IV.3), which always lead to zero total momentum for any momentum subsequence in the stack.

Computationally, the most expensive operation in our DiagMC algorithm is the calculation of vertex functions $V\left(p^{1}_{1},\ldots,p^{1}_{2v+1}\right)$ for all $v=1\ldots n_{1}/2-1$, which is necessary to obtain acceptance probabilities for all possible “Create vertex” transformations in the set IV.3. From (IV.1) one can infer that for each $v$ the calculation of the vertex function involves $\sim v^{2}$ floating-point operations. As a result, the CPU time needed to calculate all vertex functions scales as $n_{1}^{3}$. This CPU time is dominating the performance of the algorithm in the regime sufficiently close to the continuum limit, where subsequences with rather large lengths $n_{1}$ can appear. It is not inconceivable that this scaling can be made milder by carefully optimizing the calculation of vertex functions and re-using some data from previous Metropolis updates. We leave such developments for future work.

After calculation of vertex functions, the next most expensive and frequently used part of the algorithm is the generation of random momenta with probability $\pi_{0}\left(p\right)$. To this end one can implement any method for generating continuous variables with definite probability distribution, for example, the biased Metropolis algorithm Bazavov and Berg (2005). Note that this allows to work directly in the infinite volume limit, which is one of the attractive features of DiagMC approach in general Van Houcke et al. (2010).

In this work we have still chosen to work with finite but sufficiently large spatial lattice size $L_{0}\times L_{1}$, which is much larger than the correlation length for all values of $\lambda$ which we have used. In this case we generate the momenta $p$ by randomly choosing among a finite number of discrete values with given (in general, unequal) probabilities. By creating an indexed “lookup table” for the discrete values of momenta this procedure can be performed in CPU time which scales only logarithmically with volume.

The only reason for working in a finite volume is that this simplifies the histogramming of correlators of $g_{x}$ variables (original $U\left(N\right)$-valued variables in the partition function (2)), see Subsection V.1 for a more detailed discussion. The performance of the algorithm, however, practically does not depend on the lattice size as long as it is much larger than the correlation length. Rather, the dynamically generated mass gap in the principal chiral model (2) effectively sets the infrared scale for the algorithm. Needless to say, finite lattice size $L_{0}$ in the Euclidean time direction $\mu=0$ is also necessary to study the theory at finite temperature (see Subsection V.4).

The non-convergence of the re-summed perturbation theory of Section II for principal chiral model on finite-size one-dimensional lattice suggests that the series obtained on finite-size lattices might have worse convergence properties than in the infinite-volume limit. We expect that these problems will appear only at very high orders of perturbation theory. Indeed, replacing continuum integrals over Brillouin zone in the infrared-finite weak-coupling expansion of Section II by discrete sums of the form $\frac{1}{L^{2}}\sum\limits_{p}$ with $p_{\mu}=\frac{2\pi m_{\mu}}{L}$, $m_{\mu}=0\ldots L-1$ results in corrections to the continuum integrals which decrease with $L$ approximately as $e^{-mL}$, but increase with diagram order. For example, for a chain of $l$ one-loop diagrams, such corrections depend on order $l$ and lattice size $L$ as $\left(1+ae^{-mL}\right)^{l}$, with some constant $a$. Thus if the expansion is limited to some maximal order $M$ which also limits the maximal number of loop momenta, finite-volume corrections can always be made smaller by using sufficiently large lattice volume $V$. In Subsection V.3 we explicitly demonstrate the smallness of finite-volume corrections by comparing the results of simulations at $108\times 108$ and $256\times 256$ lattices.

It seems also that our DiagMC algorithm does not allow for a straightforward parallelization. The only potential for parallelization which we currently see is to relegate the generation of random momenta with probability distribution $\pi_{0}\left(p\right)$ to a separate MPI or OpenMP thread. However, as our algorithm can operate directly in the infinite-volume limit, parallelization might not be really necessary. Having a multi-CPU computer at hand, one can speed up simulations by gathering statistics independently on different CPUs.

In the next Section V we discuss in detail the numerical results which we have obtained using the algorithm IV.2 with transformation set IV.3, as well as the details of data collection and processing. The production code with the transformations IV.3 in the Metropolis algorithm IV.2 is available as a part of GitHub repository Buividovich (2017).

V Numerical results of the DiagMC simulations of the large-$N$ $U\left(N\right)\times U\left(N\right)$ principal chiral model

V.1 Correlators of $g_{x}$ variables in terms of $\phi_{p}$ fields

While the DiagMC algorithm described in Subsection IV.3 stochastically samples the correlators of $\phi_{p}$ fields in momentum space, physically one is typically interested in correlators of $U\left(N\right)$-valued $g_{x}$ variables which enter the partition function (2).

The first non-trivial correlator is the expectation value $\langle\,\frac{1}{N}{\rm Tr}\,g_{x}\,\rangle$ of the trace of $g_{x}$. It should vanish if the $U\left(N\right)\times U\left(N\right)$ symmetry of the principal chiral model (2) is not broken. By virtue of the Mermin-Wagner theorem which prohibits spontaneous symmetry breaking in $D=2$ dimensions, this should be the case at any finite value of $N$. Exact solution of the $O\left(N\right)$ sigma model at $N\rightarrow\infty$ suggests also that spontaneous symmetry breaking does not happen for nonlinear sigma models in $D=2$ in the large-$N$ limit.

However, the infrared-finite weak-coupling expansion sampled by our DiagMC algorithm is built around the vacuum state with $g_{x}=1$, $\langle\,\frac{1}{N}{\rm Tr}\,g_{x}\,\rangle=1$, which explicitly breaks the $U\left(N\right)\times U\left(N\right)$ symmetry of the principal chiral model. To ensure that our DiagMC algorithm produces physically consistent results, it is thus important to check that the expectation value $\langle\,\frac{1}{N}{\rm Tr}\,g_{x}\,\rangle$ approaches zero as the maximal order of the weak-coupling expansion $M$ is extrapolated to infinity. We demonstrate this on Fig. 10 below.

Using the definition (3), it is straightforward to expand the expectation value $\langle\,\frac{1}{N}{\rm Tr}\,g_{x}\,\rangle$ in powers of the fields $\phi_{x}$. Since our perturbative series (13) was obtained from the same expansion, it is natural to extend our formal power counting scheme based on the auxiliary parameter $\xi$ also to correlators of $g_{x}$ variables, which leads to

[TABLE]

where in the first line we have used translational invariance and $\xi$ should be set to $\xi=1$ at the end of the calculation.

Yet another observable which we consider in this work is the two-point correlator $\langle\,\frac{1}{N}{\rm Tr}\,\left(g^{{\dagger}}_{x}\,g_{y}\right)\,\rangle$ of $g_{x}$ variables. This is the simplest nontrivial correlator from which one can infer the mean energy, the mass gap and the static screening length in the large-$N$ principal chiral model (2). Similarly to (V.1), we can use the definition (3) to express this two-point correlator in terms of $\phi$ fields as

[TABLE]

In order to perform a controllable extrapolation of the expansion order to infinity, we now limit our expansion in powers of $\xi$ to some maximal order $M$, which will serve as extrapolation parameter. Furthermore, we use the momentum-space field correlators (IV.1) which are sampled in our DiagMC simulations and explicitly take into account the momentum conservation, which implies translational invariance for the correlators of $g_{x}$ variables. This allows us to express the expectation values (V.1) and (V.1) as

[TABLE]

[TABLE]

In these expressions, we use the subscript $M$ for the correlators of $g_{x}$ variables to denote that our weak-coupling expansion (as defined in Section II) was truncated at order $M$.

While the algorithm IV.2 with the transformation set IV.3 can in principle sample correlators $\langle\,\left[p_{1},\ldots,p_{2k}\right]\,\rangle_{m}$ to arbitrarily high orders $m$, we have found that in practice one can achieve better precision by sampling only field correlators $\langle\,\left[p_{1},\ldots,p_{2k}\right]\,\rangle_{m}$ with $k+m\leq M$. Indeed, only these field correlators contribute to the expectation values (V.1) and (V.1). In order to limit Monte-Carlo sampling to $k+m\leq M$, we set to zero all matrix elements $A\left(X_{n+1}|X_{n}\right)$ for transformations which lead to $k+m>M$ in the Metropolis algorithm IV.2. In our simulations we have been able to measure the first $M=12$ coefficients of the infrared-finite weak-coupling expansion of Section II with reasonably small statistical error. Note also that once the field correlators $\langle\,\left[p_{1},\ldots,p_{2k}\right]\,\rangle_{m}$ are measured for all values of $k+m<M_{max}$ up to some maximal value $M_{max}$ of $M$, we can calculate at once the truncated correlators $\langle\,\frac{1}{N}{\rm Tr}\,g_{x}\,\rangle_{M}$ and $\langle\,\frac{1}{N}{\rm Tr}\,\left(g^{{\dagger}}_{x}g_{0}\right)\,\rangle_{M}$ for all $M\leq M_{max}$. This is very convenient in practice for the extrapolation of numerical results to the infinite expansion order $M\rightarrow\infty$.

In order to measure the field correlators (V.1) and (V.1), it is of course not necessary to histogram the sequences (69) with multiple momentum variables, which would require extremely large amount of RAM memory. The expectation value (V.1) involves only the sum $\sum\limits_{p_{1}\ldots p_{2k}}\langle\,\left[p_{1},\ldots,p_{2k}\right]\,\rangle_{m}$, which is proportional to the total probability to get any momentum sequence $\left\{\left\{p_{1},\ldots,p_{2k}\right\}\right\}_{m}$ of length $2k$ with order $m$. Thus we only need to histogram the non-negative integer variables $k$ and $m$, which can only take values with $k+m\leq M$. This requires a histogram with $M\left(M+1\right)/2$ bins, which takes a negligibly small amount of RAM memory for realistic values of $M$ ($M\leq M_{max}=12$ in our case).

Likewise, in order to measure the two-point correlator (V.1) at some fixed distance $x$, we have to histogram only the variables $k$ and $m$, additionally weighting each momentum sequence $\left\{\left\{p_{1},\ldots,p_{2k}\right\}\right\}$ with the factor $\Gamma\left(x;\,p_{1},\ldots,p_{2k}\right)$. Thus in order to measure the two-point correlator (V.1) we again need the histogram with $M\left(M+1\right)/2$ bins for every value of $x$ we are interested in. In order to avoid frequent calculation of the function $\Gamma\left(x;\,p_{1},\ldots,p_{2k}\right)$ which in practice takes significant amount of CPU time, we have calculated the two-point correlator $\sum_{p}e^{ipx}\langle\,\frac{1}{N}{\rm Tr}\,\left(g^{{\dagger}}_{x}g_{0}\right)\,\rangle_{M}$ in momentum space. This Fourier transformation turns the function of $\Gamma\left(x;\,p_{1},\ldots,p_{2k}\right)$ into the collection of $\delta$-functions

[TABLE]

which can be more easily incorporated into histogramming process (in deriving (V.1) we have taken into account the invariance of correlators under the replacement $x_{\mu}\rightarrow L_{\mu}-x_{\mu}$). After the momentum-space correlator $\sum_{p}e^{ipx}\langle\,\frac{1}{N}{\rm Tr}\,\left(g^{{\dagger}}_{x}g_{0}\right)\,\rangle_{M}$ is measured, it can be easily transformed back to coordinate space by an inverse Fourier transform. The implementation of this faster and simpler histogramming procedure was the main reason for working with finite lattice size in our simulations.

In order to translate the sign-weighted probabilities obtained by histogramming $k$ and $m$ into the physical values of correlators (IV.1), we have to multiply them with the normalization factor (68). While equation (68) suggests that this normalization factor can be calculated by measuring the expectation value $\langle\,\mathcal{N}\left(X_{n}\right)\,\rangle$ over the Monte-Carlo history, in practice we have found that $\mathcal{N}\left(X_{n}\right)$ is a very noisy observable which is likely to have a heavy-tailed probability distribution. Thus directly using the equation (68) to obtain the normalization factor $\mathcal{N}_{w}$ results in large statistical errors. We have used a more convenient practical strategy and fixed $\mathcal{N}_{w}$ to be the ratio of the weighted probability in the histogram bin with $k=1$, $m=1$ and the first perturbative correction to the two-point function $\sum_{p_{0},p_{1}}\langle\,p_{0}p_{1}\,\rangle_{1}$, which we have explicitly calculated using (57). As a cross-check, we made sure that using the value of $\mathcal{N}_{w}$ calculated in this way leads to the correct estimate of the next perturbative correction $\sum_{p_{0},p_{1}}\langle\,p_{0}p_{1}\,\rangle_{2}$, which we have also explicitly calculated.

V.2 Simulation setup, sign problem, and restoration of $U\left(N\right)\times U\left(N\right)$ symmetry

In our simulations we perform two scans over parameters of the principal chiral model (2). In the first scan, discussed in Subsection V.3, we fix the lattice size to $L_{0}\times L_{1}=108\times 108$ and study the dependence of physical observables on the t’Hooft coupling $\lambda$ in the range $\lambda=3.012\ldots 5.0$, which includes the scaling region of the principal chiral model (2) where the continuum physics sets in, as well as the point of the large-$N$ phase transition Gross and Witten (1980) between the weak- and the strong-coupling phases at $\lambda_{c}=3.27$ Campostrini et al. (1995). To check the finite-volume effects, we have also performed a single simulation with $L_{0}\times L_{1}=256\times 256$ lattice at $\lambda=3.1$. In the second scan, discussed in Subsection V.4, we fix the t’Hooft coupling to $\lambda=3.012$ and study the temperature dependence of physical observables by varying the temporal lattice size $L_{0}$ in the range $L_{0}=4\ldots 90$. As discussed above, we work on finite-volume lattices only to simplify the histogramming procedure, and the speed of our DiagMC algorithm depends on the lattice size at most logarithmically.

We perform histogramming of observables every time when our DiagMC algorithm reaches a configuration with only one momentum subsequence in the stack. One can also explicitly take into account the large-$N$ factorization (IV.1) and sample from multiple subsequences in the stack, however, this does not give a large advantage since in this case the samples appear to be rather strongly correlated. In order to avoid potential problems with autocorrelations, for each parameter set we run around $6\cdot 10^{4}$ simulations of $2\cdot 10^{8}$ Metropolis updates each with statistically independent seeds for the ranlux random number generator Lüscher (1994). The outcome of each of these simulations is treated as statistically independent data point, and the statistical error for physical observables is calculated as the standard error for this data set. The total number of Metropolis updates for each data set is thus around $10^{13}$, which takes around $10^{4}$ CPU-hours (single core). The speed of our algorithm can be probably increased by an order of magnitude by carefully optimizing the calculation of vertex functions and generation of random momenta. Let us also note that achieving a comparable precision for the large-$N$ extrapolation of the results of conventional Monte-Carlo simulations takes a comparable amount of CPU time Buividovich and Valgushev (2017).

In each of statistically independent simulations, during the first $5\cdot 10^{4}$ Metropolis updates we tune the probability $p_{+}$ of the “Add index” update in the Metropolis algorithm IV.2 in order to maximize the average acceptance rate. The optimal value of $p_{+}$ and the resulting acceptance rate take values in the range $p_{+}\simeq 0.4\ldots 0.6$, $\langle\,\alpha\,\rangle\simeq 0.5\ldots 0.65$ and depend rather weakly on the t’Hooft coupling and temperature. The average length of the sequence $\left\{X_{n},\ldots,X_{0}\right\}$ of generalized indices in the Metropolis algorithm IV.2 is in the range between $\langle\,n\,\rangle=17.3$ (for $\lambda=3.012$) and $\langle\,n\,\rangle=18.1$ (for $\lambda=5.0$). The average number of subsequences in the momentum stack on which the transformations IV.3 operate is in the range between $\langle\,r\,\rangle=1.2$ (for $\lambda=5.0$) and $\langle\,r\,\rangle=1.4$ (for $\lambda=3.012$). The dependence of both $\langle\,n\,\rangle$ and $\langle\,r\,\rangle$ on temperature is much weaker, and is completely within the ranges given above.

As the most pessimistic estimate for the autocorrelation time of our DiagMC algorithm, we can use the average number of Monte-Carlo steps between the “Restart” updates in the Metropolis algorithm IV.2. Indeed, the “Restart” update completely resets the state of the algorithm, and the sequence of generalized indices generated after the “Restart” update is completely statistically independent from the sequence which was generated before that. On Fig. 8 we illustrate the dependence of this “mean return time” on the t’Hooft coupling $\lambda$. It appears that the “Restart” updates are very rare for all values of $\lambda$. On average, they are separated by $\sim 10^{4}\ldots 10^{5}$ Monte-Carlo steps. The mean return time strongly increases towards larger values of the t’Hooft coupling $\lambda$. Together with the $\lambda$ dependence of $\langle\,n\,\rangle$, this observation shows that for larger $\lambda$ higher orders in the infrared-finite weak-coupling expansion of Section II become more important, and the algorithm spends more time sampling correlators $\langle\,p_{1}\ldots p_{2k}\,\rangle_{m}$ with larger values of $k+m$. Correspondingly, returns to the initial state of the algorithm with $k=1$, $m=0$ are less frequent. In the scan over finite temperatures, the mean return time practically does not depend on temperature, and only at the highest temperatures which correspond to $L_{0}=4$ and $L_{0}=8$ it exhibits a quick increase from approximately $3.5\cdot 10^{4}$ to $5.3\cdot 10^{4}$.

On Fig. 9 we illustrate the probability distribution of the length $k$ of the topmost subsequence in the stack and the order counter $m$ in our DiagMC algorithm for $\lambda=3.1$. The largest fraction of time (almost $20\%$) is spent on sampling the two-point correlator $\langle\,p_{0}p_{1}\,\rangle_{m}$ with the maximal value $m=M_{max}-1$ ( $=11$ in our case). The least probable configuration is $m=0$, $k=12$, which is the disconnected correlator with the maximal number of external legs which contributes to expectation values (V.1) and (V.1). In total, the algorithm spends most time sampling the higher-order diagrams. This is a completely reasonable behaviour, as Monte-Carlo integration over loop momenta at high values of $m$ requires a lot of statistics, and also sign cancellations are most important for higher-order Feynman diagrams.

Let us now present our numerical results. We start by demonstrating that the expectation value $\langle\,\frac{1}{N}{\rm Tr}\,g_{x}\,\rangle_{M}$ indeed tends to zero as the maximal order $M$ in the truncated weak-coupling expansion (V.1) becomes larger. This is an indication that in the limit of large $M$ the global $U\left(N\right)\times U\left(N\right)$ symmetry of the principal chiral model (2) is restored.

On Fig. 10 we show the dependence of $\langle\,\frac{1}{N}{\rm Tr}\,g_{x}\,\rangle_{M}$ on $1/M$ for several values of t’Hooft coupling and several temperatures. In all cases $\langle\,\frac{1}{N}{\rm Tr}\,g_{x}\,\rangle_{M}$ exhibits a nonlinear but monotonically decreasing dependence on $1/M$. For the largest t’Hooft couplings $\lambda=5.0,\,4.0$ and for the highest temperature $1/T=L_{0}=4$ the tendency of $\langle\,\frac{1}{N}{\rm Tr}\,g_{x}\,\rangle_{M}$ to vanish is rather clear. In contrast, it seems that at small $\lambda$ the dependence on $M$ should become very nonlinear at small $1/M$ for $\langle\,\frac{1}{N}{\rm Tr}\,g_{x}\,\rangle_{M}$ to extrapolate to zero. Indeed, the slope of $\langle\,\frac{1}{N}{\rm Tr}\,g_{x}\,\rangle_{M}$ becomes larger towards smaller values of $1/M$. These observations support the expectation that $\langle\,\frac{1}{N}{\rm Tr}\,g_{x}\,\rangle_{M}$ vanishes at $M\rightarrow\infty$.

As already mentioned in Subsection IV.3, not all coefficients $A\left(X_{n+1}|X_{n}\right)$ are positive for the transformation set IV.3. We thus have to use the sign reweighting, as discussed in Subsection IV.2. There are two potential sources of the sign problem which might hinder our DiagMC simulations: sign cancellations within the expectation values of the correlators $\langle\,\left[p_{1},\ldots,p_{2k}\right]\,\rangle_{m}$ of $\phi$ variables, and sign cancellations within the observables (V.1) and (V.1) themselves. In order to quantify sign cancellations in the correlators $\langle\,p_{1},\ldots,p_{2k}\,\rangle_{m}$, on Fig. 11 we plot the expectation value of the sign variable $\sigma$ in the Metropolis algorithm IV.2 as a function of the t’Hooft coupling $\lambda$. For all values of $\lambda$ the mean sign $\langle\,\sigma\,\rangle$ appears to be rather small, around $10^{-4}\ldots 10^{-5}$, thus sign problem is indeed important in our simulations. A somewhat encouraging feature is that the sign problem becomes milder towards the weak-coupling limit in which the theory approaches the continuum limit. However, from Fig. 10, as well as from the exactly solvable example of the 2D $O\left(N\right)$ sigma model (see Subsection III.2) we can conclude that at the same time the infrared-finite weak-coupling expansion of Section II converges less rapidly at small $\lambda$, so that more terms in the series should be sampled. The mean sign $\langle\,\sigma\,\rangle$ also depends rather weakly on temperature, and exhibits a rapid drop from $\langle\,\sigma\,\rangle\sim 5\cdot 10^{-5}$ to $\langle\,\sigma\,\rangle\sim 3\cdot 10^{-5}$ only at highest temperatures ($L_{0}=4$ and $L_{0}=8$). On Fig. 14 in the next Subsection we will also illustrate the sign problem for a particular physical observable (mean link).

V.3 Two-point correlators at zero temperature

In this Subsection we consider the simplest nontrivial two-point correlation function $\langle\,\frac{1}{N}{\rm Tr}\,\left(g^{{\dagger}}_{x}g_{0}\right)\,\rangle$ of $U\left(N\right)$-valued fields. From this correlator one can extract such important physical properties of the principal chiral model (2) as the energy density and the mass gap in the principal chiral model (2). It has been also extensively studied using conventional Monte-Carlo simulations of both $SU\left(N\right)\times SU\left(N\right)$ and $U\left(N\right)\times U\left(N\right)$ principal chiral models at several large values of $N$ Rossi and Vicari (1994a); Campostrini et al. (1995); Rossi and Vicari (1994b); Buividovich and Valgushev (2017). Weak-coupling expansion of $\langle\,\frac{1}{N}{\rm Tr}\,\left(g^{{\dagger}}_{x}g_{0}\right)\,\rangle$ was given up to three first orders in $\lambda$ at arbitrary values of $N$ in Rossi and Vicari (1994b). In this Subsection we will compare the large-$N$ extrapolations of these results with the results of our DiagMC simulations. This allows us to compare the performance of our algorithm with conventional Monte-Carlo simulations, as well as the convergence of our infrared-finite weak-coupling expansion with that of the standard lattice perturbation theory.

A particularly simple and important observable is the mean link $\langle\,\frac{1}{N}{\rm Tr}\,\left(g^{{\dagger}}_{1}g_{0}\right)\,\rangle$, which is also related to the energy density of the principal chiral model. On Fig. 12 we illustrate the convergence of the truncated series (V.1) for the mean link by plotting $\langle\,\frac{1}{N}{\rm Tr}\,\left(g^{{\dagger}}_{1}g_{0}\right)\,\rangle_{M}$ as a function of inverse truncation order $1/M$. The dependence on $1/M$ is nonlinear with the curvature increasing towards larger $M$, which requires some non-linear extrapolation procedure to reach the limit $M\rightarrow\infty$ ($1/M\rightarrow 0$). We have found that the dependence of the truncated series for the mean link on $1/M$ is well described by the following two extrapolating functions:

[TABLE]

We perform extrapolations by fitting the numerical data with these functions using the $\chi^{2}$ weight which includes the data covariance matrix to account for correlations between numerical values at different $M$. The fit parameters are $A$, $B$ and $C$, with $A$ obviously corresponding to the result of extrapolation to $M\rightarrow\infty$. Upon fitting, both extrapolating functions (83) and (84) lead to roughly the same values of $\chi^{2}$ over the number of degrees of freedom in the fit. These extrapolations are shown on Fig. 12 as red and blue strips for the functions (83) and (84), respectively. The strip width corresponds to the confidence interval set by the statistical uncertainty in the optimal values of $A$, $B$ and $C$. The difference between the two extrapolations (83) and (84) can be used as an estimate of a systematic error of our extrapolation procedure in the absence of other prior information which helps to choose between the two fits (see below).

To compare our data with the results of conventional Monte-Carlo simulations, on Fig. 12 we also show the mean link values obtained in Monte-Carlo simulations of both $SU\left(N\right)\times SU\left(N\right)$ and $U\left(N\right)\times U\left(N\right)$ principal chiral models as functions of $1/N^{2}$ (rescaled as $10/N^{2}$ for better readability of the plots). These data is taken from Rossi and Vicari (1994a); Campostrini et al. (1995); Rossi and Vicari (1994b) for t’Hooft couplings $\lambda=3.100,\,3.230,\,4.000$ and $N=6,\,9,\,15,\,21$ and from Buividovich and Valgushev (2017) for $\lambda=3.012$ and $N=6,\,9,\,12,\,18$. The finite-$N$ data is extrapolated to the large-$N$ limit using the linear fits of the form $\langle\,\frac{1}{N}{\rm Tr}\,\left(g^{{\dagger}}_{1}g_{0}\right)\,\rangle=A+B/N^{2}$, which are shown on Fig. 12 as straight solid green lines. We have to note that the Monte-Carlo data for $\lambda=3.100,\,3.230,\,4.000$ shows noticeable deviations from the expected $A+B/N^{2}$ behavior which significantly exceed the very small statistical errors in these data. Thus our $N\rightarrow\infty$ extrapolation might be also not completely accurate. Also the Monte-Carlo data for the finite-$N$ $U\left(N\right)\times U\left(N\right)$ principal chiral model, taken from Campostrini et al. (1995), does not seem to follow the $1/N^{2}$ scaling of finite-$N$ corrections.

For the three values $\lambda=3.0120,\,3.100$ and $\lambda=3.230$ the extrapolation of the finite-$N$ data for the $SU\left(N\right)\times SU\left(N\right)$ principal chiral model clearly approaches the result of the $M\rightarrow\infty$ extrapolation with the extrapolating function (83). While the differences of the $N\rightarrow\infty$ and $M\rightarrow\infty$ extrapolations lie within the statistical uncertainty of the $M\rightarrow\infty$ extrapolation of the DiagMC data, extrapolating function (83) has a small but noticeable tendency for over-estimating the mean link. On the other hand, the second extrapolating function (84) underestimates the mean link rather significantly, and we thus conclude that the first extrapolating function (83) yields more realistic results in the $M\rightarrow\infty$ limit.

To illustrate the ability of our DiagMC algorithm to work at arbitrarily large volumes without significant slow-down, on Fig. 13 we also compare the convergence of the truncated expansions (V.1) and (V.1) for the mean link $\langle\,\frac{1}{N}{\rm Tr}\,\left(g_{1}^{{\dagger}}g_{0}\right)\,\rangle$ and the single-field expectation value $\langle\,\frac{1}{N}{\rm Tr}\,g_{x}\,\rangle$ for lattice sizes $L_{0}\times L_{1}=108\times 108$ and $L_{0}\times L_{1}=256\times 256$. Both simulations took approximately the same CPU time, and had the same strength of the sign problem.

In order to compare the convergence of our massive lattice perturbation theory constructed in Section II with the convergence of the standard lattice perturbation theory constructed using massless propagators, on Fig. 12 we also show the one-loop, two-loop and three-loop ($M=1,\,2,\,3$, respectively) results of Rossi and Vicari (1994b) for the mean link at $N\rightarrow\infty$. It appears that conventional perturbation theory yields the results which are, at least for the first three orders, closer to the physical result and also seem to converge faster. Thus the price which we have to pay for the absence of infrared divergences is the slower convergence of our infrared-finite weak-coupling expansion.

For the largest value $\lambda=4.00$ which we consider our massive perturbation theory with both extrapolations (83) and (84) yields the result which is significantly larger than the large-$N$ extrapolation of the Monte-Carlo data. Also the standard lattice perturbation theory Rossi and Vicari (1994b) seems to converge to a larger value. This disagreement between weak-coupling expansions and the physical result is an indication of the large-$N$ transition between the weak-coupling and the strong-coupling regimes Gross and Witten (1980), which in the principal chiral model (2) happens at $\lambda_{c}=3.27$ Campostrini et al. (1995). After this transition, the physical values should be rather calculated using the strong-coupling expansion.

As we have already discussed in Subsection V.2 above, our simulations are limited to around $M\sim 10^{1}$ first orders of the weak-coupling expansion due to more and more important sign cancellations between different Feynman diagrams. To illustrate this sign problem in more detail, we have considered the sign cancellations within contributions of order $\left(-\frac{\lambda}{8}\right)^{m+k}$ with fixed $m+k=n$ to the mean link $\langle\,\frac{1}{N}{\rm Tr}\,\left(g^{{\dagger}}_{1}g_{0}\right)\,\rangle_{M}$ given by (V.1). We have characterized these sign cancellations by the ratio $\langle\,\sigma\,\rangle_{n}=\langle\,\frac{I_{n}^{+}-I_{n}^{-}}{I_{n}^{+}+I_{n}^{-}}\,\rangle$, where $I_{n}^{+}$ and $I_{n}^{-}$ are the sums of the absolute values of diagrams which contribute at order $n=k+m$ to the mean link (V.1) with positive and negative re-weighting signs $\sigma$, respectively, and $\langle\,\ldots\,\rangle$ denotes averaging over Monte-Carlo history. On Fig. 14 we illustrate that $\langle\,\sigma\,\rangle_{n}$ decays seemingly exponentially with $n$ for several values of the t’Hooft coupling $\lambda$. Since our algorithm spends most time sampling diagrams with large $n\sim M$, this explains the overall strength of the sign problem which we have previously illustrated on Fig. 11.

Finally, on Fig. 15 we compare the $x$ dependence of the two-point correlator $\langle\,\frac{1}{N}{\rm Tr}\,\left(g^{{\dagger}}_{x}g_{0}\right)\,\rangle_{M}$ at $\lambda=3.1$ with Monte-Carlo data obtained in Buividovich and Valgushev (2017) at $N=6$ and $N=9$. For DiagMC results the truncation order $M$ is encoded in color, from pure blue for $M=1$ to pure red for $M=12$. While at distances of a few lattice spacings the DiagMC data clearly approaches the Monte-Carlo data, at large distances the convergence of our massive perturbation theory becomes slower. Since the expression (V.1) for $\langle\,\frac{1}{N}{\rm Tr}\,\left(g^{{\dagger}}_{x}g_{0}\right)\,\rangle_{M}$ contains the constant $x$-independent contribution, it is not surprising that it decays down to some finite “plateau” value, rather than exponentially decaying to zero as the full correlators for $SU\left(6\right)$ and $SU\left(9\right)$ do. The height of this “plateau” is, however, monotonically decreasing with $M$. It seems that one has to go to much higher values of $M$ than we use to really see how the exponential decay of the $\langle\,\frac{1}{N}{\rm Tr}\,\left(g^{{\dagger}}_{x}g_{0}\right)\,\rangle$ emerges.

It is also interesting to compare the correlators obtained from our infrared-finite weak-coupling expansion with the results of the standard lattice perturbation theory, for which the first two orders in coordinate space were given in Rossi and Vicari (1994b). These results are shown on Fig. 15 as solid black (for the first order) and grey (for the second order) lines. As in the case of mean plaquette, we observe that at short distances the standard perturbative expansion is closer to the first-principle Monte-Carlo data than those from our infrared-finite weak-coupling expansion. The results of standard perturbative expansion show, however, a completely different behavior at large distances. At distance $x\sim 10$ they become negative, and further grow logarithmically with distance. In this sense in the far infrared they behave very differently from the exponentially decaying physical result. Furthermore, the second-order result of Rossi and Vicari (1994b) starts growing logarithmically towards positive infinity at $x\gtrsim 3\cdot 10^{3}$. While such large distances are certainly not interesting for numerical analysis, this discussion illustrates the infrared divergences of the standard lattice perturbation theory which become only stronger at higher orders.

V.4 Finite-temperature phase transition

As a further application of our DiagMC algorithm, in this Subsection we consider the principal chiral model (2) at finite temperature. Not much is known about the finite-temperature physics of this model. As for other asymptotically free quantum field theories, on general grounds one can expect that at sufficiently high temperature the $SU\left(N\right)$-singlet degrees of freedom are effectively “deconfined”. However, for the principal chiral models with $N>2$ this expectation has never been explicitly tested. One of the subtle points which make such tests difficult is the absence of a local order parameter for this “deconfinement” transition in the principal chiral model. From general physical intuition in such situations one expects either a first-order phase transition or a crossover. Understanding the “deconfinement transition” in the principal chiral model can be important for understanding the physics of one-dimensional compactifications of this model, which have been actively used in recent years for constructing resurgent trans-series Cherman et al. (2014). An explicit analytic calculation Cubero (2015) of the correlation functions in the large-$N$ $SU\left(N\right)\times SU\left(N\right)$ principal chiral model based on the Bethe ansatz techniques, unfortunately, does not allow so far to study the thermodynamics of this model.

In a Monte-Carlo study Buividovich and Valgushev (2017) co-authored by one of us the numerical results for the energy density and the correlation length at different temperatures are presented for the principal chiral model (2) at quite large but finite $N$, for the same t’Hooft coupling $\lambda=3.012$ and the spatial lattice size $L_{1}=108$ as in the present study. The study Buividovich and Valgushev (2017) indicates a weak enhancement of correlation length at some critical temperature, which very slowly tends to become sharper towards larger $N$. No jump or hysteresis of the free energy indicative of a first-order phase transition was observed.

In order to get further insights into the nature of this finite-temperature phase transition or crossover in the principal chiral model, we have performed DiagMC simulations at several different temperatures, which correspond to different lattice sizes $L_{0}=4\ldots 108$ in the Euclidean time direction at fixed t’Hooft coupling $\lambda=3.012$ (and hence fixed lattice spacing, which only depends on $\lambda$). As a physical observable which is most sensitive to long-range correlations of $g_{x}$ fields we have chosen the value of the two-point correlator $\langle\,\frac{1}{N}{\rm Tr}\,\left(g^{{\dagger}}_{x}g_{0}\right)\,\rangle_{M}$ at $x=\left(0,L_{1}/2\right)$, which is the point most distant from the source at $\left(0,0\right)$ on a periodic lattice of spatial size $L_{1}$. This point is also the “midpoint” of the correlator where it takes the minimal value. To get a clearer signal, we subtract the zero-temperature data from the finite-temperature data. The result is shown on Fig. 16 for the truncation orders $M=2,\,6,\,10$. The spatial lattice size is again $L_{1}=108$.

As the truncation order $M$ is increased, the difference between the finite-temperature and zero-temperature correlator seems to grow. However, the statistical errors also grow with $M$ and for most values of $L_{0}$ the differences lie within statistical errors. However, for the smallest values $L_{0}\simeq 4,\,8$ the difference clearly decreases, indicating the decrease of correlation length at very high temperatures which is probably related to Debye screening. Furthermore, for $L_{0}$ roughly between $35$ and $45$ the difference seems to grow with $M$ beyond statistical uncertainty, resulting in a sort of peak structure around $L_{0}\simeq 40$. This enhancement of correlations was also reported in Buividovich and Valgushev (2017) at the same value of $L_{0}$. It can be interpreted is an indication of a weak finite-temperature phase transition. Since in DiagMC simulations the correlator $\langle\,\frac{1}{N}{\rm Tr}\,\left(g^{{\dagger}}_{x}g_{0}\right)\,\rangle_{M}$ quickly grows with $M$, from our data we cannot exclude the scenario of finite-order phase transition at which the correlation length would diverge.

VI Extension to gauge theories

In this Section we discuss possible extension of the approaches outlined in this work to non-Abelian gauge theories at finite density, which is the ultimate motivation of this work. The most straightforward extension would be to construct the infrared-finite weak-coupling expansion of Section II for pure $U\left(N\right)$ lattice gauge theory with the partition function

[TABLE]

Using the Cayley map parameterization (3) and the $U\left(N\right)$ integration measure (II) in terms of $N\times N$ Hermitian matrices $A_{x,\mu}$, one can calculate the few first orders of the expansion of the action in powers of the t’Hooft coupling $\lambda$:

[TABLE]

where $F_{x,\mu\nu}=\nabla_{\mu}A_{x,\nu}-\nabla_{\nu}A_{x,\mu}-2i\alpha\left[A_{x,\mu},A_{x,\nu}\right]$ is the lattice field strength tensor constructed from the non-compact field $A_{x,\mu}$ and $\nabla_{\mu}$ is the forward lattice derivative. Again we have to choose $\alpha^{2}\sim\lambda$ to arrive at the canonical normalization of the kinetic terms in the action. Higher-order terms in the action of the $A_{x,\mu}$ gauge fields are either suppressed by higher powers of $\alpha$ or contain higher lattice derivatives which are suppressed in the naive continuum limit. Again we see that the “gluon” field $A_{x,\mu}$ has acquired the mass term proportional to the t’Hooft coupling $\lambda$. The action (VI) is of course still invariant under gauge transformations $g_{x,\mu}\rightarrow\Omega_{x}g_{x,\mu}\Omega^{{\dagger}}_{x+\hat{\mu}}$, which, however, now have quite a non-trivial nonlinear form in terms of $A_{x,\mu}$ variables. In particular, they are different from the conventional continuum gauge transformations $A_{x,\mu}\rightarrow A_{x,\mu}+\nabla_{\mu}\phi_{x}-i\left[A_{x,\mu},\phi_{x}\right]$. This is why the nonzero gluon mass in (VI) does not contradict the gauge invariance of the theory. It is interesting here to draw parallels with the gluon mass which one can obtain from non-perturbative solution of (truncated) Schwinger-Dyson equations and Ward identities of continuum theory Cornwall (1982).

The inclusion of finite-density quarks would require sampling of both fermionic and bosonic correlators. Since fermionic propagators at finite density are in general complex-valued, this would require additional re-weighting and make the sign problem worse. For perturbative expansions with truncation error which scales exponentially $\epsilon_{M}\sim\alpha^{M}$, $|\alpha|<1$ with the truncation order $M$ this still allows to obtain the result in a computational time which scales as a polynomial of the required error, provided one can sample Feynman diagrams in a time which grows only exponentially with expansion order Rossi et al. (2017). In the large-$N$ limit, the number of Feynman diagrams grows exponentially with their order, thus even the direct summation will satisfy this requirement. However, for our infrared-finite weak-coupling expansion (13), the truncation error scales as $\epsilon_{M}\sim M^{-\gamma}$, $\gamma>0$ (see (83), (84)), thus computational time will grow faster then polynomial in the required error, similarly to what happens for conventional Monte-Carlo simulations with sign problem. An important direction for further work is thus to find the re-summation strategy which would lead to exponentially fast convergence of the series, with bold diagrammatic expansions being one of the options.

A less straightforward extension of this work to non-Abelian gauge theories might be based on the lattice strong-coupling expansion, combined with hopping expansion if dynamical quarks are included and the Veneziano large-$N$ limit is taken. Inside the convergence radius of the strong-coupling and hopping expansions, their truncation error should scale down exponentially with truncation order. At the same time, in the large-$N$ limit the number of strong-coupling expansion diagrams should grow not faster than exponentially with order O’Brien and Zuber (1985). Furthermore, for hopping expansion adding chemical potential does not introduce any additional sign cancellations between the diagrams, but merely modifies the hopping coefficients in the time direction Unger and de Forcrand (2011); de Forcrand et al. (2014); Unger (2014); de Forcrand et al. (2016). Thus even despite the sign cancellations which happen between strong-coupling expansion diagrams at higher orders, DiagMC algorithms based on strong-coupling expansion can be expected to solve, at least formally, the computational complexity problem (in the terminology of Rossi et al. (2017)) for finite-density gauge theories in the Veneziano large-$N$ limit.

The Metropolis algorithm IV.2 allows to systematically sample the strong-coupling expansion for large-$N$ non-Abelian gauge theories. To this end it should be applied to solve the gauge-invariant Schwinger-Dyson equations for lattice Wilson loops $\langle\,\frac{1}{N}{\rm Tr}\,\left(g_{x_{1},\mu_{1}}g_{x_{2},\mu_{2}}\ldots g_{x_{n},\mu_{n}}\right)\,\rangle$ (Migdal-Makeenko loop equations Makeenko and Migdal (1981); Eguchi and Kawai (1982)). These equations are very similar in structure to equations (IV.1) in Subsection IV.1.

In fact, we have already tried to solve the Migdal-Makeenko loop equations in large-$N$ pure gauge theory using the algorithm IV.2, sampling the loops on the lattice with probabilities proportional to Wilson loops $\langle\,\frac{1}{N}{\rm Tr}\,\left(g_{x_{1},\mu_{1}}g_{x_{2},\mu_{2}}\ldots g_{x_{n},\mu_{n}}\right)\,\rangle$. Unfortunately, this straightforward attempt did not allow to go beyond the first few orders of strong-coupling expansion, which could anyway be computed manually. The reason is the “zigzag symmetry” Polyakov (1998); Drukker et al. (1999) between Wilson loops which differ by the insertion of arbitrary path traversed forward and immediately backward (see Fig. 17), which stems from unitarity of link matrices $g_{x,\mu}g_{x+\mu,-\mu}\equiv g_{x,\mu}g^{{\dagger}}_{x,\mu}=1$. Most of the computer time goes into sampling numerous redundant loops with multiple “zigzag” insertions, and nontrivial loops which contribute to higher orders of strong-coupling expansion are entropically disfavoured. Possible solution to this problem is to rewrite the Migdal-Makeenko loop equations in terms of “irreducible” loops without zigzag insertions, which makes the equation structure significantly more nontrivial. We hope that this “gauge-fixing” in loop space would allow to go to sufficiently high orders of strong-coupling expansion. Work in this direction is in progress.

Another challenge is to extend this approach to finite $N$, where the combinatorial structures involved in the strong-coupling expansion coefficients become more complicated Unger (2014); Samuel (1980), and $1/N$ expansion might not work properly Samuel (1980). Importantly, the baryon bound states only exist at finite $N$. Nevertheless, even the finite-density gauge theory in the large-$N$ Veneziano limit might exhibit some universal features for densities below the baryon condensation threshold and is therefore physically interesting Alho et al. (2014); Gürsoy (2016). Needless to say, large-$N$ limit is also important in the context of AdS/CFT holographic duality between gauge theories and gravity.

VII Conclusions

In this work we have explored the advantages and disadvantages of the Diagrammatic Monte-Carlo (DiagMC) approach for studying asymptotically free non-Abelian lattice field theories in the weak-coupling regime. To this end we have combined two largely independent tools: the Metropolis algorithm for solving infinitely-dimensional Schwinger-Dyson equations, and the partial resummation of lattice perturbation theory which leads to massive bare propagators.

While our primary motivation are non-Abelian lattice gauge theories, eventually coupled to finite-density fermions, in this exploratory study we have considered an example of the large-$N$ $U\left(N\right)\times U\left(N\right)$ principal chiral model. This model has significantly simpler weak-coupling expansion, while having most of the non-perturbative features of the large-$N$ pure gauge theories on the lattice.

The construction of Metropolis algorithm IV.2 from Schwinger-Dyson equations is very general, and can be applied to practically any quantum field theory. Its main conceptual advantages are:

• •

The applicability of algorithm to any expansion, even without explicit diagrammatic rules. In particular, it can be used to sample both weak- and strong-coupling expansions in lattice gauge theory.

• •

The ability to work directly in the infinite-volume limit. Of course, this does not imply the absence of slowing down towards the continuum limit, since the number of expansion coefficients which should be sampled to reach a given precision might grow with correlation length (see e.g. Subsection III.2), and sign problem might make this growth faster (see Subsection V.2).

• •

The ability to work directly in the large-$N$ limit, where the number of diagrams contributing to either strong- or weak-coupling expansions becomes significantly smaller. Systematic calculation of the coefficients of $1/N$-expansion is also possible Buividovich (2016). Large-$N$ limit is not directly accessible in most other DiagMC simulations of non-Abelian gauge theories Unger and de Forcrand (2011); de Forcrand et al. (2014); Vairinhos and de Forcrand (2014); Unger (2014); Gattringer and Marchis (2016); de Forcrand et al. (2016).

Our algorithm is particularly suitable for exploring the critical behavior in large-$N$ quantum field theories with manifestly positive weights of Feynman diagrams, for example, the planar $\phi^{4}$ model with negative coupling constant or the Weingarten model of random planar surfaces on the lattice Weingarten (1980).

The resummation strategy used in this work also seems to be quite general for non-Abelian field theories. Its particular advantages are:

• •

The absence of infrared divergences which are characteristic for the standard weak-coupling expansion of massless quantum field theories. While these divergences cancel in physical observables, they should be regulated at intermediate steps, which enhances the sign problem upon removing the infrared regulator due to cancellations between numerically large summands of opposite signs. In more conventional approaches to lattice perturbation theory such as Numerical Stochastic Perturbation Theory (NSPT) Bali et al. (2014); Bali and Pineda (2016) dealing with infrared divergences also requires additional infrared regulators.

• •

The absence of factorial divergences in the weak-coupling expansion in the large-$N$ limit, which arise at sufficiently high orders in conventional lattice perturbation theory. The general arguments of Section II, however, does not exclude the possibility of exponential divergences. Our arguments in favour of convergence of our infrared-finite weak-coupling expansion are only based on the numerical evidence - for up to $M=500$ orders for the large-$N$ $O\left(N\right)$ sigma model, and for up to $M=12$ orders for the principal chiral model.

Of course, these advantageous features come at a certain price, and in some respects the resummed perturbative expansion behaves worse than the bare one:

• •

Sampling our weak-coupling expansion using a DiagMC algorithm leads to the sign problem due to in general non-positive weights of Feynman diagrams, even though the conventional Monte-Carlo simulations do not have any sign problem. Furthermore, for lattice gauge theory this sign problem is expected to become worse if finite-density fermions are added.

• •

The convergence is slower than for the conventional lattice perturbation theory, see Figs. 12 and 15.

• •

The expansion starts from the wrong vacuum state with explicitly broken $U\left(N\right)\times U\left(N\right)$ symmetry of the principal chiral model. While this symmetry seems to be gradually restored at high orders, at finite expansion order this explicit symmetry breaking results in constant contribution to the correlators of $g_{x}$ variables, which in reality are exponentially decaying.

Summarizing these advantages and disadvantages, it seems that the DiagMC approach based on the weak-coupling expansion is comparable in terms of computational time to conventional simulation methods as applied to non-Abelian lattice field theories in the large-$N$ limit. In particular, for Fig. 12 the finite-$N$ data from conventional Monte-Carlo and the DiagMC data at infinite $N$ were produced within the comparable amount of CPU time. It is also conceivable, for example, that applying Numerical Stochastic Perturbation Theory (NSPT) Bali et al. (2014); Bali and Pineda (2016) to finite-$N$ principal chiral model and then extrapolating to large-volume and large-$N$ limits would lead to the comparable quality of numerical data within the comparable CPU time. It remains to be seen how efficient the sampling of strong-coupling expansion might be.

Probably the most promising direction for further improvement is the Bold Diagrammatic Monte-Carlo (BDMC) approach, which can efficiently capture non-perturbative information by combining perturbative expansion with a self-consistent mean-field-like approximation Pollet et al. (2010); van Houcke et al. (2012); Davody (2013); Pfeffer and Pollet (2017). Unfortunately, currently most of BDMC algorithms for bosonic fields are based on certain truncations of Schwinger-Dyson equations even for the simplest case of scalar field theory with $\phi^{4}$ interactions Davody (2013); Pfeffer and Pollet (2017).

One possibility to construct the bold diagrammatic expansion would be to work with Schwinger-Dyson equations for $n$-particle-irreducible field correlators Berges (2004). Another possibility, more specific for non-Abelian field theories, is to introduce the matrix-valued Lagrange multiplier enforcing the unitarity constraint $g^{{\dagger}}_{x}g_{x}=1$, and integrate out the $g_{x}$ fields, which would automatically ensure the $U\left(N\right)\times U\left(N\right)$-invariant vacuum state with $\langle\,\frac{1}{N}{\rm Tr}\,g_{x}\,\rangle=0$. For the large-$N$ $O\left(N\right)$ sigma model, a similar treatment immediately yields the exact non-perturbative solution, with the mass being equal to the saddle-point value of the Lagrange multiplier field. In fact, this mean-field solution resums our double series (III.1) with coupling and logs of coupling into a single exponential $\exp{\left(-\frac{4\pi}{\lambda}\right)}$, which is the only term in the true trans-series for the mass gap in the large-$N$ $O\left(N\right)$ sigma model. One could thus expect that if such approach works also for the principal chiral model, it might re-sum many terms in our expansion into exponentials, thus turning it into the true trans-series of the form (1) with optimal convergence properties and reducing the sign problem due to cancellations between Feynman diagrams. Implementing this approach in practice turns out to be non-trivial, and would be explored in future work.

Acknowledgements.

This work was supported by the S. Kowalevskaja award from the Alexander von Humboldt Foundation. The authors are grateful to the Mainz Institute for Theoretical Physics (MITP) for its hospitality and its partial support during the completion of this work. The calculations were performed on “iDataCool” at Regensburg University, on the ITEP cluster in Moscow and on the LRZ cluster in Garching. We acknowledge valuable discussions with S. Valgushev, F. Werner, B. Svistunov, T. Sulejmanpasic, L. von Smekal, N. Prokof’ev, G. Dunne, O. Costin, A. Cherman and G. Bali.

Appendix A Vertex functions in the small momentum limit

The vertex functions (IV.1) can be drastically simplified if all the momenta in their arguments are much smaller than the lattice cutoff scale. In this limit the lattice Laplacian $\Delta\left(p\right)$ behaves as $\Delta\left(p\right)\approx p^{2}$, and the vertex function $V(q_{1},\cdots,q_{2n+1})$ can be represented as

[TABLE]

The first term in the last line of the above equation can be simplified as

[TABLE]

where $\theta$ is the unit step function. Similar algebra gives for the second term in the last line of (A):

[TABLE]

Combining (A) and (A), we get

[TABLE]

which is the formula (53) in the main text.

Appendix B Integration measures on $S_{N}$ and $U\left(N\right)$ manifolds

B.1 Cayley map for $U\left(N\right)$ group

With the Cayley mapping $g=\frac{1+i\phi}{1-i\phi}=\frac{2}{1-i\phi}-1$, the group-invariant distance element $ds^{2}={\rm Tr}\,\left(dg\,dg^{{\dagger}}\right)$ on $U\left(N\right)$ can be expressed as

[TABLE]

The invariance of $ds^{2}$ under shifts $g\rightarrow ug$, $u\in U\left(N\right)$ implies, in particular, the invariance of the measure under similarity transformations $g\rightarrow ugu^{{\dagger}}$ which allow to diagonalize $g$ and hence the matrix $\phi$ in (B.1). We thus assume that the matrix $\phi_{ij}$ has diagonal form $\phi_{ij}=\phi_{i}\delta_{ij}$. The distance element $ds^{2}$ can be then written as a convolution of the coordinate differentials $d\phi_{ij}$ with the diagonal metrics:

[TABLE]

where we have taken into account that the off-diagonal components $\phi_{ij}$ are not all independent, but rather satisfy $\phi_{ij}=\bar{\phi}_{ji}$. We can now express the integration measure, up to an irrelevant normalization factor, as $\int dg=\int d\phi\sqrt{g}$, where $g$ is the determinant of the metric tensor which is equal to the product of all the diagonal metric elements appearing in (B.1):

[TABLE]

where we have denoted $f_{i}=\left(1+\phi_{i}^{2}\right)^{-1}$. Combining everything together, we arrive at

[TABLE]

where on the right hand side the integration goes over all $N\times N$ Hermitian matrices. This expression is identical to (II) in the main text, up to the trivial rescaling $\phi\rightarrow\alpha\phi$.

B.2 Exponential map for $U\left(N\right)$ group

Proceeding in the same way as for the Cayley map, for the exponential map $g=e^{i\phi}$ we can relate the differentials of $\phi$ and $g$ as

[TABLE]

The group-invariant distance element takes the form

[TABLE]

Again, we can now use the invariance of the measure under similarity transformations $g\rightarrow ugu^{{\dagger}}$, $\phi\rightarrow ugu^{{\dagger}}$, $u\in U\left(N\right)$ to diagonalize $g$ and $\phi$, and express the distance element (B.2) as

[TABLE]

We see that the metric tensor is again diagonal with our choice of group coordinates, and we finally get for the integration measure (up to an overall normalization factor)

[TABLE]

where $\mathbb{M}$ is some closed domain in the $N^{2}$-dimensional space $\mathbb{H}_{N\times N}$ of the Hermitian $N\times N$ matrices. The $U\left(N\right)$ group manifold is uniquely covered by the integration in (98) if $\mathbb{M}$ is bounded by the $N^{2}-1$-dimensional manifolds at which at least two eigenvalues of $\phi$ coincide. E.g. in the case of $U\left(2\right)$ group, $\mathbb{M}$ is the direct product of the $3$-dimensional ball and the one-dimensional circle $S_{1}$.

In order to compare with the formula (II) in the main text, it is instructive to expand the Jacobian (98) in powers of $\phi$:

[TABLE]

which illustrates the appearance of double-trace terms already in the lowest order of the expansion.

B.3 Stereographic map for the sphere $S_{N}$

In order to find the integration measure on the $N$-sphere $S_{N}$ parameterized by the coordinates $\phi_{i}$ as $n_{0}=\left(1-\phi^{2}\right)\left(1+\phi^{2}\right)^{-1}$, $n_{i}=2\phi_{i}\,\left(1+\phi^{2}\right)^{-1}$, $\phi^{2}\equiv\sum_{i}\phi_{i}^{2}$, we again start by first expressing the $O\left(N\right)$-invariant distance element $ds^{2}=dn_{a}\,dn_{a}$ on $S_{N}$ in terms of the coordinates $\phi_{i}$:

[TABLE]

The metric tensor is again diagonal in the stereographic coordinates $\phi_{i}$ (which is not surprising, as stereographic mapping is conformal), and the integration measure (up to a normalization constant) can be written in terms of the square root of the metric determinant as

[TABLE]

where the integration is over the whole $N-1$ dimensional real space $\mathbb{R}^{N-1}$. In the large-$N$ limit one can also replace the power of $-\left(N-1\right)$ in (101) by $-N$, which leads to the expression (29) in the main text after the trivial re-scaling $\phi\rightarrow\sqrt{\lambda}/2\,\phi$.