Complex Paths Around The Sign Problem
Andrei Alexandru, Gokce Basar, Paulo F. Bedaque, Neill C. Warrington
CComplex paths around the sign problem
Andrei Alexandru ∗ Department of PhysicsThe George Washington UniversityWashington, DC 20052
G¨ok¸ce Ba¸sar † Department of PhysicsUniversity of North CarolinaChapel Hill, NC 27599
Paulo F. Bedaque ‡ Department of PhysicsUniversity of MarylandCollege Park, MD 20742
Neill C. Warrington § Institute for Nuclear TheoryUniversity of WashingtonSeattle, WA 98195
The Monte Carlo evaluation of path integrals is one of a few general purpose meth-ods to approach strongly coupled systems. It is used in all branches of Physics, fromQCD/nuclear physics to the correlated electron systems. However, many systems ofgreat importance (dense matter inside neutron stars, the repulsive Hubbard model awayfrom half-filling, dynamical and non-equilibrium observables) are not amenable to theMonte Carlo method as it currently stands due to the so-called “sign-problem”. Wereview a new set of ideas recently developed to tackle the sign problem based on thecomplexification of field space and the Picard-Lefshetz theory accompanying it. Themathematical ideas underpinning this approach, as well as the algorithms so far devel-oped, are described together with non-trivial examples where the method has alreadybeen proved successful. Directions of future work, including the burgeoning use of ma-chine learning techniques, are delineated.
CONTENTS
I. The sign problem 2A. Field Theory/Many-Body Physics as a path integral 2B. Physical systems with sign problems 3C. Reweighting and the sign problem 3D. The absence of a general solution 4E. A brief survey of methods to deal with sign problems 4II. Cauchy theorem, homology classes and holomorphic flow 5A. Deformation of domain of integration: a multidimensionalCauchy theorem 5B. Holomorphic gradient flow 7C. Lefschetz Thimbles and Picard-Lefschetz theory 7III. Algorithms on or near thimbles 10 ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] A. Single Thimble Methods 10Contraction Algorithm 11HMC on thimbles 12Langevin on thimbles 13Case study: bosonic gases 14B. Generalized thimble method 15Case study: 0+1D Thirring model 16Case study: 1+1D Thirring model 18C. Trapping and tempered algorithms 19D. Algorithms for the Jacobian 20Case study: real time field theory 22E. Gauge theories 24Case study: Heavy-Dense QCD 24Case study: 2D QED 25IV. Other manifolds and the algorithms that can find them 26A. Well beyond thimbles 26B. Learnifolds 26Case study: 1+1D Thirring model revisited 28C. Path optimization 29Case study: 2+1D Thirring model 29V. Conclusion and prospects 30 a r X i v : . [ h e p - l a t ] J u l Acknowledgments 31A. Computation of the Jacobian 31B. Another definition for thimbles 32References 32
I. THE SIGN PROBLEM
Monte-Carlo methods have been used with great successto study problems ranging from classical systems of particlesto studies of hadrons using lattice quantum chromodynam-ics. The usual setup is to a formulate the problem – classicalor quantum – in a way analogous to a classical statisticalsystem. Observables are then given by a multidimensionalintegrals involving a Boltzmann factor which is computednumerically by importance sampling. There are, however,important systems of great interest that cannot yet be solvedusing standard Monte-Carlo methods. These are the systemswhere the statistical weights become either complex or whosesigns oscillate. Roughly speaking, we say that the system suf-fers from a sign problem when the phase fluctuations increaseas the size of the system is increased. These fluctuations leadto delicate cancellations that preclude a stochastic evaluationof the integral. This occurs in the study of neutron matterfound in neutron stars, the repulsive Hubbard model awayfrom half-filling and all field theoretical/many-body observ-ables in real time. Not surprisingly, solving the sign prob-lem is of central importance in many fields of Physics anda number of approaches have been proposed to either solveor alleviate this problem. Some are more generic, and someare problem specific, but all approaches have fallen shortof meaningfully addressing the Physics of the systems men-tioned above.In this review we will focus on a novel set of related meth-ods relying on the analytical properties of the configurationweights. The fundamental idea is to express the partitionsum as an integral over real degrees of freedom and com-plexify each variable. The partition sum is originally an in-tegral over the real manifold in this enlarged configurationspace, however, as we will discuss, we can deform the multidi-mensional integration contour—without changing the valueof the partition function—to a manifold that has better nu-merical properties. In particular, the phase fluctuations areeither eliminated or significantly reduced. We will describethe geometry of the complex field space, its critical points,and the algorithms used to both find suitable manifolds andto integrate over them. All of these steps will be exemplifiedin simple field theories, usually in lower number of dimen-sions, that contain, however, all the properties of the theoriesof the greater physical interest.
A. Field Theory/Many-Body Physics as a path integral
The expectation value of any observable O in field theorycan be calculated by the path integral (cid:104) O (cid:105) = 1 Z (cid:90) Dφ e − S E ( φ ) O ( φ ) , Z = (cid:90) Dφ e − S E ( φ ) . (1.1)Here, φ is the generic name of the fields in the theory and S E is the euclidean (imaginary-time) action evaluated overa euclidean “time” β , equal to the inverse temperature ofthe system . The path integral in Eq. 1.1 is an integralover an infinite dimensional space. In order to evaluate itnumerically (and to properly define it), we consider a dis-cretized version where spacetime is replaced by a finite lat-tice. After discretization, the path integral becomes a fi-nite dimensional integral, albeit one over a very large num-ber of dimensions, proportional to the number of spacetimepoints composing the lattice. This is equivalent to a classi-cal statistical mechanics problem in four spatial dimensions,where the state of the system is described by the field φ defined on the entire four dimensional grid, and the prob-ability of each state is controlled by the Boltzmann factorexp[ − S E ( φ )]. Using Monte-Carlo methods, a set of n con-figurations { φ (1) , . . . , φ ( n ) } is generated with the probabilitydistribution exp[ − S E ( φ )] /Z . The observables and their er-rors are then estimated using (cid:104) O (cid:105) = 1 n (cid:88) a O ( φ ( a ) ) , (cid:15) O = (cid:115) n ( n − (cid:88) a [ O ( φ ( a ) ) − (cid:104) O (cid:105) ] . (1.2)Numerous algorithms have been developed to obtain config-urations φ ( a ) distributed according to e − S E [ φ ] in an efficientway. The cost of the sampling process increases with a mod-erate power of the spacetime volume V (between 1 and 2),despite the fact that the Hilbert space dimension of the cor-responding quantum system grows exponentially with thespace volume. This is the great advantage of Monte Carlomethods over direct diagonalization procedures. Similar expressions are obtained for the partition function Z =tr e − βH of non-relativistic quantum systems by discretizing bothspace and time, then using the Trotter formula. There is no assumption that the theory is relativistic. In fact, non-relativistic systems in the second quantized form are frequently stud-ied within this formalism.
B. Physical systems with sign problems
Many theories of interest in theoretical physics have signproblems in all currently known formulations. In fact, sys-tems that cannot be fully understood because a sign problemhinders the use of Monte Carlo simulations are pervasive inall subfields of Physics (and Chemistry). Among those somehave become “holy grails” in their respective field, problemswhose solutions would have a revolutionary impact.For instance, in nuclear physics, QCD at finite baryon den-sity has a sign problem. This prevents the understandingfrom first principles of both neutron stars and supernovae.Extensive work has been expended to evade this sign prob-lem (see, for isntance, the following reviews and the refer-ences therein (Aarts, 2016; de Forcrand, 2010; Karsch, 2000;Muroya et al. , 2003; Philipsen, 2007). Quantum Monte Carlo(QMC) studies of nuclei using “realistic nucleon-nucleon in-teractions” also suffer from the sign problem (Carlson et al. ,2015; L¨ahde et al. , 2015; Wiringa et al. , 2000). The “con-strained path algorithm” (Zhang et al. , 1995, 1997) is awidely-used approximate method to address these sign prob-lems . Lattice Field Theory studies of nuclei have simi-lar behavior; sign problems appear in studies of nuclei withdifferent proton and neutron numbers, and when repulsiveforces become sufficiently large (Elhatisari et al. , 2017; Epel-baum et al. , 2014; Lee, 2009; Lee et al. , 2004). Furthermore,lattice and QMC studies of nuclear matter encountered inastrophysics suffer from the sign problem. This includesspin polarized neutron matter (Fantoni et al. , 2001; Gan-dolfi et al. , 2014; Gezerlis, 2011) and lattice EFT studies ofnuclear matter beyond leading order (Lu et al. , 2019a) .Many cold atom systems, when formulated with latticeor QMC methods, exhibit sign problems as well. Both spinand mass imbalanced spin 1/2 fermions have a sign problem(Braun et al. , 2013; Roscher et al. , 2014). This sign-problemmakes it prohibitively difficult, for example, to conclusivelydemonstrate the existence of a number of conjectured phases(like the “LOFF” phases) in more than 1 + 1 dimensions.Bosonic non-relativistic systems exhibit sign problems aswell. This includes bosons under rotation (Berger et al. ,2020) and coupled to spin-orbit interactions (Attanasio andDrut, 2020). For a review see (Berger et al. , 2019). The constrained path algorithm is a generalization of the “fixed-nodeapproximation”, a similar approximate technique for avoiding thesign problem (Anderson, 1975) Unpolarized neutron matter, however, can be formulated free of thesign problem (Chen and Kaplan, 2004; Lee and Sch¨afer, 2005) ”Wigner SU (4)” symmetric approximations to pionless EFT have nophase oscillations and have been profitably used (Lee, 2007; Lu et al. ,2019b; Wigner, 1937). A wide variety of lattice-supersymmetric models sufferfrom a sign problem too (for a review see Ref. (Schaich,2019)). In particular, first-principles tests of the gauge-gravity duality conjecture, even in the simplest case of re-produce supergravity black hole thermodynamics from D0-brane quantum mechanics, can only claim to be bona-fidecontrolled tests of the duality if the phase fluctuations areunder control (Berkowitz et al. , 2016; Hanada et al. , 2011).Sign problems are found in condensed matter physics aswell. A particularly well-known example is the Hubbardmodel away from half filling (Hubbard, 1963; Loh et al. , 1990;White et al. , 1989) thought to model essential characteris-tics of high T c superconductors. Path integral formulationsof fullerene exhibit the sign problem as well (Ostmeyer et al. ,2020). Furthermore, some models of frustrated magnetismon triangular and kagom´e lattices, of interest for their con-jectured spin-liquid ground states, exhibit the sign problem(Lacroix et al. , 2011; Sindzingre et al. , 1994). As a resultthere is uncertainty in the zero-temperature properties ofthese models. C. Reweighting and the sign problem
The standard workaround for sampling complex actions isto use reweighting . The idea is to split the integrand into apositive part that is used for Monte-Carlo sampling, usuallythe absolute value of the integrand, and a fluctuating partthat is included in observables. Using the absolute value asa sampling weight, we have the following identity (cid:104) O (cid:105) = (cid:104) Oe − i Im S E ( φ ) (cid:105) (cid:104) e − i Im S E ( φ ) (cid:105) , (cid:104) O (cid:105) = (cid:90) Dφ e − Re S E ( φ ) Z O ( φ )(1.3)and Z ≡ (cid:82) Dφ e − Re S E ( φ ) . The idea, then, is to use the phase quenched action Re S E to sample configurations, andtake into account the imaginary part of the action when com-puting observables. From a numerical point of view, this pro-cedure works when the phase fluctuations are mild and wecan estimate the phase average, (cid:104) e − i Im S E ( φ ) (cid:105) , with enoughaccuracy; this means that the error estimate for this averageshould be significantly smaller that its mean. Since the mag-nitude of the phase for each configuration is one, to resolvethe mean accurately we require a number of configurations n (cid:29) / (cid:104) e − i Im S E ( φ ) (cid:105) . When the average phase is very small,reweighting requires a very large number of samples and be-comes impractical. For many systems at finite density, thephase average goes to zero exponentially fast in the spatialvolume/inverse temperature. This is because the phase av-erage is the ratio of two partition functions: (cid:104) e − i Im S E ( φ ) (cid:105) = ZZ = e − βfV e − βf V = e − βV ∆ f , (1.4)where ∆ f = f − f > phase quenched system. In this case, the numerical effort grows exponen-tially as we increase the volume and/or lower the tempera-ture. This is what is usually defined to be the sign problem .An even worse problem arises when calculating real time cor-relation functions. In that case, we are interested in integralsof the form (cid:104)O(cid:105) = (cid:90) Dφ e iS ( φ ) O , (1.5)where S is the real time (Minkowski space) action of the sys-tem . Since there is no damping of the magnitude of the in-tegrand and the value of the field φ ( t, x ) (for any t, x ) grows,the average phase is strictly zero, even for small sized sys-tems. A similar argument applies to observables, like partondistribution functions, defined on the light cone.We should note that, the existence of a sign problem doesnot necessarily preclude numerical study. There are caseswhere the sign problem is mild enough that most relevantinformation about the system in the region of interest canbe extracted before the sign fluctuations become an obsta-cle. For example, when studying the phase diagram of asimple heavy-dense quark model for QCD (see below), theendpoint of the first order phase transition can be studiedvia reweighting for system sizes as large as 100 even toughthe model has a sign problem (Alford et al. , 2001). We men-tion this study to point out that, from a practical point ofview, methods which merely reduce sign fluctuations, with-out completely eliminating them, are also important. D. The absence of a general solution
It is of theoretical, if not practical, interest to know if ageneric solution to the sign problem exists. If one takes anexponentially vanishing average sign in the system size as thedefinition of the sign problem, then there are definitely mod-els in which the sign problem can be solved. For instance,for many systems it is possible to rewrite the path integral In thermal equilibrium at non-zero temperature, real time corrella-tors can be computed from path integrals defined in the closed-timecontour in complex time (Keldysh, 1964; Schwinger, 1961). See sec-tion III.D. using a different set of states and obtain an expression free ofphase fluctuations. This was accomplished, for example, forthe two-component scalar theory using dual variables (En-dres, 2007; Gattringer and Kloiber, 2013), and by reorganiz-ing the summation over configurations for the heavy-densesystem mentioned earlier (Alexandru et al. , 2018d; Alford et al. , 2001). Similarly, there is a class of fermionic mod-els that, when formulated in terms of fermion bags (Alford et al. , 2001; Ayyar et al. , 2018; Chandrasekharan, 2012, 2013;Chandrasekharan and Wiese, 1999; Hann et al. , 2017; Huff-man and Chandrasekharan, 2016, 2020, 2014) have strictlypositive Boltzmann weights even though other formulationshave a severe sign problem. As it turns out, a solution ofthis kind is unlikely to work for all systems.There is an often-cited, general argument implying that ageneric solution to the sign problem, applicable to all sys-tems, is extremely unlikely to exist. It relies on the NP (cid:54) =Pconjecture from computational theory. NP decision prob-lems are problems that can be solved on a non-deterministic
Turing machine in a time that increases only polynomiallywith the system size, whereas P problems are the ones thatcan be solved in polynomial time in a deterministic way.While no proof exists, it is widely believed that there are NPproblems that are not P. In connection to this question, animportant subset of NP problems are the NP-hard or NP-complete problems. If any of these NP-hard problems can besolved in polynomial time on a classical computer, then allNP problems can, invalidating the conjecture. There are spinglass-like systems with a sign problem that can be mappedinto NP-hard problems (Troyer and Wiese, 2005). Using thechain of arguments above, a generic solution to the sign prob-lem that would solve this problem, would imply
N P = P ,which is considered highly unlikely. E. A brief survey of methods to deal with sign problems
As mentioned above, some of the most physically interest-ing models in particle, nuclear and condensed matter physicshave sign problems. Given the interest in these problems, itis not surprising that a variety of approaches have been triedto either solve or circumvent the sign problem. In this reviewwe will focus on Lefschetz thimble inspired methods, but wewant to point out some approaches attempted through theyears to understand the phase diagram of QCD and otherrelativistic theories.A first set of methods uses simulations in the parameterregion where the action is real; the result is then extrapolatedin the region of interest. One version of this idea is to relyon results from imaginary chemical potential. Monte Carlosimulations can be used either directly to infer features of thephase diagram for real chemical potential or to compute ob-servables and fit them using a polynomial ansatz or a Pad´eapproximations and then analytically continue these func-tions to real values of µ (Bellwied et al. , 2015; Bonati et al. ,2015; Borsanyi et al. , 2020; Cea et al. , 2014; D’Elia and Lom-bardo, 2003, 2004; de Forcrand and Philipsen, 2002, 2003).Another approach is to compute the derivatives of thermo-dynamic observables with respect to µ at µ = 0, then useTaylor expansions to extend these results to µ > et al. , 2019; Bonati et al. , 2018; Endrodi et al. , 2011; de For-crand et al. , 2000; Kaczmarek et al. , 2011; Miyamura, 2002).Yet another method is to use multiparameter reweighting bycombining simulations from different temperatures at µ = 0to determine the phase transition line and critical point inQCD (Fodor and Katz, 2002).Another class of methods attempt to alleviate the signproblem by a clever rewriting the path integral in terms ofnew variables. One possibility is to reorganize the sum overthe configurations in subsets that have either only positivesign contributions to the partition function, thus solving thesign problem, or a much reduced sign problem (Alexandru et al. , 2018d; Alford et al. , 2001; Bloch et al. , 2013; Chan-drasekharan and Wiese, 1999; Karsch and Mutter, 1989;Rossi and Wolff, 1984). Another direction is to reformu-late the problem in terms of dual variables in which the signproblem is absent (Endres, 2007; Gattringer and Kloiber,2013). It turns that for QCD, the use of the canonical en-semble partition function (as opposed to the grand canonicalensemble) makes the sign fluctuations milder and it can beused to investigate small enough systems (Alexandru et al. ,2005; Alexandru and Wenger, 2011; Barbour et al. , 1988;de Forcrand and Kratochvila, 2006; Hasenfratz and Tous-saint, 1992; Kratochvila and de Forcrand, 2006; Li et al. ,2011, 2010). Finally, Fermi bags are enough to completelyeliminate the sign problem in some low dimensional mod-els (Alford et al. , 2001; Ayyar et al. , 2018; Chandrasekha-ran, 2012, 2013; Chandrasekharan and Wiese, 1999; Hann et al. , 2017; Huffman and Chandrasekharan, 2016, 2020,2014). These methods are very model dependent and requireinsight to be applied in each new class of models.Recently a proposal based on the density of states methodwas explored as a way to alleviate sign fluctuations (Fodor et al. , 2007; Garron and Langfeld, 2016, 2017; Gattringer andT¨orek, 2015; Langfeld and Lucini, 2014).Finally, there is a significant effort to simulate QCD at fi-nite density using the complex Langevin approach (Klauder,1983; Parisi, 1983) , based on the idea of stochastic quan- See (Berger et al. , 2019) for a recent review of complex Langevinapproach. tization (Parisi and Wu, 1981). This method shares withthe thimble methods its starting point: the configurationspace of N real degrees of freedom is extended to a N di-mensional complex one. The important difference is thatcomplex Langevin approach sets up a stochastic process thatmoves freely in this enlarged space of 2 N real degrees of free-dom, whereas the methods we discuss in this review samplean N dimensional manifold. Results show that, while insta-bilities are present in complex Langevin QCD simulations,for heavy quark masses credible results can be obtained fortemperatures above the deconfinement transition. In thehadronic phase, the simulations become unstable and un-reliable (Aarts, 2009; Aarts et al. , 2013, 2011, 2010; Aartsand Stamatescu, 2008; Fodor et al. , 2015; Seiler et al. , 2013;Sexty, 2014). II. CAUCHY THEOREM, HOMOLOGY CLASSES ANDHOLOMORPHIC FLOWA. Deformation of domain of integration: a multidimensionalCauchy theorem
The well known Cauchy theorem for functions of one com-plex variable states that for an analytic function f ( z ) theintegral over a closed loop vanishes: (cid:73) C f ( z ) = 0 . (2.1)This can be used to “deform” the contour of integration from,say, the real line, to a different contour on the complex plane,as long as the initial and final points of the contours coin-cide. In many applications the contour starts and/or ends ata point on the infinity and the issue becomes whether mov-ing these ending points may cross a “singularity of f ( z ) atinfinity”. For instance, take the integral (cid:90) dφ e − φ (2.2)over different contours on the complex plane starting/endingat different points at the infinity. Since there are no singu-larities at any finite values of z , Cauchy’s theorem allows usto deform the contour of integration as long as no singularity“at infinity” is crossed. The integral in Eq. 2.2 is well-defined(it converges) if and only if the initial and final asymptoticdirections of the contour are in the regions A, . . . , D shownin Fig. 1. The integral over two different contours whoseends lie on the same regions have, on account of Cauchy’stheorem, the same value. For instance, the real line, contour1, is equivalent to contour 2 since both start in region A andend in region B . The integral over contour 3 is not even
14 32A BCD - - - - FIG. 1: Several contours of integration for the integral inEq. 2.2. Contour 1 (the real line) and 2 produce the sameresult. Contour 4 a different result while the integral overcontour 3 is divergent. The gray areas show direction in thecomplex plane (“good” regions) where the integrandvanishes fast enough so the integral converges.well-defined as it diverges, while the value for the integralover contour 4 is different from the value on contours 1 or 2.In fact, imagine starting from the real line and continuouslydeforming it towards contour 4. At some point the integralwill cease to be well-defined as its end point leaves region B and the integral becomes divergent. As the end point en-ters region C the integral becomes finite again but acquiresa different value than on the real line.In fact, there are only three independent classes of con-tours (known as “homology classes”) on which the integralin Eq. 2.2 may be evaluated: those that start in region A and end in region B , C or D , denoted A → B , A → C and A → D , respectively. Any other contour with different aasymptotic behavior, for instance B → C , can be expressedas a linear combination of contours (with integer coefficients)belonging to one of these three classes. Cauchy’s theoremguarantees that any contour that lies in one of these classescan be smoothly deformed to some other contour in the sameclass without changing the value of the integral. In contrast,as explained above, it cannot be deformed to a contour thatlies in a different class. In short, all possible domains overwhich the integral Eq. 2.2 is well-defined can be classified asa linear combinations of three discrete classes of contours.Each class contains a continuous family of “equivalent” con-tours that can be smoothly deformed to one another withoutchanging the value of the integral. As we will see below, thereason that there are three classes is that the function φ inthe exponent is a quartic polynomial which in general has FIG. 2: Above is a schematic of a multi-dimensionaldeformation. The original domain of integration, M ⊂ R N , is deformed to M ⊂ C N . This deformationsweeps out a manifold B ⊂ C N whose boundary is ∂ B = −M ∪ M .three saddle points. All the observations above generalize to higher dimensions.Instead of integrals over one dimensional paths we will con-sider integrals over N -cycles, orientable manifolds with noboundary with real dimension N immersed in the 2 N di-mensional space. The integral over a cycle M is defined by (cid:90) M f ( φ ) dφ ∧ · · · ∧ dφ N = (cid:90) M f ( φ ( ζ )) det J ( ζ ) dζ . . . dζ N , (2.3)where φ i = Φ i ( ζ , . . . , ζ N ) is a parametrization of the N -dimensional manifold M by N real coordinates ζ , . . . ζ N , M is the region of R N used to parametrize M anddet J ( ζ ) = ∂ ( φ ...φ N ) ∂ ( ζ ...ζ N ) is the determinant of the Jacobian ofthe parametrization, which is in general a complex number. φ stands for all φ , . . . , φ N (and similarly for ζ ).Assume that we have two such cycles M and M that canbe smoothly deformed into one another. The space swept bythe deformation will be denoted with B and the two cyclesform the boundary ∂ B = M − M where the minus signmeans oriented in opposite way (see Fig. 2). By Stokes’ Note our simple example actually a degenerate case where all threesaddle points are at φ = 0, but it is easy to lift the degeneracy byadding a term (cid:15)φ to the exponent. theorem we have : (cid:90) ∂ B f ( φ ) dφ ∧ · · · ∧ dφ N = (cid:90) B df ( φ ) ∧ dφ ∧ · · · ∧ dφ N , (2.4)where df = ∂f∂φ i dφ i + ∂f∂ ¯ φ i d ¯ φ i ( ¯ φ is the complex conjugate of z ).Since f ( z ) is assumed to be holomorphic we have ∂f∂ ¯ φ i = 0. Inthe sum ∂f∂φ dφ + . . . ∂f∂φ N dφ N every term is proportional toone of the terms in dφ ∧ · · · ∧ dφ N so df ∧ dφ ∧ · · · ∧ dφ N = 0since dφ i ∧ dφ i = 0. We arrive then at (cid:90) ∂ B f ( φ ) dφ ∧ · · · ∧ dφ N == (cid:90) M − M f (Φ( ζ )) det J ( ζ ) dζ . . . dζ N = 0 . (2.5)which is the generalization of the Cauchy theorem we areinterested in . This theorem can be used to deform themanifold of integration without altering the value of the in-tegral just as we discussed above for the one dimensionalcase. In fact, our discussion of contour deformation readilygeneralizes to the multidimensional case. For manifolds ap-proaching the infinity along certain “directions” (in reality, N -dimensional planes) the integral is convergent and well-defined (“good regions”); for others it is not. Furthermore,it can be shown, assuming the integrand is well behaved in asense discussed below, that the manifolds for which the inte-gral converges are separated in discrete equivalence classes:those with the same asymptotic properties lead to the sameintegral. A continuous deformation of manifolds of integra-tion from one equivalent class to another, that is, from one“good region” to another necessarily goes through manifoldswhere the integral diverges. Such deformations are the ana-logue of deformations crossing a “singularity at the infinity”in the one-dimensional case. All this is in close analogy tothe familiar one-dimensional case. A detailed discussion ofthe mathematical details can be found in (Pham, 1983). B. Holomorphic gradient flow
We will be interested in deforming integrals from R N (thereal cycle) to some other N -cycle without altering the valueof the integral but alleviating the sign problem in integrals Readers not familiar with the formalism of differential forms maytake the right side of Eq. 2.3 as the definition of an integral over N -dimensional manifolds embedded in C N . We will use this definitionextensively in this paper. We thank Scott Lawrence for a discussion on this point. of interest in field theory, which are typically of the form (cid:90) R N e − S ( φ ) O ( φ ) dφ i , (2.6)where S is the action of the theory and O some observable.One way of performing this deformation is with the help ofthe holomorphic flow. The holomorphic flow is defined forevery action S by the differential equations: dφ i dt = ∂S∂φ i . (2.7)For every point φ in R N and a fixed flow time T , the solutionof Eq. 2.7 with the initial condition φ ( t = 0) = ζ defines apoint ˜ φ = F T ( ζ ) in C N . By flowing all points of R N in thismanner we obtain the flowed manifold M T = F T ( R N ) .The holomorphic flow has two important properties: ddt S R = 12 (cid:20) dSdt + dSdt (cid:21) = ∂S∂φ i ∂S∂φ i | ≥ , (2.8) ddt S I = 12 i (cid:20) dSdt − dSdt (cid:21) = 12 i (cid:20) ∂S∂φ i ∂S∂φ i − ∂S∂φ i ∂S∂φ i (cid:21) = 0 , that is, the imaginary part S I is constant along the flow whilethe real part of the action S R increases monotonically (thatis why Eq. 2.7 is also called upward flow). The fact that S R increases along the flow means that the integrand van-ishes along asymptotic directions even faster in the flowedmanifold M T than in R N , leading to the convergence of theintegral at all T . By the arguments exposed above, thismeans that M T is equivalent to R N for the purpose of com-puting the integral, that is, it is in the same homology classas R N , as in the one dimensional example explained in thebeginning of this section. C. Lefschetz Thimbles and Picard-Lefschetz theory
Even though M T is equivalent to R N , evaluating the pathintegral on M T rather than R N is computationally advan-tageous in controlling the sign problem. Before we explainwhy this is, we first introduce the necessary mathematicalbackground (for a different perspective, see Appendix B).We begin by focusing on the stationary points of the flow,namely the critical points of the action φ c where ∂S/∂φ i | φ c = Other flows to generate manifolds were proposed in (Tanizaki et al. ,2017). This can also be seen by noting that the holomorphic flow is thegradient flow of S R and the hamiltonian flow for the “hamiltonian” S I .
0. The
Lefschetz thimble T attached to a critical point φ c isdefined as the set of initial conditions φ (0) ∈ C N for whichthe downward flow dφ i dt = − ∂S∂φ i . (2.9)asymptotically approaches the critical point. Similarly, the dual-thimble K is the set of all point for which the upwardflow asymptotes to φ c . For a constructive definition for T ,we begin by linearizing the flow around φ c : dφ i dt = ∂ S∂φ i ∂φ j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) φ = φ c (cid:124) (cid:123)(cid:122) (cid:125) H ij ( ¯ φ j − ¯ φ cj ) (2.10)whose solution can be written as φ ( t ) − φ c = N (cid:88) a =1 c a ρ ( a ) e λ a t , (2.11)where c a are real and ρ ( a ) are the solutions to the modifiedeigenvector problem (“Takagi vectors” (Takagi, 1924)) H ij ρ ( a ) j = λ a ¯ ρ ( a ) i . (2.12)The modified eigenvalues λ a can be chosen to be real,and then the eigenvalues/eigenvectors come in pairs( λ a , ρ ( a ) ) , ( − λ a , iρ ( a ) ). The set of N vectors ρ ( a ) which de-fine the directions around a critical point where the flowmoves away from the critical point forms a basis (with realcoefficients) for the tangent space of T at φ c . Likewise theset of N vectors iρ ( a ) which define the directions arounda critical point where the flow moves towards the criticalpoint forms a basis for the tangent space of K at φ c . Thesetwo tangent spaces together span the tangent space of C N .With this knowledge, in the infinitesimal neighborhood ofthe critical point, we can solve for the vanishing cycle v ( (cid:15) )as S ( φ ) − S ( φ c ) ≈ z i H ij z j = (cid:15) which is an N − T . The thimble can beconstructed by taking the vanishing cycle as the initial con-dition and flowing by upward flow: T = ∪ ≤ T < ∞ F T ( v ( (cid:15) )),when (cid:15) →
0. In other words we can build the thimble sliceby slice by using the flow. We can further use the fact thatthe flow defines a one-to-one map between the initial pointand the flowed point and instead consider an infinitesimallysmall N dimensional ball, B , in the tangent plane. B is al-ready a small portion of the thimble near φ c . If we take B as the initial condition, its image under upward flow with T → ∞ is the thimble: T = F T →∞ ( B ). This is the mainidea behind the “contraction algorithm” that is a method to - - ( z )- ( z ) FIG. 3: Thimbles (blue), dual thimbles (yellow), criticalpoints (blue dots), their pre-image under the flow (orangestars), and the flowed real line (dashed red) for G = 1 . e i . , p = 1 , µ = 0 . m = i .
1. The arrowsindicate the direction of the upward flow.simulate path integrals on a given thimble (see section III.A).For a concrete illustration of these ideas, considerFig. 3, where the action is taken to be S ( φ ) = φ /G − log (cid:2) ( p + iµ ) + ( φ + m ) (cid:3) . S can be thought of as a toymodel for the action of a fermionic model coupled to an aux-iliary field φ , after the fermions have been integrated out.Notice that e − S is a holomorphic function, even though S isnot; this is a feature common to theories with fermions. Thistheory has three critical points, attached to which are thim-bles and dual-thimbles. Only thimbles 1 and 2 contributeto the integral. The real line, evolved by the holomorphicflow by a time T = 1 . N -cycles identified by their asymptoticbehavior. The thimbles are representatives of these equiva-lence classes, each thimble representing a different class .More concretely, let us assume that there are finitely manycritical points, φ cα indexed by α and Im S ( φ cα ) (cid:54) = Im S ( φ cβ )for α (cid:54) = β . Attached to each critical point there exists a In this review we only consider integration domains with no bound-aries. The generalization of thimbles with boundaries are studiedextensively in (Delabaere and Howls, 2002). These assumptions ensures that no two critical point is connected byflow since the flow conserves the imaginary part which is known asthe
Stokes phenomenon . We will discuss Stokes phenomenon brieflyin chapter III. thimble, T α , and a dual thimble K α . As explained above,different thimbles do not intersect each other (they carrydifferent values of Im S ) and T α intersects K β if and only if α = β . In other words (cid:104)K α , T β (cid:105) = δ αβ where (cid:104) , (cid:105) denotes theintersection number between two cycles. The intersectionoccurs at φ cα . Since Re S is bounded from below on a thim-ble, the integral (2.3) is guaranteed to be well-defined whenevaluated on a thimble T α . In fact, the set of all thimblesforms a complete basis for the space of equivalence classes of“good domains” (i.e. the homology group) and any domain,say M , over which (2.3) is well-defined is equivalent to aunique linear combination of thimbles (Pham, 1983): M ≡ (cid:88) α n α ( M ) T α , n α ( M ) = (cid:104)K α , M(cid:105) . (2.13)Here the integer coefficients n α are given by the number ofintersections between M and the dual thimble K α . The signdepends on the relative orientations of K α and M . Noticethat some of the n α may vanish; it is said then that thoseparticular thimbles do not contribute to the integral. A sim-ple example of this is shown in Fig. 3.Thimbles are the multi-dimensional generalization of theconcept of “steepest descent” or “stationary phase” contourfrom the theory of complex functions of one variable. Nat-urally, they are useful in studying the semi-classical expan-sion of path integrals in field theory (Cherman et al. , 2014;Dunne and ¨Unsal, 2016) and their asymptotic analysis (see(Aniceto et al. , 2019) for a recent review of the new devel-opments related to “resurgent transseries”). Also, thimbleshave been used in attempts at defining ill defined path inte-grals by defining the relevant partition function as an inte-gral over one or more thimbles instead of over R N (Harlow et al. , 2011; Witten, 2010, 2011). For our purposes, the rel-evant property of the thimbles is that the imaginary part ofthe action and, consequently, the phase of the integrand ofthe partition function, is constant on the thimble. There-fore instead of evaluating the path integral on R N wherethe phase is a rapidly oscillating function, evaluating it inon the equivalent thimble decomposition where the phaseis piecewise constant can provide significant practical advan-tage. This fact by itself, however, is not quite enough to solvethe sign problem. As can be seen from Eq. 2.3, the phase ofthe integrand depends also on the phase of the Jacobian (the“residual phase”). The Jacobian will have a rapidly oscillat-ing phase if the shape of the manifold of integration oscillatesquickly along real and imaginary directions. For theories inthe semi-classical regime this does not happen because theparts of the thimble with significant statistical weight areclose to the critical point. Experience shows that the resid-ual phase in many strongly coupled models introduces a very mild sign problem (see below for many examples) .An important question that naturally arises then is: whichthimble, or combination of thimbles, is equivalent to the R N ?We can answer this question by considering the manifold M T obtained by taking every point of R N as an initial conditionand flowing them by a “time” T . Since the real part of theaction grows monotonically with T the integral remains con-vergent at all T and, by the arguments above, the value of theintegral remains the same. Since R N and the dual thimble ofany critical point are N dimensional spaces they will generi-cally intersect on isolated points, if they intersect at all. If wecall each of those points ζ c we have φ c = F T →∞ ( ζ c ). Start-ing from one of these intersection points ζ c the flow leads tothe critical point on a trajectory lying on the dual thimble K (see Fig. 3 and Fig. 4). The trajectory starting at pointsnear ζ c initially approaches the critical point but then veersalong the unstable directions of the critical point slowly ap-proaching the thimble (see Fig. 4). Points in R N far fromthe intersection points take a more direct route towards infin-ity (or some other point where the action diverges). There-fore, all points in R N flow, at large times, to points near aset of thimbles that, together, are equivalent to R N (or topoints where the action diverges). Furthermore every thim-ble is counted as many times as there are intersection pointsbetween the corresponding dual thimble. Consequently thethimble decomposition of R N can explicitly be obtained asthe limit, M T →∞ = (cid:88) α n α ( R N ) T α , where M T =0 = R N . (2.14)It is worth stressing that even though the thimble decompo-sition is obtained as the infinite flow time limit, the value ofthe integral remains unchanged during the deformation and M T is equivalent to R N for any finite value of T : Z = (cid:90) R N dφ e − S ( φ ) = (cid:90) M T dφ e − S ( φ ) = (cid:88) α n α ( R N ) (cid:90) T α dφ e − S ( φ ) (2.15)It should be noted that in theories where more than onethimble contribute to the partition function, there is a possi-bility that the contributions from different thimbles comewith phases exp( − i Im S eff ) (constant over each separatethimble) which induces a sign problem. This kind of sign One can construct examples of extremely strongly coupled theorieswhere the residual phase introduces a severe sign problem (Lawrence,2020) ζ c ∈ R N flow to the critical points φ c ∈ C N . The points in the neighborhood of each ζ c approach the thimble but eventually veer off. In this figurewe show two such thimbles T α and T β for illustrativepurposes.problem is not helped by integrating over thimbles. How-ever, in order for the contributions from different thimblesto (nearly) cancel an (approximate) symmetry is requiredrelating the contribution of different thimbles. Monte Carlomethods can be adapted to situations like that by samplingpoints related by the symmetry at the same time.In field theories, where the dimensionality of the integralis large, it is extremely difficult to find the thimbles – it isin fact equivalent to classifying all complex solutions of theequations of motion– and even harder to find their intersec-tion numbers n α . The discussion of the previous paragraphwill be useful however, in establishing an algorithm to solvethis problem numerically and “on-the-fly” during a MonteCarlo run. It also clarifies the fact that there is nothing spe-cial about thimbles as opposed to other manifolds obtainedfrom flowing R N by a finite time T . These other manifoldsdo not improve the sign problem as much as the thimbles dobut still give the correct result for the integral and can beadvantageous for numerical/algorithmic reasons. III. ALGORITHMS ON OR NEAR THIMBLESA. Single Thimble Methods
Early simulations using complex manifolds focused onsampling the path integral contribution from the “main”thimble, the thimble associated to the critical point withthe smallest value of S R ( φ c ). This was based on the hopethat in the relevant continuum/thermodynamic limits thepath integral would be dominated by the contribution ofa single thimble or that a regularization can be defined for relevant QFTs in term of a single thimble path inte-gral (Cristoforetti et al. , 2012a; Di Renzo et al. , 2019). Al-though there is no evidence that this conjecture is valid, algo-rithms to sample a single thimble are obvious stepping stonestowards multi-thimble integration. We will discuss in thissection the algorithms proposed to sample the integral alonga single thimble: the contraction algorithm, a Metropolisbased algorithm (Alexandru et al. , 2015), a Hybrid Monte-Carlo algorithm (Fujii et al. , 2013), and the Langevin algo-rithm (Cristoforetti et al. , 2012a).As discussed earlier finding the thimble decomposition forthe path integral is a very hard problem which was only at-tempted for quantum mechanical systems (Fujii et al. , 2015).However, in many cases it is feasible to find the “main”thimble even for realistic systems using the symmetry of theproblem. The problem of finding the critical point is usu-ally reduced to a “gap” equation to be solved analyticallyor numerically. For the algorithms discussed in this section,we assume that we have identified this critical point and wewant to sample configurations on the corresponding thimble.Another important challenge facing any algorithm for theMonte Carlo evaluation of integrals over thimbles is to re-strict sampling to the thimble manifold. For most systemsthere is no known method that can identify points on thethimble based on the local behavior of the action. Rather, apoint has to be transported though the reverse flow (Eq. 2.9)to decide whether it approaches the critical point or not. Thethimble attached to this point can then be constructed byintegrating the upward flow equations starting in the neigh-borhood of the critical point. As the thimble on the neigh-borhood of the critical point is approximated by the tangentspace spanned by the Takagi vectors with positive eigenval-ues (in Eq. 2.12) we can take points on the tangent plane(close enough to the critical point) as the initial conditionsof the holomorphic flow Eq. 2.7 to find points lying on thethimble. This “backward-and-forward” procedure then al-lows us to find points on the thimble nearby other pointson the thimble, as required by Monte Carlo procedures, atthe expense of integrating the flow equations. This processprovides a map between the N dimensional neighborhood ofthe critical point to the thimble attached to it. It is an es-sential ingredient for all single thimble algorithms discussedhere. For a given parametrization of the tangent space nearthe critical point φ c : φ n = φ c + N (cid:88) a =1 ζ a ρ ( a ) , ζ a ∈ R , (3.1)integrating the upward flow for a time T produces a map φ n → φ f = F T ( φ n ). Here φ n is a point near φ c and φ f ismoved far by the flow. For large enough T , this will map1a small neighborhood of the critical point into a manifoldvery close to the thimble and the larger the value of T , thecloser the manifold generated by the φ n → φ f mapping isto the thimble. As a practical method of determining anappropriate value for T , simulations can be carried out forincreasing values of T until the results converge.Having chosen an appropriate T , we have now the meansto parametrize the thimble using the tangent plane close tothe critical point. We can then approximate the integral overthe thimble as (cid:90) T dφ f e − S ( φ f ) ≈ (cid:90) U dφ n det J ( φ n ) e − S ( φ f ( φ n )) , (3.2)where J ij = ∂ ( φ f ) i /∂ ( φ n ) j is the Jacobian of the map and U is the region around φ c in the tangent plane that is mappedto the manifold approximating the region of the thimble thatdominates the integral. For the special case where the tan-gent plane is in the same homology class as the thimble, theregion U can be extended to the entire tangent plane andthe relation above becomes exact for all flow times T . Forthe case when the tangent plane is not in the same homol-ogy class, the relation only becomes exact in the limit oflarge T . In practice the region U is generated implicitly inthe simulations: we start in the neighborhood of the criticalpoint and the proposed updates move smoothly, or in smalldiscrete steps, through the configuration space and the po-tential barriers force the simulation to stay in the relevantregion. To fix terminology we will refer to the region U inthe tangent plane as the parametrization manifold and theimage under the map F T ( U ) as the integration manifold.The goal of the algorithms presented here is to samplethe integration manifold according to the Boltzmann fac-tor exp( − S ). Since the action and the integration measureare complex, we need to use a modified Boltzmann factorfor sampling. The probability density we will sample corre-sponds to P ( φ f ) | dφ f | = 1 Z e − Re S ( φ f ) | dφ f | , Z ≡ (cid:90) T | dφ f | e − Re S ( φ f ) . (3.3)The final result for observables will have to include the phase (cid:104)O(cid:105) = (cid:104)O e iϕ (cid:105) (cid:104) e iϕ (cid:105) , e iϕ ≡ e − i Im S ( φ f ) dφ f | dφ f | (3.4)Since we are sampling the configurations from a single thim-ble, or from a manifold that is very close to it, the imaginarypart of the action is constant (or nearly so.) The only fluctu-ation come from the residual phase associated with the phaseof the measure dφ f . If we view this as an integral over the parametrization manifold, then the probability measure is P ( φ n ) = 1 Z e − Re S eff ( φ n ) , S eff = S ( φ f ( φ n )) − ln det J ( φ n ) . (3.5)The complex phase in this case is exp( − i Im S eff ) and thefluctuations of this phase are dominated by the Jacobianphase which correspond to the residual phase. Note thatto compute the effective action for a point φ n in theparametrization space, we have to integrate the upward flowdifferential equation with initial condition φ n for a time T to get φ f . Then S ( φ f ) is the action contribution. The othercontribution comes from the Jacobian. As explained in Ap-pendix A the Jacobian matrix can be computed by integrat-ing the matrix differential equation dJdt = H ( φ ( t )) J ( t ) , (3.6)where H ( φ ( t )) is the Hessian of S along the flow and theinitial condition J (0) is a matrix whose columns form an or-thonormal basis in the tangent to the parametrization spaceat φ n . This equation flows a basis in the tangent space at φ n to a basis in the tangent space at φ f . Since our parametriza-tion space is a hyperplane the basis for the tangent space at φ n can be chosen to be the same at all points in U , for exam-ple the positive Takagi vectors or any other basis spanningthis tangent space.This equation can also be used to map a single infinitesimaldisplacement represented by a vector v n in the tangent spaceat φ n to a displacement represented by a vector v f in thetangent space on the thimble at φ f . In the equation above J ( t ) is then replaced with v ( t ) the column vector representingthe displacement. The initial condition is v (0) = v n and thefinal result, v ( T ) = v f , is a vector in the tangent space at φ f . Because of this we will sometime call this equation the vector flow . Contraction Algorithm
Several sampling algorithms are based on the mapping be-tween the tangent plane and the (approximate) thimble. Themost straightforward is the contraction algorithm (Alexandru et al. , 2015, 2016a), which is generates configurations in theparametrization manifold based on the probability P usingthe Metropolis method (Metropolis et al. , 1953) based on theeffective action Re S eff . The basic process is detailed below.1. After a critical φ c point is identified, the tangent spaceof its thimble is compute by solving Eq. 2.12 and find-ing the ρ ( a ) corresponding to positive λ ( a ) .22. Start with a point φ n = φ c + (cid:80) Na =1 ζ a ρ ( a ) on the tan-gent space. Evolve φ n by the holomorphic flow bya time T to find φ f , compute the Jacobian J ( φ n ) byintegrating the flow equation for the basis, and thencompute the action S eff ( φ n ).3. Propose new coordinates ζ (cid:48) = ζ + δζ , where δζ is arandom vector chosen with symmetric probability func-tion, that is P ( δζ ) = P ( − δζ ). Evolve φ (cid:48) n = (cid:80) a ζ (cid:48) a ρ ( a ) by the holomorphic flow by a time T to find φ (cid:48) f , com-pute J ( φ (cid:48) n ), and S eff ( φ (cid:48) n ).4. Accept/reject ζ (cid:48) with probability min { , e − S (cid:48) eff + S eff } .5. Repeat from step 3 until a sufficient ensemble of con-figurations is generated.To make the updating effective, we have to account for thefact that the map F T is highly anisotropic. If we considerthe flow close to the critical point, we see that displacementsin the direction of the Takagi vector ρ ( a ) are mapped intovectors that have their magnitude increased by exp (cid:0) λ ( a ) T (cid:1) .Even small differences in the eigenvalues λ ( a ) lead to largedifferences as T increases. If the parametrization space pro-posals δζ are isotropic then the update process becomes in-efficient. Ideally we would like to generate proposals thatare isotropic on the integrations manifold, but since the mapchanges from point to point, this requires care to ensure thatthe detailed balance is preserved. As it turns out this is pos-sible but we will discuss this point later. An easy fix for thisproblem is to adjust the size of displacement for proposalbased on the flow around the critical point. The proposalis then δζ a = exp (cid:0) − λ ( a ) T (cid:1) δ with δ a random variable cho-sen with uniform probability in the interval [ − ∆ , ∆]. Thestep size ∆ is tuned to get reasonable acceptance rates. Ifthe distortions induced by the map F T vary little from φ c to the points sampled by the process, then this algorithm iseffective.By far the most computationally expensive part of thecontraction algorithm—and most other thimble algorithms—is the computation of the Jacobian (even for most bosonicsystems the cost scales with N and N is proportional to thespacetime volume.) Methods to deal with this problem arediscussed in section III.D.Another Metropolis based method was proposed to sam-ple single thimble configurations (Mukherjee et al. , 2013) andwas tested for a single plaquette U (1) problem. In this pro-posal the Jacobian is not included in the sampling and itis to be included via reweighting in the observable measure-ment. This reweighting will fail for most systems that havemore than a few degrees of freedom since for this systemsthe Jacobian fluctuates over many orders of magnitude. HMC on thimbles
A more sophisticated algorithm based on Hybrid MonteCarlo (Duane et al. , 1987) was proposed and tested for the φ model (Fujii et al. , 2015). In principle, a straightforwardextension of HMC could be applied to the action Re S eff onthe parametrization manifold. The problem with such an ap-proach is that it would require the calculation of the deriva-tives of det J , or some related quantity, which is quite cum-bersome. Of course this could be side-stepped by neglectingthe Jacobian in the sampling (Ulybyshev et al. , 2020a), butthis requires reweighting it in the observables which fails forlarge systems. The proposal is then to use HMC as definedby the Hamiltonian in the larger C N space, where the motionis confined to be on the thimble via forces of constraint (Fu-jii et al. , 2015). This has the advantage that the Jacobianis accounted for implicitly, but the algorithm requires solv-ing implicit equations to project back to the thimble. Forthe cases where the thimble is relatively flat/smooth, theseequations can be solved robustly via iteration, as is the casewith the φ system in the parameter range investigated.The basic idea is to integrate the equations of motion gen-erated by the Hamiltonian H ( π, φ f ) = 12 π † π + Re S ( φ f ) , (3.7)subject to the constraint that φ f ∈ T . Forces of constraintperpendicular to the thimble keep the system confined onits surface. The momentum π is in the tangent space at φ f , so it is a real linear combination of columns of J ( φ n ).The perpendicular force has to be a real linear combinationof the columns of iJ ( φ n ), since this forms a basis in thespace perpendicular (according to the scalar product (cid:104) v | w (cid:105) ≡ Re v † w ) to the thimble.For a practical implementation we need to provide an in-tegrator for these equations of motion for finite time steps.A symplectic integrator for this problem is provided by thefollowing method π / = π − ∂ φ f Re S ( φ f ) ∆ t iJ ( φ n ) λ ,φ (cid:48) f = φ f + π / ∆ t ,π (cid:48) = π / − ∂ φ f Re S ( φ (cid:48) f ) ∆ t iJ ( φ (cid:48) n ) λ (cid:48) . (3.8)The map ( π, φ f ) → ( π (cid:48) , φ (cid:48) f ) is symplectic and time reversible,thus satisfying the requirements for HMC. Note that thismap requires the determination of λ and λ (cid:48) , two sets of N real numbers which encode the effect of the constraint forcesacting perpendicular on the thimble. λ is determined by therequirement that φ (cid:48) f ∈ T and λ (cid:48) by requiring that π (cid:48) is in3the tangent space at φ (cid:48) f . For small enough ∆ t , these require-ments lead to unique “small” solutions (which vanish in the∆ t → λ s. A solution for λ (cid:48) can be computed ina straightforward way, via the projection method we discussbelow. Computing λ is more difficult and the current pro-posal is to use an iterative method (Fujii et al. , 2015). Thisiteration is guaranteed to converge for small enough ∆ t , butfor a fixed size ∆ t no guarantees can be made even for theexistence of a solution.With these ingredients in hand, the basic steps of HMCare the following:1. At the beginning of each “trajectory” an isotropic gaus-sian momentum π is generated in the tangent space at φ f , P ( π ) ∝ exp (cid:0) − π † π/ (cid:1) .2. The equations of motion are integrated by repeatedlyiterating the integrators steps above for a t/ ∆ t times,where t is the “trajectory” length.3. At the end of trajectory the proposed ( π (cid:48) , φ (cid:48) f ) are ac-cepted with a probability determined by the change inHamiltonian P acc = min { , exp( −H + H (cid:48) ) } .One important ingredient for this and other algorithms wewill discuss later, is the projection to the tangent space at φ f .If we have the Jacobian matrix in hand J ( φ n ), its columnsform a real basis of the tangent space and the columns of iJ ( φ n ) form a basis for the orthogonal space. Every vector v ∈ C N can then be decomposed in its parallel, P (cid:107) ( φ f ) v , andperpendicular component, P ⊥ ( φ f ) v , using standard algebra.This step is required to find λ (cid:48) in the symplectic integrator.It can also be used to find the starting momentum, at thebeginning of the trajectory: we generate a random vector in C N with probability P (˜ π ) ∝ exp (cid:0) − ˜ π † ˜ π/ (cid:1) and then projectit to the tangent plane π = P (cid:107) ( φ f )˜ π .The projection discussed above can be readily imple-mented when we have the Jacobian matrix J ( φ n ). How-ever, calculating this matrix is an expensive operation thatis likely to become a bottle-neck for simulations of systemswith large number of degrees of freedom. One solution forthis problem is the following (Alexandru et al. , 2017a): weuse the map v → J ( φ n ) v , that maps the tangent space at φ n on the parametrization manifold to the tangent space at φ f on the thimble. This calculation can be implementedefficiently, by solving the vector flow equation, Eq. 3.6, fora single vector v . We extend this to arbitrary vectors thatare not included the tangent space. For a generic vector v we split it into v = P (0) (cid:107) v and iv = P (0) ⊥ v . Here P (0) (cid:107) isthe projection on the tangent space of the parametrizationmanifold, the space spanned by the Takagi vectors, and P (0) ⊥ its orthogonal complement. Both v and v belong to the tangent space at φ n , so J ( φ n ) v , can be computed usingthe vector flow equations. This defines then a map fromany vector v to J ( φ n ) v = J ( φ n ) v + iJ ( φ n ) v , which re-quires two integrations of the vector flow. Using this mapwe can then compute J − ( φ n ) v using an iterative method,such as BiCGstab. It is then straightfoward to prove that P (cid:107) ( φ f ) v = J ( φ n ) P (0) (cid:107) J ( φ n ) − v . Langevin on thimbles
The Langevin algorithm was proposed as possible sam-pling method for single thimble manifolds (Cristoforetti et al. , 2013, 2012a,b). The idea is to sample the thimble man-ifold T with probability density proportional to exp( − Re S )with respect to the Riemann measure induced by embedding T in C N . The residual phase of the measure is taken intoaccount via reweighting. The imaginary part of the action isconstant over the thimble and will not contribute to averages.The Langevin process simulates the evolution of the sys-tem via a drift term due to the action and a brownian motionterm. The discretized version of the process is given by thefollowing updates: φ (cid:48) f = φ f − ∂ φ f Re S ( φ f )∆ t + η √ t (3.9)where the vector η is a random N dimensional vector, in thetangent space at the thimble at φ f .Two details are important here: how the vector η is chosenand how the new configuration φ (cid:48) f is projected back to thethimble. The proposal is to chose η isotropically at φ f bygenerating a gaussian ˜ η unconstrained in C N and then pro-jecting it to the tangent space at φ f using a procedure similarto the projection outlined in the section above, η = P (cid:107) ( φ f )˜ η .This ensures an isotropic proposal in the tangent space andthe norm of the vector is adjusted such that it follows the χ -distribution with N degrees of freedom (Cristoforetti et al. ,2013).At every step we start with φ f on the thimble and wemove along the tangent direction, since both the drift andthe random vector lie in the tangent plane. Unless the thim-ble is a hyperplane, this shift will take us out of the thimble.A projection back to the thimble is required. The meth-ods proposed rely on evolving the new configuration in thedownward flow toward the critical point, projecting there tothe thimble and flowing back (Cristoforetti et al. , 2012a,b).This proposal was found to be unstable (Cristoforetti et al. ,2012b). The only simulations that we are aware of that em-ploy this algorithm involve simulations on the tangent planeto the thimble (Cristoforetti et al. , 2013). In this case theupdates do not require any projection since the manifold is4 Re ϕ Im ϕμ < μ c h ∈ ϕ - ϕ + ϕ Re ϕ Im ϕμ > μ c h ∈ ϕ - ϕ ϕ + Re ϕ Im ϕμ > μ c h ∉ ϕ - ϕ ϕ + FIG. 5: Projections of the thimbles (blue) T , T + , and T − and dual thimbles (red) K , K + , and K − onto the1-complex-dimensional subspace constant fields. The intersection of the original domain of integration with this subspacecorresponds to the real line. The typical arrangement of thimbles varies with the chemical potential.flat. To make this algorithm practical for the general case arobust projection method is needed.A final note about Langevin algorithm: for a finite ∆ t themethod is not exact. Simulations have to be carried out fordecreasing ∆ t and then extrapolated to ∆ t = 0 to removethe finite step-size errors. For other Langevin methods, anaccept/reject step can be used to remove the finite step-sizeerrors, but this has not been developed for thimble simula-tions.While both Langevin method and HMC algorithm performupdates directly on φ f with drift (or force) term evaluatedlocally, it is worth emphasizing that the updates still requirethe integration of the flow equations. This is because theprojection of the shift to the tangent plane to the thimbleand the required projection back to the manifold after theupdate, can only be currently done by connecting φ j with itsimage under the flow φ n in the infinitesimal neighborhoodof the critical point. The advantage of these methods overMetropolis, assuming that a practical projection method isavailable, is that the updates can lead to large change inaction leading to small autocorrelation times in the Markovchain. Case study: bosonic gases
We presently consider the relativistic Bose gas at finitedensity for an application of these algorithms to bosonic sys-tems with sign problems. The continuum Euclidean action of this system is S = (cid:90) d x (cid:2) ∂ φ ∗ ∂ φ + ∇ φ ∗ · ∇ φ + ( m − µ ) | φ | + µ ( φ ∗ ∂ φ − φ∂ φ ∗ ) (cid:124) (cid:123)(cid:122) (cid:125) j ( x ) + λ | φ | (cid:3) , (3.10)where φ = ( φ + iφ ) / √ U (1) symmetry which breaks spontaneously at high density.In Euclidean space, the current j is complex and causes asign problem .This system was studied with the contraction algorithmin (Alexandru et al. , 2016a), HMC method (Fujii et al. ,2013), and the Langevin process (Cristoforetti et al. , 2013).The following lattice discretization of Eq. 3.10 was used S = (cid:88) x,a (cid:34)(cid:16) m (cid:17) φ x,a φ x,a − (cid:88) ν =1 φ x,a φ x +ˆ ν,a − cosh µ φ x,a φ x +ˆ0 ,a + i sinh µ (cid:15) ab φ x,a φ x +ˆ0 ,b + λ (cid:0) φ x,a φ x,a ) − h ( φ x, + φ x, ) (cid:35) , (3.11) This is most readily seen in Fourier space in the continuum: (cid:82) d x j ( x ) = (2 π ) − (cid:82) d p ( − ip ) | φ ( p ) | is purely imaginary. (cid:15) ab is the antisymmetric tensor and (cid:15) = 1. Thislattice action will be used for the remainder of this discus-sion. The final term must be included in the lattice theoryto obtain a well-defined thimble decomposition and we take h small.To apply the contraction algorithm, it is first necessaryto find critical points (extrema) of the action Eq. 3.11. Re-stricting attention to those critical points which are constantin spacetime, the following extremum condition is obtained:(2 + m ) φ − µφ + 2 λ | φ | φ = h . (3.12)Three extrema exist and we denote them φ , φ + , φ − . Thecorresponding Lefschetz thimbles will be denoted T , T + , T − .Depending on the parameters of the theory, different combi-nations of thimbles contribute to the path integral. To thisend, the one-dimensional projections of T , T + , T − depictedin Fig. 5 are useful.For µ < µ c = cosh − (1+ m / T contributes to thepath integral. This is because S R ( φ ± ) < S R ( φ ) for any φ onthe original integration manifold, and therefore no point canflow to φ ± by the upward flow. This is sufficient to eliminate T ± as contributing thimbles.For µ > µ c , the contributing thimbles changes. As seen inthe center of Fig. 5, when h ∈ R , there are flow trajectoriesconnecting both φ − and φ + to φ . This feature, called Stokesphenomenon , introduces complications into the decomposi-tion of the path integral into an integer linear combination ofthimbles. We avoid Stokes phenomenon altogether by sim-ply introducing a complex h ; for a detailed discussion of ourprocedures see (Alexandru et al. , 2016a).Since our purpose is to illustrate the Contraction Algo-rithm, let us consider only the µ > µ c case. As an example,let m = λ = 1 . , h = 0 . i/
10) and µ = 1 .
3. Withthese choices, T + contributes most to the path integral. Theresults obtained on flowed manifolds are plotted in Fig. 6.The variance of S I decreases as a function of flow time; thisdemonstrates that the integral over T + indeed has reducedphase fluctuations relative to R N . Furthermore, the con-vergence of observables as a function of flow time stronglysuggests convergence to T + . B. Generalized thimble method
The main limitation of the methods discussed so far isthat they are capable of computing the integral over onlyone thimble. However, the integral over the real variablesis generically equivalent to the integral over a collection ofthimbles. Finding these collection of thimbles is a daunt-ing process; integrating over all of them an even hardertask. Fortunately, there is a way of bypassing this difficulty T 〈 R e s i du a l 〉 - - - - - - - S I FIG. 6: The imaginary part of the action and the residualphase computed on F T ( T + ) using the contractionalgorithm. The horizontal line denotes the value ofIm S ( φ + ).based on what we learn in section II: the generalized thimblemethod.Recall that if every point of R N (the integration regionof the path integral) is taken to be the initial condition forthe Eq. 2.7 that is then integrated for a time T , we obtaina manifold M T = F T ( R N ) that is equivalent to the initial R N manifold (in the sense that the path integral over R N and M T are the same). In addition, for large enough valuesof T , M T approaches exactly the combination of thimblesequivalent to R N . It is important to understand how thethimbles are approached. In the large T limit an isolatedset of points in R N , let us call each of them ζ c , approachthe critical points φ c of the relevant thimbles. Points nearthem initially approach the critical points but, when close tothem, move along the unstable directions, almost parallel tothe thimble but slowly approaching it (see Fig. 4). Pointsfar from ζ c run towards a point when the action diverges,either at infinity or at a finite distance (in fermionic theoriesthimbles meet at points where the action diverges as exem-plified by the Thirring model discussed below). This meansthat the correct combination of thimbles equivalent to theoriginal path integral can be parametrized by points in R N .This is an advantage over the contraction method where only6 ⊙ ⊙ ⊙ ⊙ - - - A A T = T = T = FIG. 7: Complex A = 1 /V (cid:80) x A ( x ) plane for theThirring model. The blue squares are critical points, theblue lines the thimbles. The dashed line is the tangentspace to the “main” thimble while the other solid lines arethe manifolds M T obtained by flowing the tangent spaceby T = 0 . , .
05 and 0 .
5. Notice how M T approaches thecorrect combination of thimbles as T is increased.one thimble at a time could be parametrized. We have then (cid:90) R N dφ e − S ( φ ) = (cid:90) M T d ˜ φ e − S ( ˜ φ ) det J ( ˜ φ )= (cid:90) R N dζ e − S [ F T ( ζ )] det J ( ζ ) . (3.13)The generalized thimble method consists in using aMetropolis algorithm on R N with the action Re S eff where S eff ( ζ ) = S [ F T ( ζ )] − (log det J ( ζ )). Generalized Thimble Algorithm (GTA)
1. Start with a point ζ in R N . Evolve it by the holomor-phic flow by a time T to find φ f = F T ( ζ ).2. Propose new coordinates ζ (cid:48) = ζ + δζ , where δζ is arandom vector drawn from a symmetric distribution.Evolve it by the holomorphic flow by a time T to find φ (cid:48) f = F T ( ζ (cid:48) ).3. Accept ζ (cid:48) with probability P acc = min { , e − ∆ Re S eff } .4. Repeat from step 2 until a sufficient ensemble of con-figurations is generated.Methods to speed up—or bypass—the frequent computationof the Jacobian J are an improvement of the method andwill be discussed below (see III.D).While the algorithm above is exact, the practical appli-cability of the GTA depends on the landscape induced byexp( − Re S eff ) on M T . At large T , the points ζ that are - - - μ < δ ψψ > FIG. 8: The difference in the value of the chiral condensatebetween the exact result and the one obtained by thecontraction method (with T = 2) shown in red and thegeneralized thimble method (with T = 0, that is,integration over the tangent space. ) The parameters are N = 8, m = 1 and g = 1 / M T lie onsmall, isolated regions. This explains why the phase of theintegrand fluctuates less on M T than on R N . The imag-inary part of S [ F T ( ζ )] on points on M T are the same asthe imaginary parts of the action S ( ζ ) in a little regionaround ζ c , the only region with significant statistical weightexp( − Re S eff [ F T ( ζ )]).In between the regions around the different ζ c lie areaswith small statistical weight exp( − Re S eff [ F T ( ζ )]) that aremapped to points where the action (nearly) diverges, as wediscussed in II.C. A probability landscape of this form maytrap the Monte Carlo chain in one of the high probabilityregions, breaking ergodicity. A trapped Monte Carlo chainis effectively sampling only one of the thimbles contributingto the integral (more precisely, it is an approximation to aone thimble computation). This problem can be alleviatedby making T small. In that case M T will be farther awayfrom the thimbles, the phase oscillations are larger and theoriginal sign problem may not be controlled. The usefulnessof the GTA relies then in being able to find a value of T such that the sign problem is sufficiently ameliorated whilethe trapping of the Monte Carlo chain is not a problem. Inseveral examples discussed below, over a large swatch of pa-rameter space, it is not difficult to find a range of values of T for which the GTA is useful. Still, one should perform duediligence and try to diagnose trapping signs in every calcu-lation, as it is always the case in Monte Carlo calculations. Case study: 0+1D Thirring model
We will use the finite density/temperature Thirring modelin 0+1 , A S I A S I FIG. 9: Histogram of average field A (left) and imaginarypart of the action (right) in a Monte Carlo sampling in the1D Thirring model using the GTM with T = 0 (top line)and T = 0 . g a = 1 / , N = 32 calculation.In the T = 0 calculation the phase e − iS I fluctuates toowildly and the result has large uncertainties. On the T = 0 . T = 0 . L = ¯ ψ a ( i∂/ + m + µγ ) ψ a + g N F ¯ ψ a γ µ ψ a ¯ ψ a γ µ ψ a , (3.14)where φ is a spinor for the appropriate spacetime dimensionand a indexes the N F different flavors of fermions. Thistheory is, in 1 + 1 dimensions, asymptotically free. The N F case is identical to the Gross-Neveu model and its groundstate breaks a discrete symmetry spontaneously and, in thisrespect, resembles QCD. For N F > d dimensions is: S = (cid:88) x,ν N F g (1 − cos A ν ( x )) + (cid:88) x,y ¯ ψ a ( x ) D xy ψ a ( y ) , (3.15) with D Wxy = δ xy − κ (cid:88) ν =0 , (cid:2) (1 − γ ν ) e iA ν ( x )+ µδ ν δ x + ν,y + (1 + γ ν ) e − iA ν ( x ) − µδ ν δ x,y + ν (cid:3) (3.16)with 1 /κ = 2 m + 4 d or D KSxy = mδ xy + 12 (cid:88) ν =0 , (cid:2) η ν ( x ) e iA ν ( x )+ µδ ν δ x + ν,y − η † ν ( x ) e − iA ν ( x ) − µδ ν δ x,y + ν (cid:3) (3.17)with η ( x ) = 1 , η = ( − x , η = ( − x + x and the flavorindex goes from 1 to N F in the Wilson fermion case butfrom 1 to N F / A ν ( x ) leads to a discretized version of Eq. 3.15,showing their equivalence. Integration over the fermion fieldsleads to purely bosonic action more amenable to numericalcalculations: S = N F (cid:32) g (cid:88) x,ν (1 − cos A ν ( x )) − γ log det D ( A ) (cid:33) , (3.18)with γ = 1 (Wilson) or γ = 1 / N F Dirac fermions in the continuum. Thepresence of the chemical potential µ renders the fermion de-terminant complex and is the origin of the sign problem inthis model.The 0 + 1 dimensional case can be solved exactly with thelattice action in Eq. 3.18 and it has been used as a checkon several methods designed to handle sign problems (Fujii et al. , 2017; Li, 2016; Pawlowski and Zielinski, 2013). Itsthimble structure is known. In the A ( x ) = constant sectorit is shown in Fig. 7. There is one purely imaginary criticalpoint that has the smallest value of the real part of the ac-tion, therefore called the “main critical point”. Therefore, inthe semiclassical limit it should dominate the path integral.Thimbles touch each other at points where the fermion de-terminant vanishes and the effective bosonic action diverges(shown as blue squares in Fig. 7). The tangent space to themain thimble ( T ) is just the real space shifted in the imag-inary direction (dashed red line in Fig. 7). The integrationover the tangent space is no more expensive than over thereal space since no flowing is required and the Jacobian of thetransformation is one. The tangent space, lying parallel tothe real space, has the same asymptotic behavior as R N andis equivalent to it for the computation of the integral. Thefigure also shows the result of “flowing” the tangent space bydifferent values of T ; the larger the value of T , the closer theresulting manifold(s) approach the thimbles. Starting from8the tangent space and using a flow time T = 2 the manifold F T ( T ) obtained is nearly indistinguishable from the thim-bles.In (Alexandru et al. , 2015) the model was studied using thecontraction algorithm. The results, shown on Fig. 8 indicatethat the fermion condensate, for instance, is close to the ex-act result but does not agree with it, in particular for certainvalues of µ near the transition from (cid:104) ¯ ψψ (cid:105) = 0 to (cid:104) ¯ ψψ (cid:105) (cid:54) = 0.The size of the discrepancy is consistent with a semiclassicalestimate of the contributions of other thimbles (besides themain thimble). Similar behavior was seen a 1-site model offermions (Tanizaki et al. , 2016). The integration over thetangent space, however, gives the correct result. Of course,the average sign on the tangent space is smaller than the oneobtained with the contraction method. For not too low tem-peratures the sign fluctuation is, however, small enough toallow for the computation to be done on the tangent plane.But as the temperature is lowered, the sign fluctuations growand it becomes difficult to sample the correct distribution,as predicted by general arguments (see Eq. 1.4). One canthen use the generalized thimble method and integrate onthe manifold F T ( T ) for a suitable value of T . Too small a T the sign fluctuation is too large; a T too large is essentiallyan integration over one thimble and the wrong results is ob-tained. It is interesting to understand how the transitionbetween these two behaviors occur. In Fig. 9 histogramsof the imaginary part of the effective action are shown forboth T = 0 and T = 0 .
5. It is clear that for T = 0 . T = 0 (no flow) the dis-tribution is broader. Consequently, the values of the phaseexp( − i Im S ) fluctuate less when there is flow and the signproblem is minimized. On the other hand, for large enoughflow time, the probability distribution exp( − Re S ) becomesmultimodal and the trapping of Monte Carlo chains can pre-vent proper sampling. Thus, the GTM trades the sign prob-lem by a the problem of sampling a multimodal distribution.This trade is not without profit: in many cases one can findvalues of T such that the sign problem is sufficiently allevi-ated but trapping has not set in yet. These values of T canbe determined by trial and error. As T is increased trappingoccurs, quite suddenly, and it is not difficult to detect it bynoticing a jump on the values of the observables. Also, thereare well studied ways to deal with trapping, as explainedin the next section. Still trapping is a source of concern inGTM calculations and other, more general techniques, havebeen developed to avoid it (see section IV). μ / m f R e 〈 e - i S I 〉 N T = T = 〈 n 〉 / m f FIG. 10: Fermion density (top) and average sign (bottom)of the 1 + 1dimensional Thirring model on a 10 ×
10 lattice, g = 1 , m = − .
25 (lattice units). The sign problem isstrongly suppressed and one moves the path integrationfrom R N to the tangent plane T and from that to theflowed manifold F T ( T ) allowing for precise measurements ofthe density and other observables (Alexandru et al. , 2017b). Case study: 1+1D Thirring model
The lessons learned in applying the generalized thimblemethod to the 0 + 1 dimensional Thirring model carry on tothe more interesting 1 + 1 dimensional case. Extensive calcu-lations on the finite density/temperature 1 + 1 dimensionalThirring model with two flavors were made over a range ofparameters in the strong coupling region (Alexandru et al. ,2017b) with both Wilson and staggerred fermions. The thim-ble structure of the 1 + 1 models is more complex than the0 + 1 case. Still, all critical points/thimbles present in the0 + 1 dimensional case have analogues in 1 + 1 dimensions(which has many others without a 1 + 1 dimensional ana-logue). It is still true that the closest critical point to the realspace (the “main critical point”) is a constant shift of A ( x )by an imaginary amount and that its tangent space is justa translation of R N by an imaginary amount (see Fig. 7).The path integration over R N has a bad sign problem forall values of the chemical potential larger than the fermion(renormalized) mass ( µ > m f ), that is, for all values of µ forwhich there is an appreciable number of fermion-antifermion9unbalance . The integration over the tangent space of themain thimble can be accomplished at no extra cost by sim-ply shifting the variables of integration by a constant imag-inary amount. This step, by itself, improves the sign prob-lem considerably. The reason is that the tangent space is a(rough) approximation to the main thimble, specially the re-gion near the critical point that dominates the path integralin the semiclassical regime. Still, for larger volumes, smallertemperatures and higher chemical potential, the shift to thetangent space is not enough to control the sign fluctuation.It was determined that flow times of the order of T = 0 . et al. , 2017b) it wasalso demonstrated that the same method works well as thecontinuum and thermodynamic limits are approached. C. Trapping and tempered algorithms
The landscape induced by exp( − Re S eff ) on theparametrization manifold changes as a function of the flowtime T . For small T the landscape is typically flat, while forlarger T the landscape is steeper. When the sign problem issevere enough to require large flow times, the landscape ofexp( − Re S eff ) has high peaks and low valleys and the prob-ability distribution can become multi-modal. The purposeof this section is to detail several algorithms addressing thisdifficulty.We first discuss the method of tempered transitions (Neal,1996). Designed to combat trapping, a tempered pro-posal is a composite proposal assembled from small stepswhich, taken together, more rapidly cover phase space thana standard proposal. A tempered proposal is constructedas follows. First, let p ( φ ) , p ( φ ) , ..., p n ( φ ) be a sequenceof increasingly relaxed probability distributions such that p ( φ ) ≡ p ( φ ) is the distribution of interest and p n ( φ ) is sig-nificantly more uniform. Next, for every i , let ˆ T i be a tran-sition probability satisfying detailed balance with respect to p i , that is p i ( φ ) ˆ T i ( φ → φ (cid:48) ) = p i ( φ (cid:48) ) ˆ T i ( φ (cid:48) → φ ) . (3.19)Then a tempered update ˆ T is executed by first generating a We note here that, contrary to other approaches, the thimble methodtrivially reproduces the “Silver Blaze” phenomena, the fact that thesystem is trivial at small temperatures and chemical potentials smalerthan the mass of the lightest fermionic excitation (Cohen, 2003). sequence of 2 n configurations φ → φ → ... → φ n ≡ φ (cid:48) n → φ (cid:48) n − → ... → φ (cid:48) , (3.20)using transition probabilities ˆ T , ˆ T , . . . , ˆ T n , ˆ T n , . . . , ˆ T , fol-lowed by an accept/reject step with probability: P acc ( φ → ... → φ (cid:48) ) = min { , F ( φ ) /F ( φ (cid:48) ) } . (3.21)where F ( φ ) ≡ p ( φ ) p ( φ ) p ( φ ) p ( φ ) · · · p n − ( φ n − ) p n − ( φ n − ) p n ( φ n − ) p n − ( φ n − ) . (3.22)What is gained by using tempered proposals is enhancedergodicity. Since the distributions p i are increasingly uni-form, the corresponding transition probabilities ˆ T i may growin support without decreasing the acceptance probability.To apply this general framework to simulations trapped byholomorphic gradient flow, suppose the flow time T is largeenough that the probability distribution of interest p ( ζ ) = p ( ζ ) = e − Re S eff ( ζ ) Z (3.23)is multi-modal. Consider a sequence of flow times T < T <... < T n such that T = T and T n (cid:28) T . This defines asequence of probability distributions p ( ζ ) , p ( ζ ) , . . . , p n ( ζ )which are decreasingly multi-modal; we use this sequence toperform tempered proposals.Applying this method to the (0+1) dimensional ThirringModel at finite density (Alexandru et al. , 2017c), severelytrapped simulations have been liberated. Certain thermo-dynamic parameters exist for which at least five thimblescontribute non-negligibly to the path integral. Trapping toa single thimble, however, can become arbitrarily severe: forexample, at T = 0 .
5, the multi-modality of p ( ζ ) is so se-vere that over the course of a Metropolis with 10 steps nota single transition occurred. Tempered proposals free thesetrapped MCs however; this is demonstrated in Fig. 12 whereproper sampling of the T = 0 . parallel tempering , was proposed to helpsample such from multi-modal distributions (Earl and Deem,2005; Geyer, 1991; Swendsen and Wang, 1986). Parallel tem-pering involves simulating n replicas of the system of in-terest, each having a particular value of the tempering pa-rameter. Each stream evolves separately and swaps betweenreplicas are added satisfying detailed balance. The swap-ping of configurations between adjacent replicas leads to en-0FIG. 11: Here we show the evolution of field space sampledin the (0+1) dimensional Thirring model as a function offlow time. At small flow times the distribution is relativelyuniform and much of phase space is sampled. Thedistribution sharpens as the flow time increases, and atsufficiently large flow times, the shoulder thimbles centeredabout ± . - - - A FIG. 12: Plotted is A = N − t (cid:80) t A ,t at after eachtempered transition. The five heavily visited positions infield space correspond to five thimbles contributing to thepath integral. This distribution is to be compared with thesharpest distribution of Fig. 11, where only one thimble issampled.hanced ergodicity relative to the single chain case. Fukumaet. al. have developed the “Tempered Lefschetz ThimbleMethod” (TLTM), an application of parallel tempering tomulti modal distributions generated by flow (Fukuma andUmeda, 2017). As with tempered transitions, in this methodthe flow time is chosen as a tempering parameter. The TLTMmethod has been successfully applied to the (0+1) Thirringmodel (Fukuma and Umeda, 2017) where trapping due to flow times as large as T = 2 . . The au-thors also studied how to pick the flow times optimally anddevised a geometric method for this optimization (Fukuma et al. , 2018). More recently, the TLTM has been applied tothe Hubbard model away from half filling on small lattices(Fukuma et al. , 2019). D. Algorithms for the Jacobian
The most computationally expensive part of many algo-rithms involving deformation of contours in field space – likethe contraction or the generalized thimble method – is thecalculation of the Jacobian J related to the parametrizationof the manifold of integration. For bosonic systems wherethe Hessian can be computed efficiently the calculation timeis dominated by the matrix multiplication in the flow equa-tion and its computation complexity is O ( N ) where N isproportional to the spacetime volume of the theory. Thecalculation of det J has also similar computational complex-ity. This prohibitive cost prevents the study of all but thesmallest models.Fortunately, there are ways of bypassing this large cost.In Ref. (Cristoforetti et al. , 2014) a stochastic estimator wasintroduced to compute the phase , Φ( φ n ) = arg det J ( φ n ).The main idea stems from the observation that the Jaco-bian can be expressed as J = U R for some unitary ma-trix U and some real, upper-triangular matrix R , a prop-erty that follows from the fact that J † J ∈ R ; thereforearg det J = arg det U . Note that since J and U are relatedby a real matrix, this corresponds to a change in basis in thetangent plane, so the columns of U form a basis of the tan-gent space too, an orthonormal basis. Moreover U satisfies d log det U ( t ) /dt = − i Im Tr (cid:0) U T ( t ) H ( t ) U ( t ) (cid:1) . The trace canbe estimated stochastically by using random vectors ξ ∈ R N with (cid:104) ξ i ξ j (cid:105) = δ ij , where the average is taken over the randomsource; if we generate N R vectors we haveTr (cid:0) U T ( t ) H ( t ) U ( t ) (cid:1) ≈ N R N R (cid:88) r =1 ( ξ ( r ) ) T U T ( t ) H ( t ) U ( t ) ξ ( r ) . (3.24)Now η ( r ) = U ξ ( r ) is a random vector in the tangent plane,isotropically distributed and its length, with respect to thereal Euclidean metric, satisfies (cid:104) (cid:10) η ( r ) (cid:12)(cid:12) η ( r ) (cid:11) (cid:105) = N . We cangenerate such vectors without computing U : we generate a Because the thermodynamic parameters used in (Fukuma andUmeda, 2017) do not match those in (Alexandru et al. , 2017c) itis currently not possible to compare the efficacy of tempered transi-tions and the TLTM. A comparison would, however, be useful. η isotropically in C N with (cid:104) ˜ η † ˜ η (cid:105) = 2 N , forexample using a gaussian distribution P (˜ η ) ∝ exp (cid:0) − ˜ η † ˜ η/ (cid:1) ,and then project it to the tangent space η = P (cid:107) ( φ f )˜ η usingthe same procedure presented when we discussed the HMCalgorithm. Using i Φ( t ) = log det U ( t ), the phase can then beestimated fromΦ( T ) ≈ Φ(0) − Im (cid:90) T dt N R N R (cid:88) r =1 Tr (cid:16) η ( r ) ( t ) T H ( t ) η ( r ) ( t ) (cid:17) , (3.25)whose computational cost scales as O ( N × N R ). By compar-ing this stochastic estimation algorithm by explicit computa-tion for a complex φ theory, Ref. (Cristoforetti et al. , 2014)presented numerical evidence that this algorithm indeed pro-vides a nontrivial speedup for the computation of the residualphase in relatively large systems. However, its applicabilityis limited to the phase of the Jacobian; the GTM requiresthe magnitude also.For methods that require the Jacobian, we can substitutethem with computationally cheap estimators . The idea isto use the estimators during the generation of configurationsand correct for the difference when computing the observ-ables. Two estimators for log det J have been introducedin (Alexandru et al. , 2016b). They are given by W = (cid:90) T dt (cid:88) i ρ ( i ) † H ( t ) ρ ( i ) , W = (cid:90) T dt Tr H ( τ )(3.26)where ρ ( i ) are the Takagi vectors of H ij (0) with positiveeigenvalues. The first estimator, W , is equal to log det J forquadratic actions. The second estimator is equal to ln det J when the Jacobian is real along the flow. As such, it is ex-pected to be a good estimator for Jacobians which are mostlyreal. The bias introduced by the use of estimators insteadof the Jacobian is corrected by reweighting the differencebetween them when computing observables with the help of: (cid:104)O(cid:105) = (cid:104)O e − ∆ S (cid:105) Re S (cid:48) eff (cid:104) e − ∆ S (cid:105) Re S (cid:48) eff (3.27)where S (cid:48) eff = S − W , and ∆ S = S eff − Re S (cid:48) eff . The estimatoris useful when ∆ S has small fluctuations over the sampledthe field configurations, that is, if W , “tracks” log det J well.For theories where the Hessian can be computed efficiently,for example for bosonic theories with local actions, W esti-mator has computational cost of O ( N ) and W has O ( N )complexity, a significant improvement over O ( N ) for thefull Jacobian. In order to use Eq. 3.27 the correct Jacobian J needs to be computed. This has to be done, however, onlyon field configurations used in the average in Eq. 1.2. Typ- ically, configurations obtained in subsequent Monte Carlosteps are very correlated and only one configuration out oftens or hundreds of steps are used in Eq. 3.27. The ideais then to use the cheaper Jacobian estimators, like W , W during the collection of configurations and to compute theexpensive Jacobian J only when make measurements, whichcheapens the calculation by orders of magnitude. This strat-egy was used, for instance, in the φ model in 3 + 1 dimen-sions (Alexandru et al. , 2016a) and the Thirring model in1 + 1 dimensions (Alexandru et al. , 2017b), both at finitedensity. However for other class of problems, such as realtime systems, the estimators W , W do not provide a signif-icant improvement.A rather more robust algorithm for the Jacobian havebeen introduced in Ref. (Alexandru et al. , 2017a). The keyidea is to modify the proposal mechanism in such a wayas to incorporate the Jacobian as part of the effective ac-tion. As an added bonus, the procedure leads to isotropicproposals on the integration manifold. As in the contrac-tion algorithm, the goal is to generate a distribution on theparametrization manifold with probability proportional toexp[ − Re S eff ( φ n )]. This is a Metropolis method, so we needto make a proposal and then accept/reject it. For updateproposals, we generate a random complex vector in the tan-gent plane at φ f , uniformly distributed with normal distri-bution P ( η ) ∝ exp (cid:0) − η † η/δ (cid:1) . The parameter δ controls thestep-size and is tuned to optimize the acceptance rate. Thevector η is generated using the projection discussed earlier:a ˜ η ∈ C N sampled from a Gaussian distribution and then η = P (cid:107) ( φ f )˜ η using the vector flow projection. The update inthe parametrization space is φ (cid:48) n = φ n + (cid:15) where (cid:15) = J − ( φ n ) η is a vector in the tangent space at φ n . Here we take advan-tage of the fact that the parametrization space is flat and φ (cid:48) n does not need to be projected.Since the proposals are not symmetric, the accept/rejectstep has to be slightly modified to satisfy detailed bal-ance. The added factor does not cancel the Jacobian,unless the proposal satisfies an implicit equation that isnot easy to solve. A better alternative is based an algo-rithm by Grady (Grady, 1985): the ratio of Jacobians istaken into account implicitly using a stochastic generatedvector. The vector is generated with probability P ( ξ ) ∝ exp (cid:2) − ξ † ( J (cid:48)† J (cid:48) ) ξ (cid:3) , where J (cid:48) = J ( φ (cid:48) n ) and the proposal is ac-cepted with probability (Alexandru et al. , 2017a): P acc = min { , e − Re[ S (cid:48) − S ]+ ξ † ∆ Jξ − (cid:15) † ∆ J(cid:15) } , (3.28)where ∆ J = ( J (cid:48)† J (cid:48) ) − ( J † J ). We stress that ξ is a complex random vector with 2 N independent components, whereas (cid:15) has only N independent components.The highlight of this method is that by construction is2samples the probability distribution e − Re S | det J | withoutan explicit computation of | det J | . It only requires the com-putation of J − η and J(cid:15) both of which scale as O ( N ) formost bosonic theories.A simplified algorithm that may lead to further computa-tional speedup can be achieved when instead of J ( φ n ) we ap-proximate it with J ( φ c ). The Jacobian is then only requiredto compute the displacements (cid:15) and the accept/reject is donesimply based on the change of the action since ∆ J = 0. Forthis method J ( φ c ) − can be computed once at the start ofthe simulation. Of course, the difference between J ( φ n ) and J ( φ c ) has to be included by reweighting the observables asit was done with W , . This method should work well whenthe fluctuations of J ( φ n ) are mild. In Ref. (Alexandru et al. ,2017a) this was to shown to be the case for the real timestudy of a 1+1 dimensional φ theory even in the stronglycoupled regime. Case study: real time field theory
The generalized thimble method and the whole machineryused in dealing with the computational cost of the Jacobianwas applied to one of the most challenging sign problems: thecalculation of real time correlators in field theory. These cor-relators are the building blocks for the computation of trans-port coefficients like diffusivity, conductivity, viscosities, etc.,and are of great importance in a variety of physical contexts.Similar methods can also be used in fully non-equilibriumsituations. At the same time the available theoretical toolsto study this problem are limited. Even perturbation theoryrequires complicated resummations and in the strongly cou-pled regime the conventional lattice methods are not applica-ble as detailed below. Alternatively, stochastic quantization(or “complex Langevin”) have been utilized but it seems toconverge to the wrong result if the time separation t − t (cid:48) ismore than the inverse temperature β (Berges et al. , 2007).The central objects of interest here are time dependent correlation functions of the form (cid:104)O ( t ) O ( t (cid:48) ) (cid:105) β = Tr(ˆ ρ O ( t ) O ( t (cid:48) )) (3.29)where ˆ ρ is the density matrix which reduces to the fa-miliar Boltzmann factor, e − βH /T r ( e − βH ), in equilibrium.Time dependent correlation functions can be generated fromthe Schwinger-Keldysh (SK) path integral (Keldysh, 1964; t max t max - i β / - i β / - i β C Re ( t ) Im ( t ) FIG. 13: The Schwinger-Keldysh contour in complex timeplane. The real part corresponds to forward and backwardstime evolution and the imaginary part corresponds to theinsertion of the equilibrium density matrix.Schwinger, 1961), (cid:104)O ( t ) O ( t (cid:48) ) (cid:105) β = Tr (cid:104) O (0) e − iH ( t − t (cid:48) ) O (0) e iH ( t − t (cid:48) + iβ ) (cid:105) = 1 Z (cid:90) D φ e iS SK [ φ ] O ( t ) O ( t (cid:48) ) , (3.30)(3.31)where the SK action is obtained by integrating the La-grangian over a complex contour, shown in Fig. 13. Thereal part corresponds to forward and backward time evo-lution and the imaginary part corresponds to the insertionof the equilibrium density matrix, e − β ˆ H / Tr e − β ˆ H . For in-stance, a discretized Schwinger-Keldysh action for a scalartheory reads: , S ( φ ) = (cid:88) t,(cid:126)x a t a
12 ( φ t +1 ,(cid:126)x − φ t,(cid:126)x ) a t + 12 (cid:88) ˆ i ( φ t,(cid:126)x +ˆ i − φ t,(cid:126)x ) a + V ( φ t,(cid:126)x ) (cid:35) ,a t = ia for 0 ≤ t < N t − ia for N t ≤ t < N t a for 2 N t ≤ t < N t + N β (3.32)from which the correlators follow (cid:104) φ t ,(cid:126)x φ t ,(cid:126)x (cid:105) = (cid:82) (cid:0) (cid:81) t,(cid:126)x dφ t,(cid:126)x (cid:1) e − S ( φ ) φ t ,(cid:126)x φ t ,(cid:126)x (cid:82) (cid:0) (cid:81) t,n dφ t,(cid:126)x (cid:1) e − S ( φ ) . (3.33) For simplicity we consider the bosonic case but the formalism can begeneralized to the fermionic case in a straightforward fashion. We alsoinclude an overall factor i so that the associated Boltzmann weightis e − S . ≤ t < N t , is pure imaginary as expected from thereal time evolution and leads to a severe sign problem. Infact, due to the fact that its pure phase with no dampingterm, it is impossible to define a “phase quenched” measureand reweigh the phase. For this reason conventional latticemethods do not work not for real time problems even if un-limited computational power is available . By contrast, onany manifold M that is obtained by flowing from R N bysome fixed flow time, Re S > et al. , 2016c;Mou et al. , 2019a) and 1+1 dimensional bosonic field theorieswith V ( φ ) = λφ /
4! potential. In Figs. 14 and 15 the twolowest spatial Fourier modes of the time-ordered correlator C ( t − t (cid:48) , p ) = T (cid:104) φ ( t, p ) φ ( t (cid:48) , p ) † (cid:105) β (3.34)where φ ( t, p ) = 1 N x N x − (cid:88) x =0 e ipx φ tx (3.35)are plotted for different values of λ (Alexandru et al. , 2017a).To ensure the validity of the method the weak coupling( λ = 0 .
1) Monte-Carlo result is compared with the zeroth,first and second order perturbation theory calculations per-formed analytically. In the strong coupling regime which liesoutside of the domain of perturbation theory (see Fig. 16)the method works as well as it does in the weak couplingregime without any problems. In the quantum mechanicalcase a similar cross check has been performed which showedagreement between the Monte-Carlo results and the exactresult obtained from numerically solving the Schr¨odingerequation (Alexandru et al. , 2016c). In Refs. (Mou et al. ,2019a,b) the 1+1 dimensional model was studied with a non-equilibrium density matrix.The sign problem in the real time problem gets more severewhen the time interval between the operators, | t − t (cid:48) | , in unitsof inverse temperature is increased. This is because the realpart of C that generates the pure phase contribution to thepath integral becomes larger. Therefore a larger flow time isneeded to handle larger | t − t (cid:48) | . Currently with the help of In principle, it is possible to extract the real time correlator (3.30)from a purely Euclidean time correlator by analytic continuation.The extrapolation is, however, numerically unstable and requires ex-ponentially accurate precision in Euclidean time. + + + + + + + + +× × × × × × × × ×* * * * * * * * * - - - - t' I m C ( t - t ', ) + + + + + + + + +× × × × × × × × ×* * * * * * * * * λ = λ = λ = + O ( λ ) × O ( λ ) * O ( λ )- R e C ( t - t ', ) FIG. 14: The Monte Carlo computation of time ordercorrelation function defined in Eq. (3.34) for p = 0 and λ = 0 . , . ,
1. The λ = 0 . O ( λ ), O ( λ )and O ( λ ) which are offset in the x axis for visual clarity. λ = λ = λ = - - - t' I m C ( t - t ', π / N x ) - - R e C ( t - t ', π / N x ) FIG. 15: The Monte Carlo computation of time ordercorrelation function defined in Eq. (3.34) for p = 2 π/N x and λ = 0 . , . , | t − t (cid:48) | = 4 β has been achieved on a lattice with N t = 8, N β = 2, and N x = 8. Extending to larger time separations seem to behindered by trapping in a local minima in the Monte-Carlo(Metropolis-Hastings) evolution which calls for alternativesampling methods to be utilized.4 λ R e C ( , π / N x ) R e C ( , ) FIG. 16: The comparison between the perturbativecalculation and the Monte-Carlo result.The dotted, dashedand solid lines denote O ( λ ), O ( λ ) and O ( λ ) calculationsrespectively E. Gauge theories
The thimble structure of gauge theories is more compli-cated due to the fact that critical points are not isolatedpoints but they are continuous manifolds formed as gaugeorbits. The thimbles attached to the critical points carrythe same degeneracy due to the gauge symmetry. One mightenvision instead of working with the degenerate field space M fixing the gauge and working with the quotient space M / G where the critical points would be isolated and Picard-Lefschetz theory can be used as usual. As we discuss below,this is possible for an abelian gauge fields, however not fornon-abelian gauge fields. The reason for that is that somecritical points have nontrivial stabilizers and become singularpoints on M / G . In this case Picard-Lefschetz theory hasto be modified to accommodate these complications whichhas been discussed in the context of Chern Simons theory in(Witten, 2011)Lattice gauge theory remains largely unexplored from theperspective of Picard-Lefschetz theory at the moment . We For example the zero field configuration is such a critical point sta-bilized by the whole gauge group. Though exploratory studies exist (Di Renzo and Eruzzi, 2018;Pawlowski et al. , 2020). review a few exploratory examples from the literature in thenext two sections. Before doing so we first discuss somegeneralities. In lattice gauge theory the fundamental degreesof freedom are gauge links, U i where i ≡ ( x, µ ) collectiveindex for the link variable U xµ ≡ U ( x + ˆ µ, x ) ≡ U i . Thederivative with respect to the link variable is defined as D ai f ( U ) ≡ ∂∂t f (cid:16) e itT a U i (cid:17)(cid:12)(cid:12)(cid:12) t =0 . (3.36)As usual we consider the complexification of the Lie groupwhere the link variables can be parameterized as U = e iξ a T a where T a are the group generators and ξ a are complex vari-ables. For example the complexification of SU ( N ) leads to SL ( N ). The holomorphic flow equation reads dU i dτ = i (cid:88) a ( T a D ai S ( U )) U i (3.37)and it satisfies the desired properties d Re S ( U ) /dt = |D ai S ( U ) | > d Im S ( U ) /dt = 0. Unlike ordinaryderivatives, the groups derivatives do not commute,[ D ai , D bj ] = − f abc δ ij D cj , [ ¯ D ai , ¯ D bj ] = − f abc δ ij ¯ D cj , [ D ai , ¯ D bj ] = 0 , (3.38)where f abc are the structure constants of the gauge groupsuch that [ T a , T b ] = if abc T c . Therefore the flow equationfor the tangent space generated by e ai = D ai S is modified as de ai dτ = e bj D bj D ai S − f abc e bi D ai S . (3.39)
Case study: Heavy-Dense QCD
Some exploratory work towards implementing the thimblemethod in QCD has been done within the so-called “heavy-dense QCD” that is QCD with heavy quarks in the highdensity limit (Zambello and Di Renzo, 2018). As opposedto the heavy mass limit where the quarks decouple from thetheory, in the simultaneous high-mass, high-density limit m → ∞ , µ → ∞ , e µ /m : fixed , (3.40)the quarks remain in the picture and the theory has a non-trivial phase structure controlled by µ . Just like QCD,heavy-dense QCD also exhibits a sign problem. At the sametime it is not as computationally demanding as full QCDwhich makes it a fruitful arena for testing new approaches tothe sign problem (Aarts et al. , 2016; Zambello and Di Renzo,2018). In this limit the fermion determinant simplifies quite5dramatically as (Bender et al. , 1992; Blum et al. , 1996)det D f → (cid:89) (cid:126)x det (1 + γP (cid:126)x ) det (cid:0) γP − (cid:126)x (cid:1) (3.41)where γ ≡ (2 e µ /m ) N t and ˜ γ ≡ (2 e − µ /m ) N t and P (cid:126)x = (cid:81) N t − t =0 U ( (cid:126)x, t ) is the Polyakov loop. Eq. (3.41) has a simplephysical interpretation: in the infinite mass limit, quarks arepinned to their spacial location and do not move. Therefore aquark (anti-quark) at a spatial point (cid:126)x is simply described bythe Polyakov loop P (cid:126)x ( P − (cid:126)x ). Furthermore due to the highdensity limit, the anti-quark contribution is negligible (i.e. γ (cid:29) ˜ γ ) and one can neglect the second determinant in theright hand side of Eq. (3.41). Since the fermion determinanthas no dependence on the spatial links, U µ (cid:54) =0 , one can obtainthe effective action for heavy-dense QCD by integrating outthe spatial degrees of freedom in the QCD path integral, Z QCD = (cid:90) DU µ e β Nc (cid:80) p ( Tr U p +Tr U † p ) N f (cid:89) f =1 det D f → (cid:90) DU e − S HD ( U ) (3.42)where we assumed all N f quarks are heavy and has identicalchemical potentials for simplicity. The effective action forthe heavy-dense QCD in this case is S HD ≡ S gauge − N f (cid:88) (cid:126)x log [det (1 + γP (cid:126)x )] ,S gauge ≈ − (cid:18) β (cid:19) N t (cid:88) (cid:104) (cid:126)x(cid:126)y (cid:105) (cid:16) Tr P (cid:126)x Tr P − (cid:126)y + Tr P (cid:126)y Tr P − (cid:126)x (cid:17) (3.43)The leading order pure gauge action, including the coefficient( β/ N t , follows from the character expansion of the orig-inal gauge action (Langelage et al. , 2014). In the low tem-perature limit where N t (cid:29) β ∼ O (1) the pure gaugecontribution, S gauge , can further be neglected and S HD sim-ply reduces to S HD ≈ − N f (cid:88) (cid:126)x log det (1 + γP (cid:126)x ) (3.44)In particular for N c = 3 the determinant over the gaugegroup reduces todet (1 + γP (cid:126)x ) = 1 + γ Tr P (cid:126)x + γ Tr P − (cid:126)x + γ (3.45)Higher order corrections to S HD is given in powers of thehopping parameter, κ ≡ /m , and can be found in Ref.(Zambello and Di Renzo, 2018). Furthermore Ref. (Zam-bello and Di Renzo, 2018) focuses on µ ≈ µ c ≡ m = − log(2 κ ) where nuclear phase transition occurs at zero tem-perature. It is possible and convenient to work in the tem-poral gauge which eliminates all the links in P (cid:126)x but one forfixed (cid:126)x (say t = 0), so that P (cid:126)x = U ( (cid:126)x, t = 0) ≡ U (cid:126)x . Theholomorphic gradient flow equation (3.37) in the temporalgauge reads dU (cid:126)x dt = i (cid:80) a ( T a D a S HD [ U ]) U (cid:126)x = − γ (cid:80) a T a Tr( T a U (cid:126)x )det(1+ γU (cid:126)x ) U (cid:126)x . (3.46)The critical points satisfy Tr( T a U crx ) = 0 and therefore areelements of the center: U cr(cid:126)x = e iω (cid:126)x , ω x ∈ (cid:26) πnN c (cid:12)(cid:12)(cid:12)(cid:12) n = 0 , . . . , N c − (cid:27) . (3.47)Since ω x can take one of these three values at each lattice site,the number of critical points exponentially grows with thevolume as ( N c ) V . However they contribute to the path in-tegral with different weights. Ref. (Zambello and Di Renzo,2018) studied this model with N c = 3 in small spatial vol-umes up to 3 − and in a parameter range where only afew critical points, hence thimbles, contribute significantly tothe path integral and estimated their semiclassical weights, e − S HD [ U cr ] , by importance sampling. Furthermore they per-formed the Monte-Carlo computations of the charge density (cid:104) n (cid:105) and the Polyakov loop (cid:104) P (cid:105) = 1 /N (cid:80) (cid:126)x Tr (cid:104) U (cid:126)x (cid:105) over thethimbles with one and two lattice sites. The results showthe expected behavior in the cold limit near µ = µ cr , namely (cid:104) n (cid:105) sharply changing from 0 to 1 (i.e. the Silver Blaze be-havior (Cohen, 2003)) and (cid:104) P (cid:105) having a narrow peak around µ cr . Furthermore the contribution of three thimbles is nec-essary to obtain this expected result. Case study: 2D QED
Another example of a gauge theory, two-dimensional QEDwith the lattice action S = 1 g (cid:88) r (1 − cos P r ) − (cid:88) a ln det D ( a ) , (3.48)with D ( a ) xy = m a δ xy + 12 (cid:88) ν ∈{ , } (cid:2) η ν e iQ a A ν ( x )+ µδ ν δ x +ˆ ν,y − η ν e − iQ a A ν ( x ) − µδ ν δ x,y +ˆ ν (cid:3) . (3.49) The saturation of the fermion density is due to finite volume. et al. , 2018a) by using thegeneralized thimble method. Here D ( a ) denotes the (Kogut-Susskind) fermionic matrix for flavor a , and P r denotes theplaquette, P r ≡ A ( r ) + A ( r + ˆ x ) − A ( r + ˆ t ) − A ( r ) , (3.50)and ˆ t and ˆ x are the unit vectors in time and space direc-tion. In the case of abelian gauge theories it is convenientto work with the complexified gauge field, A µ ( x ) ∈ C N ,where N is the number of lattice sites, instead of the gaugelinks. This way, degeneracies to due gauge redundancy canbe addressed in a straightforward fashion. For any point x , the gauge orbit is generated by gauge transformations, A µ ( x ) → A µ ( x ) + α ( x + ˆ µ ) − α ( x ), which forms an N − subspace. The original real field space can there-fore be locally expressed as a direct product R N = M × G where M is the space of gauge inequivalent field configura-tions and G is the gauge orbit. The flow leaves G invariantbecause the gradient ∂S∂A is orthogonal to the gauge orbits,and it commutes with gauge transformations. Therefore themiddle-dimensional manifold (i.e. a manifold with real di-mension N ) obtained by flowing R N by an amount T , canbe decomposed as M T × G where M T is the result of flow-ing M by T . Furthermore, the critical points on the gaugefixed slice are now isolated and the thimble decompositionfollows straightforwardly as the limit T → ∞ . This argu-ment illustrates conceptually how the generalized thimblemethod works in the presence of an abelian gauge field. Forthe actual lattice computations, however, there is no needto fix the gauge as the Markov chain will randomly sample G which has no effect on the results as long as only gaugeinvariant observables are evaluated.In Ref. (Alexandru et al. , 2018a) the generalized thim-ble method has been put into action in a U (1) gauge theorywith three fermion species with charges 1 , −
2. Theircharges are chosen in such a way that a state with finitefermion density does not necessarily have a net charge, whichwould render the energy of the state infinite in the thermody-namics limit. It is shown that a computational speedup canbe achieved in computing the equation of state compared toconventional real space computation. Even though the coldlimit, which shows the Silver Blaze behavior (on a 14 ×
10 lat-tice), can be achieved this way, faster algorithms are neededto go to larger lattices. Finally, two-dimensional QED in In order to have neutral excitations a three flavor model with charges(2 , − , −
1) have been studied in Ref. (Alexandru et al. , 2018a). Thetwo flavor model with equal charges has no sign problem due to chargeconjugation symmetry. Note that α =constant is not a gauge transformation. the continuum has been studied in the mean field approxi-mation in (Tanizaki and Tachibana, 2017) where the exactthimble decomposition has been worked out for N f = 1 , , IV. OTHER MANIFOLDS AND THE ALGORITHMS THATCAN FIND THEMA. Well beyond thimbles
We have seen above that deforming the integration from R N to a proper combination of thimbles is not always de-sirable from the numerical point of view. The generalizedthimble method, for instance, uses a rough approximation ofthimbles that, while having a smaller average sign, has betterergodic properties. There is no reason, however, to be lim-ited to manifolds close to the thimbles. The condition that S I is constant is only one constraint in a 2 N dimensionalspace and, presumably, there are many manifolds of integra-tion where the sign of the integrand is fixed. In this sectionwe consider a few methods to search for other manifolds un-related to thimbles that both alleviate the sign problem andare numerically convenient.An ideal manifold of integration would i) ameliorate thesign problem significantly both because the action is nearlyreal on it and because the residual phase is small, ii) be com-putationally cheap to find and iii) be parametrized in sucha way that the associated Jacobian is also computationallycheap. All these restrictions are hard to satisfy at the sametime and only the first steps in this direction were taken. Itseems that insight into particular models will be essential toexploit this general idea profitably. We will show below that,in cases where the thimble method generate some of this in-sight, it is not difficult to improve it by allowing for moregeneral manifolds. In other theories it is an open problemto find a way to capitalize on the freedom of picking moregeneral manifolds. B. Learnifolds
Suppose a number of points on the thimble(s)—or someapproximation of it—are obtained using the computationallycostly holomorphic flow equations. It is reasonable to expectthat the sign fluctuations on a manifold that interpolatesbetween the original manifold, R N , and the thimble will besmall. A point found using the flow (in reality a complex fieldconfiguration or an element of C N ) can be viewed as a mapconnecting its real part to the corresponding imaginary part,that is for φ ∈ M T the map f takes Re φ → f (Re φ ) = Im φ .7Here we assume that the manifold does not “fold”; i.e. everyreal value of the field corresponds to an unique imaginarypart. We seek to find a manifold with a simple parametriza-tion that “interpolates” the configurations sampled on M T .The result of the interpolation problem can be thought as amap that approximates f , the map that connects the givenreal parts of coordinates to their imaginary parts. Thus,finding the interpolation of the points obtained by the flowcan be formulated as learning a general rule from a set ofexamples. This is a typical problem studied by the artifi-cial intelligence community and we can borrow some of theirtechniques to bear on it.More concretely, suppose we have a “training set”, S , thatis a number of field configurations φ ai lying on the mani-fold M T , where i = 1 , . . . , N indexes their components and a = 1 , . . . , |S| , with |S| being the size of the training set. Weparametrize the interpolating manifold L S , an approxima-tion to M T , as φ i = ζ i + i ˜ f i ( ζ ) , (4.1)where ζ i ∈ R N and ˜ f i is a real function approximating f .The function ˜ f i is represented by a feed-forward network ofthe type depicted in Fig. 17. The nodes on the left layerrepresent the input values, in our case the values of ζ i . Theresults are combined on the second layer by making linearcombinations of them, adding a bias and feeding them to anonlinear function σ ( x ) that we will take to be the of theform σ ( x ) = log(1 + e x ). The result is v j = σ (cid:32) b j + (cid:88) i w ij ζ i (cid:33) , (4.2)where j indexes the nodes of the second layer. These re-sults are then combined again, piped through σ ( x ) and fedto the next layer. At the end all results are combined in asingle number which represents ˜ f i =0 ( ζ ). By translation in-variance, the values of ˜ f i ( ζ ) for other i (cid:54) = 0 can be obtainedby translating the inputs ζ i . The feed-forward network isparametrized by the weights ( w ’s) of every link and biases b ’s of every node. These parameters are chosen in order tominimize the discrepancy between the training set and theresults of the network: C ( w, b ) = 1 |S| |S| (cid:88) a =1 (cid:12)(cid:12)(cid:12) ˜ f w,b (Re φ a ) − Im( φ a ) (cid:12)(cid:12)(cid:12) , (4.3)where ˜ f w,b (Re φ a ) is the result of applying the network, with ζ ζ ζ ζ f ( ζ ) FIG. 17: Topology of a feed-forward network with 5 layers:one input layer with 4 nodes, three intermediate layerswith 3 nodes each and one output layer with one node. Theinputs in the incoming layer (shown on the left) are the(real) values of the field. The output is the imaginary valueof the coordinate of the first point of the lattice ˜ f ( ζ ).parameters w ij and b j to Re φ a . In order to minimize C ( w, b ) a gradient descent algorithm is used. The com-putation of the gradient is efficiently done using the back-propagation algorithm. In fact, the existence of this simplealgorithm is an important motivation to use feed-forwardnetworks as opposed to networks with more complicatedtopologies. The minimization process is sped up tremen-dously by using the Adaptive Moment Estimate Algorithm(ADAM)(Kingma and Ba, 2014) (other methods are dis-cussed in (Ruder, 2016)), another borrow from the artificialintelligence literature. Since the manifold L S is defined bya network which learned how to approximate M T we call L S the “learnifold”. An example of the practical use of thismethod will be described in the next section.Three comments are worth making at this point. Firstly,while the usefulness of this method can be gauged on a caseby case basis, its correctness is guaranteed by construction.In fact, any network will define a manifold of integration inthe same homology class as R N since the mapping φ i ( ζ ) = ζ i + is ˜ f i ( ζ ) , (4.4)for 0 ≤ s ≤ R N and L S . Care must be taken on theasymptotic behavior of L S , determined the function ˜ f i ( ζ ), in Other cost functions can be used instead of C ( w, b ) used here. L S to be in the same class as M T but, in the caseof periodic ζ as in the applications below, this is automatic.Secondly, the parametrization Eq. 4.1 helps avoiding“trapping” of Monte Carlo chains as compared to theparametrization through the flow used in the (generalized)thimble method. This is explained by the Fig. 18b wherewe can see that the same manifold parametrized by a smallregion of R N (in the generalized thimble method) or a muchlarger region using Eq. 4.1. Regions of R N with large sta-tistical weights are then less separated by low probabilityregions, which facilitates the Monte Carlo sampling. Relat-edly, the Jacobian of the parametrization Eq. 4.1 fluctuatesless than the Jacobian of the flow parametrization.Finally, it is unlikely that the learnifold L S is better atcontrolling the sign problem than the manifold obtained byflow (of which the training set is taken). The usefulness ofthe method relies in the possibility of sampling the learni-fold at a cost orders of magnitude cheaper than flowing (andcomputing/estimating the Jacobian). The hope is that theincrease in statistics compensates allowed by the speed of theprocess compensates for the smaller average sign. Case study: 1+1D Thirring model revisited
The application of the learnifold method to the 1 + 1 di-mensional Thirring model is discussed in (Alexandru et al. ,2017d). The first step is to collect a number of complexconfigurations φ a = F T ( ζ a ) obtained by evolving real con-figurations ζ a by a “time” T according to the flow equationsEq. 2.7 in order to form the training set S . This step is iden-tical to what is done in the generalized thimble method. Onewould like to approximate the manifold M T particularly wellin the region that is going to be sampled the most. We alsowant to sample some of the field configurations on M T to-ward the large field value region to make sure that the profileof the learnifold matches M T in this region too. Thereforesome configurations are collected running a Metropolis chainwith weight e − S and some others with weight e − S/τ , with τ >
1. This way M T is sufficiently sampled so the interpo-lation manifold L S approximates it well in the statisticallyimportant regions and it is not radically wrong at asymp-totically far away regions. Notice that there is no particularreason for the Monte Carlo chain to be thermalized while col-lecting these configurations; al that is required is to have agood enough sampling of the statistically important regionsof M T and some sampling of the other regions. After theseconfigurations are generated they are used as the training setfor the feed-forward network. In order to enforce translationinvariance, all translations of them are added to the trainingset, resulting in a larger set, typically of the order of 10 el- (a) Parametrization using the holomorphic flow: small regions of R N map into large regions of thimbles (or M T ).(b) Parametrization of Eq. 4.1: regions of R N mapping into largeregions of L S are larger and with smaller gaps between them,which helps to prevent trapping of the Monte Carlo’s Markovchain. FIG. 18ements. This set is too large to be used in the minimizationprocess so different subsets (“minibatches”) are used at dif-ferent steps of the gradient descent (or ADAM step). Detailscan be found in (Alexandru et al. , 2017d). A comparison ofthe computational cost between the learnifold and the gen-eralized thimble methods is not straightforward because thecost of the learnifold method is divided into a “fixed” costrelated to the generation of the training set and the mini-mization of the cost function on one hand and the runningof the Monte Carlo given the optimal manifold. The sec-ond part is faster than the flowing required by the thimblemethod by orders of the magnitude but the first part maydominate the total costs. Calculations of the kind done in(Alexandru et al. , 2017b) can be done more effectively usingthe learnifold method.Similar methods have recently been applied to the solu-tion of the Hubbard model away from half-filling (Ulybyshev et al. , 2020b)9 T / m f ≃ T / m f ≃ T / m f ≃ T / m f ≃ T / m f ≃ μ 〈 ψψ 〉 FIG. 19: (cid:104) ¯ ψψ (cid:105) as a function of µ for the 2 + 1 dimensionalThirring model in a β × lattice showing the melting ofthe chiral condensate as the density is increased. The solidlines are fit to the functional form (cid:104) ¯ ψψ (cid:105) = A tanh( β ( µ − µ c )) C. Path optimization
An even more radical departure from thimbles is embod-ied in the path optimization method. The idea here is toconsider a family of manifolds M λ parametrized by a setof parameters (collectively denoted by λ ) and to maximizethe average sign within that family. The result of the max-imization defines a manifold M λ that is the best, withinthat family, at ameliorating the sign problem. The viabilityof the method rests on the observation that the gradient ofthe average sign in parameter space can be calculated witha (usually very short) sign-problem free Monte Carlo calcu-lation. Indeed, the average sign on a manifold M is givenby (cid:104) σ (cid:105) λ = (cid:82) M λ dφ e − S ( φ ) (cid:82) M λ | dφ | e − Re S ( φ ) = (cid:82) R N dζ e − S eff ( ζ ) (cid:82) R N dζ e − Re S eff ( ζ ) , (4.5)where S eff ( ζ ) = S [ φ ( ζ )] − ln det J ( ζ ) includes the determinantof the Jacobian of the M λ parametrization φ i = φ i ( ζ ). Thenumerator of Eq. 4.5 is independent of λ , due to Cauchy’stheorem. The denominator, however, being an integral of anon-holomorphic function, does depend on λ . In fact, ∇ λ (cid:104) σ (cid:105)(cid:104) σ (cid:105) = (cid:82) R N dζ e − S eff ( ζ ) (cid:0) −∇ λ Re S + Re Tr (cid:0) J − ∇ λ J (cid:1)(cid:1)(cid:82) R N dζ e − Re S eff ( ζ ) . (4.6)The average sign does not affect the direction of the vector ∇ λ (cid:104) σ (cid:105) and can be neglected during the maximization pro- cess while the right hand side term in Eq. 4.6, the average (cid:104)−∇ λ Re S +Re Tr (cid:0) J − ∇ λ J (cid:1) (cid:105) Re S eff , can be computed by theMonte Carlo method without encountering a sign problem.Knowledge of the gradient ∇ λ (cid:104) σ (cid:105) (to be more precise, knowl-edge of its direction in λ -space) allows a maximization rou-tine, like ADAM (Kingma and Ba, 2014) to find the values of λ leading to the largest possible sign within the parametrizedfamily of manifolds.A few facts make this scheme practical. First, a roughcomputation of the gradient is usually enough; some stochas-tic noise is actually useful in avoiding local minima that couldotherwise trap the maximization process. Second, the lastconfigurations obtained with one value of λ is close to beingthermalized as λ is changed to a nearby value during themaximization process, bypassing the need for long thermal-ization periods at each step of the minimization. Finally,it is imperative that the family of manifolds considered i)includes only manifolds in the same homology class as R N ,ii) contain manifolds where the sign problem is sufficientlyameliorated and iii) are parametrized in such a way the com-putation of the Jacobian J is cheap. Condition i) is relativelyeasy to satisfy but there is tension between conditions ii) andiii). Theoretical insight into the specific model of interest isrequired for the successful application of this method and itis currently sorely missed in most theories of physical signif-icance.We pause to note that contour deformations can be ap-plied to any theory with a holomorphic path integrand. Thisincludes theories with real actions but complex observables.This fact, combined with path optimization, has been used totame the signal to noise problem encountered in the calcula-tion of correlation function in simple field theories (Detmold et al. , 2020). Case study: 2+1D Thirring model
The path optimization method was applied to a one-dimensional integral in (Mori et al. , 2017), to the 1+1DThirring model at finite density in (Alexandru et al. , 2018b),to the 1+1D φ model in (Mori et al. , 2018) (with a neu-ral network parametrization of the manifolds similar to theone discussed in section IV.B ), the PNJL model in 0+1D(Kashiwa et al. , 2019a,b), 0+1D QCD (Mori et al. , 2019)and 1+1D φ (Bursa and Kroyter, 2018). Here we discussits application to the 3D Thirring model (Alexandru et al. ,2018c).The action defining the 2+1D Thirring model is in0FIG. 20: Phase diagram of the 3D Thirring model(Alexandru et al. , 2018c). The thick central line shows thelocation where (cid:104) ¯ ψψ (cid:105) µ,T = 0 . (cid:104) ¯ ψψ (cid:105) , and its width thestatistical errors. The thinner lines indicate (cid:104) ¯ ψψ (cid:105) µ,T = (0 . ± . (cid:104) ¯ ψψ (cid:105) , to help gauge the sharpnessof the transition.Eq. 3.15. The family of manifolds considered is given by A ( x ) = ζ ( x ) + i ( λ + λ cos ζ ( x ) + λ cos(2 ζ ( x ))) ,A ( x ) = ζ ( x ) ,A ( x ) = ζ ( x ) , (4.7)where λ , λ and λ are real numbers parametrizing the man-ifolds. The ansatz in Eq. 4.7 is motivated by the follow-ing considerations. Firstly, the determinant of the Jacobian J = ( ∂A ν /∂ζ i ) is trivial to compute (the cost scales withthe spacetime volume V , as opposed to V ) since the valueof A µ ( x ) depends only on ζ ( x ) evaluated at the same space-time point x :det J = (cid:89) x (1 − λ sin ζ ( x ) − λ sin(2 ζ ( x ))) . (4.8)Secondly, in the limit µ → ∞ the partition functionlim µ →∞ Z ≈ (cid:20)(cid:90) d A e g ( (cid:80) ν cos A ν )+ i A (cid:21) βV , (4.9)factorizes into a separate integral at every spacetime pointand the sign problem arises entirely from A . The ansatzin Eq. 4.7 reflects that. Thirdly, in the weak couplinglimit g → A = iα, A = A = 0, for some real constant α . The ansatz in Eq. 4.7contains manifolds that approach this thimble near its crit- ical point. Finally, the variables A ν are periodic variableswith a period 2 π (so they belong to ( S ) N , not R N ) andthe question of whether the manifolds defined by Eq. 4.7have a different asymptotic behavior is not present. Fur-thermore, by varying s from s = 1 to s = 0 in A ( x ) = ζ ( x ) + is ( λ + λ cos ζ ( x ) + λ cos(2 ζ ( x )) we see that everymember of the family of manifolds can be smoothly deformedto ( S ) N , guaranteeing the applicability of the Cauchy the-orem.This method was used in (Alexandru et al. , 2018c) in lat-tices of sizes up to 10 and action parameters near the con-tinuum. It is interesting to examine how the maximizationprocess proceeds. The parameter λ acquires very quicklya non-zero value very close to the position of the criticalpoint A = iα, A = A = 0. The corresponding manifolddoes not go through exactly through the critical point andhas a larger average sign than the space tangent to that crit-ical point. Afterwards, λ and λ settle on their preferredvalues, giving a little curvature to the manifold. More com-plicated functions of ζ in Eq. 4.7 do not seem to improvethe average sign. It seems that one is required to go beyondthe factorized form in Eq. 4.7 for further improvements.The results obtained by this method show a clear tran-sition (technically a crossover as the fermion mass breakschiral symmetry explicitly) between a phase with large chiralcondensate to another, at higher temperatures and densities,where chiral symmetry is restored and the chiral condensateis small (see Fig. 19). The resulting phase diagram is shownin Fig. 20. V. CONCLUSION AND PROSPECTS
A bird’s eye view of the developments described here re-veal some broad lessons that should not be lost amidst thetechnical details. The first is that the main idea the thim-ble approach to the sign problem is based on is sound andno fundamental flaw has been revealed, either conceptual orpractical. This is not to say one can currently use the methodto solve any sign problem. But the set of ideas exploredin this review provide a new setting where the simulationof many models can and should be attempted. This is nota trivial statement. There is a common perception amongnon-practitioners that there is a “conservation of difficulty”and that any approach to solve the sign problem will reveal,at closer inspection, a simulation cost as large as the naiveattempts. This is demonstrably untrue, as the examples dis-cussed in this review show.The second foundational lesson is that the integration overall the relevant thimbles is both necessary and possible, withseveral algorithms already proposed and tested, usually in1small scale simulations. In fact, the rapid algorithmic de-velopment in the last few years generated a problem – oran opportunity – as most calculations were aimed at demon-strating the algorithm correctness and scaling properties andnot focused on the Physics of the problem. For instance, cur-rent technology should be able to clarify the phase diagramof a variety of 1 + 1 dimensional models at finite tempera-ture/density. At a larger computational cost, 1 + 2 dimen-sional models can also be studied now. In fact, recent papersbegan setting up the path towards a solution to the repulsiveHubbard model away from half-filling, a result that wouldbe a game-changer in the field (Fukuma et al. , 2019, 2020;Hubbard, 1963; Mukherjee and Cristoforetti, 2014; Ostmeyer et al. , 2020; Saito, 2017; Ulybyshev et al. , 2019, 2020a; White et al. , 1989)An important insight arising from the research on thimble-related methods was that deformation of the integrationmanifold to other manifolds is both possible and profitable.This observation, as simple as it is, has vast consequences.Indeed, the condition that the imaginary part of the (effec-tive) action to be constant is only one constraint in a 2 N dimensional space. This leaves a 2 N − N − N − . A fewideas exist on how to explore this newfound freedom. One isto use information about the model obtained elsewhere to de-vise parametrized families of integration manifolds suitablefor that particular model. This approach provides a wayof bringing physical insight into a the Monte Carlo calcula-tion that is sometimes characterized as a brute force method.Whatever insight is brought to the model, obtained by rig-orous or intuitive, approximate methods, can then be usedto speed up a calculation, hopefully exponentially, that isguaranteed to converge to the correct answer by the MonteCarlo method. A surprising recent development is that thephysical insight into a model can be substituted by system-atic machine learning techniques. We expect the near futureto bring much more developments in this direction. ACKNOWLEDGMENTS
The authors would like to thank Hank Lamm, ScottLawrence, Greg Ridgway and Sohan Vartak, who collabo-rated with us on this topic. We also thank Tom Cohen, Notice that contrary to the multidimensional case, in the familiarcase of a single complex variable, the condition Im S eff = 0 defines aunique contour. Gerald Dunne, Tarun Grover, Sergei Gukov, Jay Sau, Chris-tian Schmidt, Luigi Scorzato, and Mithat ¨Unsal for conver-sations on this topic over the last few years. This work wassupported in part by the US DoE under contracts No. DE-FG02-93ER-40762, DE-FG02-95ER40907. AA gratefully ac-knowledges the hospitality of the Physics Department at theUniversity of Maryland where part of this work was carriedout.
Appendix A: Computation of the Jacobian
The evolution by the holomorphic flow Eq. 2.7 by time T maps initial conditions ζ into Φ( ζ ). The Jacobian of thistransformation is derived in this appendix.Begin by considering two infinitesimally close coordinates ζ and ζ (cid:48) , and let v v = φ ( ζ (cid:48) , − φ ( ζ,
0) (A1)denote the corresponding difference vector between them.Flowing for time ∆ t , both φ ( ζ (cid:48) ,
0) and φ ( ζ,
0) move, chang-ing the difference vector; we denote this time dependent dif-ference as v (∆ t ) (see Fig. 21). In the limit ∆ t → v a (∆ t ) ≡ φ a ( ζ (cid:48) , ∆ t ) − φ a ( ζ, ∆ t )= (cid:104) φ a ( ζ (cid:48) ,
0) + ∆ t ∂S∂φ a ( φ a ( ζ (cid:48) , (cid:105) − (cid:104) φ a ( ζ,
0) + ∆ t ∂S∂φ a ( φ a ( ζ, (cid:105) = (cid:2) φ a ( ζ (cid:48) , − φ a ( ζ, (cid:3) + ∆ t ∂ S∂φ a ∂φ b ( φ a ( ζ, (cid:2) φ b ( ζ (cid:48) , − φ b ( ζ, (cid:3) = v a (0) + ∆ tH ab ( φ a ( ζ, v b (0) . (A2)In other words, a vector evolves along a flow trajectory ac-cording to the differential equation dv a dt = H ab ( φ ( ζ, t )) v b ( t ) . (A3)We can use the equation above to evolve a set of N vectorsforming an orthonormal basis. Packaging these vector in thecolumns of a matrix J (0) =
1, we see that J ( t ) obeys dJdt = HJ. (A4)2FIG. 21: Two nearby points evolving by the holomorphicflow. Their difference vector, shown in blue, evolvesaccording to Eq. A3.
Appendix B: Another definition for thimbles
In this appendix we give a different perspective on Lef-schetz thimbles. We begin with focusing on the stationarypoints of the flow, namely the critical points of the action φ c , such that ∂S/∂φ i | φ c = 0. Around a critical point it isalways possible to find local coordinates { z i = x i + iy i } with i = 1 , . . . , N such that S ( φ ) − S ( φ c ) = z + · · · + z N = ( x + · · · + x N ) − ( y + · · · + y N )+2 i ( x y + · · · + x N y N ) (B1)whose existence is guaranteed by the Morse lemma. Nowconsider the N − v ( s ) defined by x + · · · + x N = s, y = · · · = y N = 0. This surface is known as the vanishing cycle as it vanishes at the critical point. It can beviewed as the level set of the action around the critical point, S − ( s + s c ) where s c = S ( φ c ). We now “move” the vanishingcycle by varying s . This can be done by taking the vanishingcycle around the critical point v ( (cid:15) ) and then flowing it. As s runs from 0 to ∞ , the vanishing cycle sweeps an N (real) In our analysis we consider only isolated, quadratic (non-degenerate)critical points. A degenerate critical point where the Hessian deter-minant of φ vanishes can be split into µ number of non-degeneratecritical points with a small deformation with µ being the Milnornumber of the critical point and a similar analysis presented in thissection follows (Pham, 1983). dimensional surface. This N dimensional surface, defined asthe union of vanishing cycles on the half line 0 ≤ s < ∞ , T = ∪ s v ( s ), is known as the Lefschetz thimble associated with thecritical point φ c . Similarly, we define an N − v D ( s ) by x = · · · = x N = 0 , y + · · · + y N = s .We shall call the union of these dual cycles on the half line0 ≤ s < ∞ , K = ∪ s v D ( s ), the dual thimble . REFERENCES
Aarts, G. (2009), Phys. Rev. Lett. , 131601, arXiv:0810.2089[hep-lat].Aarts, G. (2016), Journal of Physics: Conference Series ,022004.Aarts, G., F. Attanasio, B. J¨ager, and D. Sexty (2016), JHEP , 087, arXiv:1606.05561 [hep-lat].Aarts, G., L. Bongiovanni, E. Seiler, D. Sexty, and I.-O. Sta-matescu (2013), Eur. Phys. J. A49 , 89, arXiv:1303.6425 [hep-lat].Aarts, G., F. A. James, E. Seiler, and I.-O. Stamatescu (2011),Eur. Phys. J.
C71 , 1756, arXiv:1101.3270 [hep-lat].Aarts, G., E. Seiler, and I.-O. Stamatescu (2010), Phys. Rev.
D81 , 054508, arXiv:0912.3360 [hep-lat].Aarts, G., and I.-O. Stamatescu (2008), JHEP , 018,arXiv:0807.1597 [hep-lat].Alexandru, A., G. Basar, and P. Bedaque (2015),arXiv:1510.03258 [hep-lat].Alexandru, A., G. Basar, P. Bedaque, G. W. Ridgway, andN. C. Warrington (2016a), Phys. Rev. D94 (4), 045017,arXiv:1606.02742 [hep-lat].Alexandru, A., G. Ba¸sar, P. F. Bedaque, H. Lamm,and S. Lawrence (2018a), Phys. Rev.
D98 (3), 034506,arXiv:1807.02027 [hep-lat].Alexandru, A., G. Basar, P. F. Bedaque, and G. W. Ridgway(2017a), Phys. Rev.
D95 (11), 114501, arXiv:1704.06404 [hep-lat].Alexandru, A., G. Basar, P. F. Bedaque, G. W. Ridgway,and N. C. Warrington (2016b), Phys. Rev.
D93 (9), 094514,arXiv:1604.00956 [hep-lat].Alexandru, A., G. Basar, P. F. Bedaque, G. W. Ridgway,and N. C. Warrington (2017b), Phys. Rev.
D95 (1), 014502,arXiv:1609.01730 [hep-lat].Alexandru, A., G. Basar, P. F. Bedaque, S. Vartak, andN. C. Warrington (2016c), Phys. Rev. Lett. (8), 081602,arXiv:1605.08040 [hep-lat].Alexandru, A., G. Basar, P. F. Bedaque, and N. C. Warrington(2017c), arXiv:1703.02414 [hep-lat].Alexandru, A., P. F. Bedaque, H. Lamm, and S. Lawrence(2017d), Phys. Rev.
D96 (9), 094505, arXiv:1709.01971 [hep-lat]. Note that what we call the thimble and the dual thimble are re-ferred as downward/upward cycles referring to the fact that theweight, e − Re S , monotonically decreases/increases over them withflow(Witten, 2011). Alexandru, A., P. F. Bedaque, H. Lamm, and S. Lawrence(2018b), arXiv:1804.00697 [hep-lat].Alexandru, A., P. F. Bedaque, H. Lamm, S. Lawrence, andN. C. Warrington (2018c), Phys. Rev. Lett. (19), 191602,arXiv:1808.09799 [hep-lat].Alexandru, A., G. Bergner, D. Schaich, and U. Wenger (2018d),Phys. Rev.
D97 (11), 114503, arXiv:1712.07585 [hep-lat].Alexandru, A., M. Faber, I. Horvath, and K.-F. Liu (2005), Phys.Rev.
D72 , 114513, arXiv:hep-lat/0507020 [hep-lat].Alexandru, A., and U. Wenger (2011), Phys. Rev.
D83 , 034502,arXiv:1009.2197 [hep-lat].Alford, M. G., S. Chandrasekharan, J. Cox, and U. J. Wiese(2001), Nucl. Phys.
B602 , 61, arXiv:hep-lat/0101012 [hep-lat].Anderson, J. B. (1975), The Journal of Chemical Physics (4),1499, https://doi.org/10.1063/1.431514.Aniceto, I., G. Basar, and R. Schiappa (2019), Phys. Rept. ,1, arXiv:1802.10441 [hep-th].Attanasio, F., and J. E. Drut (2020), Phys. Rev. A , 033617.Ayyar, V., S. Chandrasekharan, and J. Rantaharju (2018), Phys.Rev. D (5), 054501, arXiv:1711.07898 [hep-lat].Barbour, I. M., C. T. H. Davies, and Z. Sabeur (1988), Phys.Lett. B215 , 567.Bazavov, A., et al. (HotQCD) (2019), Phys. Lett.
B795 , 15,arXiv:1812.08235 [hep-lat].Bellwied, R., S. Borsanyi, Z. Fodor, J. G¨unther, S. D. Katz,C. Ratti, and K. K. Szabo (2015), Phys. Lett.
B751 , 559,arXiv:1507.07510 [hep-lat].Bender, I., T. Hashimoto, F. Karsch, V. Linke, A. Nakamura,M. Plewnia, I. Stamatescu, and W. Wetzel (1992), Nucl. Phys.B Proc. Suppl. , 323.Berger, C. E., K. J. Morrell, and J. E. Drut (2020), “Thermo-dynamics of rotating quantum matter in the virial expansion,”arXiv:2004.02344 [cond-mat.quant-gas].Berger, C. E., L. Rammelm¨uller, A. C. Loheac, F. Ehmann,J. Braun, and J. E. Drut (2019), arXiv:1907.10183 [cond-mat.quant-gas].Berges, J., S. Borsanyi, D. Sexty, and I. O. Stamatescu (2007),Phys. Rev. D75 , 045007, arXiv:hep-lat/0609058 [hep-lat].Berkowitz, E., E. Rinaldi, M. Hanada, G. Ishiki, S. Shi-masaki, and P. Vranas (2016), Phys. Rev. D (9), 094501,arXiv:1606.04951 [hep-lat].Bloch, J., F. Bruckmann, and T. Wettig (2013), JHEP , 140,arXiv:1307.1416 [hep-lat].Blum, T. C., J. E. Hetrick, and D. Toussaint (1996), Phys. Rev.Lett. , 1019, arXiv:hep-lat/9509002.Bonati, C., M. D’Elia, M. Mariti, M. Mesiti, F. Negro,and F. Sanfilippo (2015), Phys. Rev. D92 (5), 054503,arXiv:1507.03571 [hep-lat].Bonati, C., M. D’Elia, F. Negro, F. Sanfilippo, and K. Zambello(2018), Phys. Rev.
D98 (5), 054510, arXiv:1805.02960 [hep-lat].Borsanyi, S., Z. Fodor, J. N. Guenther, R. Kara, S. D. Katz,P. Parotto, A. Pasztor, C. Ratti, and K. K. Szabo (2020),arXiv:2002.02821 [hep-lat].Braun, J., J.-W. Chen, J. Deng, J. E. Drut, B. Friman, C.-T. Ma,and Y.-D. Tsai (2013), Phys. Rev. Lett. , 130404.Bursa, F., and M. Kroyter (2018), JHEP , 054,arXiv:1805.04941 [hep-lat]. Carlson, J., S. Gandolfi, F. Pederiva, S. C. Pieper, R. Schiav-illa, K. Schmidt, and R. Wiringa (2015), Reviews of ModernPhysics (3), 1067.Cea, P., L. Cosmai, and A. Papa (2014), Phys. Rev. D89 (7),074512, arXiv:1403.0821 [hep-lat].Chandrasekharan, S. (2012), Phys. Rev. D , 021701,arXiv:1205.0084 [hep-lat].Chandrasekharan, S. (2013), Eur. Phys. J. A49 , 90,arXiv:1304.4900 [hep-lat].Chandrasekharan, S., and U.-J. Wiese (1999), Phys. Rev. Lett. , 3116, arXiv:cond-mat/9902128 [cond-mat].Chen, J.-W., and D. B. Kaplan (2004), Phys. Rev. Lett. ,257002.Cherman, A., D. Dorigoni, and M. Unsal (2014), arXiv:1403.1277[hep-th].Cohen, T. D. (2003), Phys. Rev. Lett. , 222001, arXiv:hep-ph/0307089 [hep-ph].Cristoforetti, M., F. Di Renzo, G. Eruzzi, A. Mukherjee,C. Schmidt, L. Scorzato, and C. Torrero (2014), Phys. Rev. D89 (11), 114505, arXiv:1403.5637 [hep-lat].Cristoforetti, M., F. Di Renzo, A. Mukherjee, and L. Scorzato(2013), Phys. Rev.
D88 (5), 051501, arXiv:1303.7204 [hep-lat].Cristoforetti, M., F. Di Renzo, and L. Scorzato (AuroraScience)(2012a), Phys. Rev.
D86 , 074506, arXiv:1205.3996 [hep-lat].Cristoforetti, M., L. Scorzato, and F. Di Renzo (2012b),
Proceed-ings, Extreme QCD 2012 (XQCD12): Washington, USA, Au-gust 21-23, 2012 , 10.1088/1742-6596/432/1/012025, [J. Phys.Conf. Ser.432,012025(2013)], arXiv:1210.8026 [hep-lat].Delabaere, E., and C. J. Howls (2002), Duke Math. J. (2),199.D’Elia, M., and M.-P. Lombardo (2003), Phys. Rev.
D67 , 014505,arXiv:hep-lat/0209146 [hep-lat].D’Elia, M., and M. P. Lombardo (2004), Phys. Rev. D , 074509,arXiv:hep-lat/0406012.Detmold, W., G. Kanwar, M. L. Wagman, and N. C. Warrington(2020), “Path integral contour deformations for noisy observ-ables,” arXiv:2003.05914 [hep-lat].Di Renzo, F., and G. Eruzzi (2018), Phys. Rev. D97 (1), 014503,arXiv:1709.10468 [hep-lat].Di Renzo, F., S. Singh, and K. Zambello (2019), in , arXiv:2002.00472[hep-lat].Duane, S., A. Kennedy, B. Pendleton, and D. Roweth (1987),Phys.Lett.
B195 , 216.Dunne, G. V., and M. ¨Unsal (2016),
Proceedings, 33rd Inter-national Symposium on Lattice Field Theory (Lattice 2015):Kobe, Japan, July 14-18, 2015 , PoS
LATTICE2015 , 010,arXiv:1511.05977 [hep-lat].Earl, D. J., and M. W. Deem (2005), Physical Chemistry Chem-ical Physics (23), 3910.Elhatisari, S., E. Epelbaum, H. Krebs, T. A. L¨ahde, D. Lee, N. Li,B.-n. Lu, U.-G. Meißner, and G. Rupak (2017), Phys. Rev.Lett. , 222505.Endres, M. G. (2007), Phys. Rev. D75 , 065012, arXiv:hep-lat/0610029 [hep-lat].Endrodi, G., Z. Fodor, S. D. Katz, and K. K. Szabo (2011), JHEP , 001, arXiv:1102.1356 [hep-lat]. Epelbaum, E., H. Krebs, T. A. L¨ahde, D. Lee, U.-G. Meißner,and G. Rupak (2014), Physical Review Letters (10),10.1103/physrevlett.112.102501.Fantoni, S., A. Sarsa, and K. E. Schmidt (2001), Physical ReviewLetters (18), 10.1103/physrevlett.87.181101.Fodor, Z., and S. D. Katz (2002), JHEP , 014, arXiv:hep-lat/0106002 [hep-lat].Fodor, Z., S. D. Katz, and C. Schmidt (2007), JHEP , 121,arXiv:hep-lat/0701022 [hep-lat].Fodor, Z., S. D. Katz, D. Sexty, and C. T¨or¨ok (2015),arXiv:1508.05260 [hep-lat].de Forcrand, P. (2010), “Simulating qcd at finite density,”arXiv:1005.0539 [hep-lat].de Forcrand, P., and S. Kratochvila (2006), Hadron physics, pro-ceedings of the Workshop on Computational Hadron Physics,University of Cyprus, Nicosia, Cyprus, 14-17 September 2005 ,Nucl. Phys. Proc. Suppl. , 62, arXiv:hep-lat/0602024 [hep-lat].de Forcrand, P., and O. Philipsen (2002), Nucl. Phys.
B642 , 290,arXiv:hep-lat/0205016 [hep-lat].de Forcrand, P., and O. Philipsen (2003), Nucl. Phys.
B673 , 170,arXiv:hep-lat/0307020 [hep-lat].de Forcrand, P., et al. (QCD-TARO) (2000),
Lattice field theory.Proceedings, 17th International Symposium, Lattice’99, Pisa,Italy, June 29-July 3, 1999 , Nucl. Phys. Proc. Suppl. , 408,arXiv:hep-lat/9911034 [hep-lat].Fujii, H., D. Honda, M. Kato, Y. Kikukawa, S. Komatsu, andT. Sano (2013), JHEP , 147, arXiv:1309.4371 [hep-lat].Fujii, H., S. Kamata, and Y. Kikukawa (2015), JHEP , 078,[Erratum: JHEP02,036(2016)], arXiv:1509.08176 [hep-lat].Fujii, H., S. Kamata, and Y. Kikukawa (2017), arXiv:1710.08524[hep-lat].Fukuma, M., N. Matsumoto, and N. Umeda (2018), JHEP ,060, arXiv:1806.10915 [hep-th].Fukuma, M., N. Matsumoto, and N. Umeda (2019), Phys. Rev.D , 114510.Fukuma, M., N. Matsumoto, and N. Umeda (2020), in , arXiv:2001.01665[hep-lat].Fukuma, M., and N. Umeda (2017), arXiv:1703.00861 [hep-lat].Gandolfi, S., A. Lovato, J. Carlson, and K. E. Schmidt (2014),Physical Review C (6), 10.1103/physrevc.90.061306.Garron, N., and K. Langfeld (2016), Eur. Phys. J. C76 (10), 569,arXiv:1605.02709 [hep-lat].Garron, N., and K. Langfeld (2017), Eur. Phys. J.
C77 (7), 470,arXiv:1703.04649 [hep-lat].Gattringer, C., and T. Kloiber (2013), Nucl. Phys.
B869 , 56,arXiv:1206.2954 [hep-lat].Gattringer, C., and P. T¨orek (2015), Phys. Lett.
B747 , 545,arXiv:1503.04947 [hep-lat].Geyer, C. J. (1991), Interface Foundation of North America.Gezerlis, A. (2011), Phys. Rev. C , 065801.Grady, M. (1985), Phys. Rev. D32 , 1496.Hanada, M., J. Nishimura, Y. Sekino, and T. Yoneya (2011),JHEP , 020, arXiv:1108.5153 [hep-th].Hann, C. T., E. Huffman, and S. Chandrasekharan (2017), An-nals Phys. , 63, arXiv:1608.05144 [cond-mat.str-el]. Harlow, D., J. Maltz, and E. Witten (2011), JHEP , 071,arXiv:1108.4417 [hep-th].Hasenfratz, A., and D. Toussaint (1992), Nucl. Phys. B371 , 539.Hubbard, J. (1963), Proceedings of the Royal Society of London.Series A, Mathematical and Physical Sciences (1365), 238.Huffman, E., and S. Chandrasekharan (2016), Phys. Rev. E (4), 043311, arXiv:1605.07420 [cond-mat.str-el].Huffman, E., and S. Chandrasekharan (2020), Phys. Rev. D (7), 074501, arXiv:1912.12823 [cond-mat.str-el].Huffman, E. F., and S. Chandrasekharan (2014), Phys. Rev. B (11), 111101, arXiv:1311.0034 [cond-mat.str-el].Kaczmarek, O., F. Karsch, E. Laermann, C. Miao, S. Mukherjee,P. Petreczky, C. Schmidt, W. Soeldner, and W. Unger (2011),Phys. Rev. D83 , 014504, arXiv:1011.3130 [hep-lat].Karsch, F. (2000), Nuclear Physics B - Proceedings Supplements , 14 , proceedings of the XVIIth International Symposiumon Lattice Field Theory.Karsch, F., and K. H. Mutter (1989), Nucl. Phys.
B313 , 541.Kashiwa, K., Y. Mori, and A. Ohnishi (2019a), Phys. Rev.
D99 (11), 114005, arXiv:1903.03679 [hep-lat].Kashiwa, K., Y. Mori, and A. Ohnishi (2019b), Phys. Rev.
D99 (1), 014033, arXiv:1805.08940 [hep-ph].Keldysh, L. (1964), Zh. Eksp. Teor. Fiz. , 1515.Kingma, D. P., and J. Ba (2014), ArXiv e-prints arXiv:1412.6980[cs.LG].Klauder, J. R. (1983), Recent developments in high-energyphysics. Proceedings, 22. Internationale Universit¨atswochen f¨urKernphysik: Schladming, Austria, February 23 - March 5,1983 , Acta Phys. Austriaca Suppl. , 251.Kratochvila, S., and P. de Forcrand (2006), Proceedings,23rd International Symposium on Lattice field theory (Lattice2005): Dublin, Ireland, Jul 25-30, 2005 , PoS
LAT2005 , 167,arXiv:hep-lat/0509143 [hep-lat].Lacroix, C., P. Mendels, and F. Mila (2011), Introduction to Frus-trated Magnetism: Materials, Experiments, Theory, SpringerSeries in Solid-State Sciences, Volume 164. ISBN 978-3-642-10588-3. Springer-Verlag Berlin Heidelberg, 2011 -1 .L¨ahde, T. A., T. Luu, D. Lee, U.-G. Meißner, E. Epelbaum,H. Krebs, and G. Rupak (2015), “Nuclear lattice simulationsusing symmetry-sign extrapolation,” arXiv:1502.06787 [nucl-th].Langelage, J., M. Neuman, and O. Philipsen (2014), JHEP ,131, arXiv:1403.4162 [hep-lat].Langfeld, K., and B. Lucini (2014), Phys. Rev. D90 (9), 094502,arXiv:1404.7187 [hep-lat].Lawrence, S. (2020),
Sign Problems in Quantum Field The-ory: Classical and Quantum Approaches , Other thesisarXiv:2006.03683 [hep-lat].Lee, D. (2007), Physical Review Letters (18), 10.1103/phys-revlett.98.182501.Lee, D. (2009), Progress in Particle and Nuclear Physics (1),117.Lee, D., B. Borasoy, and T. Schaefer (2004), Physical Review C (1), 10.1103/physrevc.70.014007.Lee, D., and T. Sch¨afer (2005), Physical Review C (2),10.1103/physrevc.72.024006.Li, A., A. Alexandru, and K.-F. Liu (2011), Phys. Rev. D84 ,071503, arXiv:1103.3045 [hep-ph]. Li, A., A. Alexandru, K.-F. Liu, and X. Meng (2010), Phys. Rev.
D82 , 054502, arXiv:1005.4158 [hep-lat].Li, D. (2016), arXiv:1605.04623 [hep-lat].Loh, E. Y., J. E. Gubernatis, R. T. Scalettar, S. R. White, D. J.Scalapino, and R. L. Sugar (1990), Phys. Rev. B , 9301.Lu, B.-N., N. Li, S. Elhatisari, D. Lee, J. E. Drut, T. A. L¨ahde,E. Epelbaum, and U.-G. Meißner (2019a), “Ab initio nuclearthermodynamics,” arXiv:1912.05105 [nucl-th].Lu, B.-N., N. Li, S. Elhatisari, D. Lee, E. Epelbaum, and U.-G.Meißner (2019b), Physics Letters B , 134863.Metropolis, N., A. W. Rosenbluth, M. N. Rosenbluth, A. H.Teller, and E. Teller (1953), The Journal of Chemical Physics (6), 1087.Miyamura, O. (QCD-TARO) (2002), Quark matter 2001. Proceed-ings, 15th International Conference on Ultrarelativistic nucleusnucleus collisions, QM 2001, Stony Brook, USA, January 15-20, 2001 , Nucl. Phys.
A698 (1-4), 395.Mori, Y., K. Kashiwa, and A. Ohnishi (2017), Phys. Rev.
D96 (11), 111501, arXiv:1705.05605 [hep-lat].Mori, Y., K. Kashiwa, and A. Ohnishi (2018), PTEP (2),023B04, arXiv:1709.03208 [hep-lat].Mori, Y., K. Kashiwa, and A. Ohnishi (2019), PTEP (11),113B01, arXiv:1904.11140 [hep-lat].Mou, Z.-G., P. M. Saffin, and A. Tranberg (2019a), JHEP ,135, arXiv:1909.02488 [hep-th].Mou, Z.-G., P. M. Saffin, A. Tranberg, and S. Woodward (2019b),JHEP , 094, arXiv:1902.09147 [hep-lat].Mukherjee, A., and M. Cristoforetti (2014), Phys. Rev. B90 (3),035134, arXiv:1403.5680 [cond-mat.str-el].Mukherjee, A., M. Cristoforetti, and L. Scorzato (2013), Phys.Rev.
D88 (5), 051502, arXiv:1308.0233 [physics.comp-ph].Muroya, S., A. Nakamura, C. Nonaka, andT. Takaishi (2003), Progress of Theoretical Physics (4), 615, https://academic.oup.com/ptp/article-pdf/110/4/615/5301553/110-4-615.pdf.Neal, R. (1996), Statistics and Computing (4), 353.Ostmeyer, J., E. Berkowitz, S. Krieg, T. A. L¨ahde, T. Luu,and C. Urbach (2020), “The semimetal-mott insulator quan-tum phase transition of the hubbard model on the honeycomblattice,” arXiv:2005.11112 [cond-mat.str-el].Parisi, G. (1983), Phys. Lett. , 393.Parisi, G., and Y.-s. Wu (1981), Sci. Sin. , 483.Pawlowski, J. M., M. Scherzer, C. Schmidt, F. P. G. Ziegler, andF. Ziesch´e (2020), in , arXiv:2001.09767 [hep-lat].Pawlowski, J. M., and C. Zielinski (2013), Phys. Rev. D (9),094503, arXiv:1302.1622 [hep-lat].Pham, F. (1983), Singularities (Arcata, Calif., 1981) , MethodsAppl. Anal. 1
Part 2 .Philipsen, O. (2007), The European Physical Journal Special Top-ics (1), 29.Roscher, D., J. Braun, and J. E. Drut (2014), Phys. Rev. A ,063609. Rossi, P., and U. Wolff (1984), Nucl. Phys. B248 , 105.Ruder, S. (2016), ArXiv e-prints arXiv:1609.04747 [cs.LG].Saito, H. (2017), Journal of the Physical Society of Japan (9),093001, https://doi.org/10.7566/JPSJ.86.093001.Schaich, D. (2019), PoS LATTICE2018 , 005, arXiv:1810.09282[hep-lat].Schwinger, J. S. (1961), J. Math. Phys. , 407.Seiler, E., D. Sexty, and I.-O. Stamatescu (2013), Phys. Lett. B723 , 213, arXiv:1211.3709 [hep-lat].Sexty, D. (2014), Phys. Lett.
B729 , 108, arXiv:1307.7748 [hep-lat].Sindzingre, P., P. Lecheminant, and C. Lhuillier (1994), Phys.Rev. B , 3108.Swendsen, R. H., and J.-S. Wang (1986), Phys. Rev. Lett. ,2607.Takagi, T. (1924), Japanese journal of mathematics :transactionsand abstracts , 83.Tanizaki, Y., Y. Hidaka, and T. Hayata (2016), New J. Phys. (3), 033002, arXiv:1509.07146 [hep-th].Tanizaki, Y., H. Nishimura, and J. J. M. Verbaarschot (2017),arXiv:1706.03822 [hep-lat].Tanizaki, Y., and M. Tachibana (2017), JHEP , 081,arXiv:1612.06529 [hep-th].Thirring, W. E. (1958), Annals Phys. , 91.Troyer, M., and U.-J. Wiese (2005), Phys. Rev. Lett. , 170201,arXiv:cond-mat/0408370 [cond-mat].Ulybyshev, M., C. Winterowd, and S. Zafeiropoulos (2019),arXiv:1906.02726 [cond-mat.str-el].Ulybyshev, M., C. Winterowd, and S. Zafeiropoulos (2020a),Phys. Rev. D101 (1), 014508, arXiv:1906.07678 [cond-mat.str-el].Ulybyshev, M. V., V. I. Dorozhinskii, and O. V. Pavlovskii(2020b), Physics of Particles and Nuclei (3), 363.White, S. R., D. J. Scalapino, R. L. Sugar, E. Y. Loh, J. E.Gubernatis, and R. T. Scalettar (1989), Phys. Rev. B , 506.Wigner, E. (1937), Phys. Rev. , 106.Wiringa, R. B., S. C. Pieper, J. Carlson, and V. R. Pand-haripande (2000), Physical Review C (1), 10.1103/phys-revc.62.014001.Witten, E. (1978), Nuclear Physics B (1), 110.Witten, E. (2010), “A new look at the path integral of quantummechanics,” arXiv:1009.6032 [hep-th].Witten, E. (2011), Chern-Simons gauge theory: 20 years after.Proceedings, Workshop, Bonn, Germany, August 3-7, 2009 ,AMS/IP Stud. Adv. Math. , 347, arXiv:1001.2933 [hep-th].Zambello, K., and F. Di Renzo (2018), Proceedings, 36th In-ternational Symposium on Lattice Field Theory (Lattice 2018):East Lansing, MI, United States, July 22-28, 2018 , PoS
LAT-TICE2018 , 148, arXiv:1811.03605 [hep-lat].Zhang, S., J. Carlson, and J. E. Gubernatis (1995), Phys. Rev.Lett. , 3652.Zhang, S., J. Carlson, and J. E. Gubernatis (1997), Phys. Rev.B55